A Numerical Solution of Heat Equation for Several Thermal Diffusivity Using Finite Difference Scheme with Stability Conditions

Wahida Zaman Loskor^1*, Rama Sarkar²

¹Department of Science and Humanities, Bangladesh Army International University of Science and Technology, Cumilla, Bangladesh.
²Department of Mathematics and Statistics, Bangladesh University of Business and Technology, Dhaka, Bangladesh.

DOI: 10.4236/jamp.2022.102034                

Abstract

The heat equation is a second-order parabolic partial differential equation, which can be solved in many ways using numerical methods. This paper provides a numerical solution that uses the finite difference method like the explicit center difference method. The forward time and centered space (FTCS) is used to a problem containing the one-dimensional heat equation and the stability condition of the scheme is reported with different thermal conductivity of different materials. In this study, results obtained for different thermal conductivity of distinct materials are compared. Also, the results reveal the well-behavior properties of the materials in good agreement.

Keywords

Heat Equation, Finite-Difference Scheme, Explicit Centered Difference Scheme, Thermal Diffusivity

1. Introduction

Heat propagation in a continuous medium, such as a solid, liquid, or gas, is described by the heat equation, which is a partial differential equation. It is the basic mathematical theory of thermal conductivity. In other words, the heat equation narrates the distribution of heat (or temperature fluctuation) in a particular area over time.

The one-dimensional heat equation solution is crucial since it appears often in numerous scientific and engineering applications. Also, the matter of simulation techniques of the heat equation has become an interesting area within the sector of numerical solution methods. Many scientific groups are involved in handling the matter of numerical simulation of heat flow [1]. [2] provides a numerical solution for transient thermal distribution in a portion where chemical, electrical, or nuclear energy is transformed into thermal energy. [3] used an exponential-sinusoidal one-dimensional analytical model to define the temperature property, demonstrating that the heat equation can still be solved analytically. [4] solved one-dimensional transient heat transport during a compound slab using a natural analytical technique. He investigated the transient response of one-dimensional multilayered compound conduction slabs to rapid temperature changes in the surrounding fluid. [1] provided a completely unique analytical approach utilized to solve the transient heat conduction equation in a one-dimensional hollow compound cylinder with time-dependent boundary conditions. In [5], the stability analysis of a one-dimensional heat equation is presented. They investigated the analytical solution of heat equation as an initial value problem in infinite space. It’s found that centered difference schemes (CNS) are more efficient and accurate solutions than explicit centered difference schemes (ECDS) in terms of time step selection for this equation. In [6], they build and validate mathematical models and numerical algorithms for the simulation of heat transfer in composite materials. The problem’s key characteristics are the heat source’s dependency on the solution, discontinuous diffusion coefficients, and nonlinear convection and radiation boundary conditions. This [7] study proposes a novel technique for solving diffusion and heat equations that combines the variation iterative method with an integral transform similar to the Sumudu transform. The method is accurate and efficient in the generation of approximate solutions for partial differential equations, it says. This [8] paper proposes a numerical technique to obtain the solution for the heat conduction equation of Copper. It says Copper is highly able to conduct heat and electrical conductivity. The results indicate that the finite difference method is highly effective for obtaining approximate solutions of thermal conductivity equations for Copper.

The flowing of heat through a thin rod formed of a homogenous material that is fully heated towards its length so that heat can only flow with its edges is considered in this work. Any position towards the rod is symbolized by x, and the length of the rod is symbolized by L in order that $0\le x\le L$. Therefore the temperature, $u\left(x,t\right)$ of the rod at any point is a function of position, x (in meters or centimeters) and time t (in second). The heat equation is also named expansion equation.

The present study involves the implementation of a finite difference scheme for heat equation for different thermal diffusivities with stability conditions and performs a computer programming for the method to compare the efficiency and perfection of the solution for seven different thermal diffusivities for seven different materials such as Silver, Gold, Copper, Aluminum, Cast Iron, Granite and Brick materials.

1.1. Mathematical Model of Heat Equation and Analytic Solution

The heat equation is resulting from the principles of conservation of energy and Fourier’s law of cooling. The heat flow is used to identify the equation in terms of temperature measurements. The one-dimensional heat equation is a partial differential equation of second order that describes how heat moves through an object over time and written as [9],

$\frac{\partial u}{\partial t}=k{\nabla }^{2}u$ (1)

where $k=\frac{K}{\rho c}$ is called the thermal diffusivity of the rod’s substance.

K = thermal conductivity of the rod’s material;

$\rho$ = mass per unit volume of the rod;

c = the rod’s specific heat.

To find the analytic solution of the one-dimensional heat equation appended by initial and boundary conditions that formulate an initial boundary value problem (IBVP) is considered. As a result, with beginning and boundary conditions, the heat equation becomes:

$\begin{array}{l}{u}_t\left(x,t\right)=k{u}_{xx}\left(x,t\right)\text{for}0\le x\le L\text{and}t>0\\ u\left(x,0\right)=f\left(x\right)\text{for}0\le x\le L\\ u\left(0,t\right)=u\left(L,t\right)=0\text{for}t>0\end{array}\}$ (2)

The exact solution can be found by applying separation of variables which is described in the following.

Consider $u\left(x,t\right)=X\left(x\right)T\left(t\right)$ then the heat equation becomes: $X\left(x\right)T^{\prime }\left(x\right)=k{X}^{″}\left(x\right)T\left(t\right)$.

Divided by $kT\left(t\right)$ and $X\left(x\right)$, it implies, $\frac{T^{\prime }\left(t\right)}{kT\left(t\right)}=\frac{X^{″}\left(x\right)}{X\left(x\right)}$.

Let $\mu$ be a constant such that $\frac{T^{\prime }\left(t\right)}{kT\left(t\right)}=\mu$ and $\frac{X^{″}\left(x\right)}{X\left(x\right)}=\mu$.

Or, $X^{″}-\mu X=0$ and $T^{\prime }-\mu kT=0$ (3)

Incident-I: If $\mu =0$, the equation becomes $X^{″}=0$. It has a solution of general form, $X\left(x\right)=Ax+B$ through boundary conditions $0=X\left(0\right)=B$ and $0=X\left(L\right)=AL+B$. So $A=B=0$, then the solution becomes a zero solution.

Incident-II: For $\mu ={\lambda }^2$, we get $X^{″}-{\lambda }^{2}X=0$, also $T^{\prime }-{\lambda }^{2}kT=0$.

Then the general solution is formulated by $X\left(x\right)=C_1{\text{e}}^{\lambda x}+C_2{\text{e}}^{-\lambda x}$ and $T=C{\text{e}}^{{\lambda }^{2}kt}$.

And applying the boundary conditions, $0=X\left(0\right)=C_1+C_2\Rightarrow C_1=-C_2$ and $0=X\left(L\right)=C_1{\text{e}}^{\lambda L}+C_2{\text{e}}^{-\lambda L}$. Therefore $0=C_1\left({\text{e}}^{\lambda L}-{\text{e}}^{-\lambda l}\right)$.

This indicates that $C_1=0$, also $C_2=0$ which implies $X\left(x\right)=0$.

Incident-III: If $\mu =-{\lambda }^2$ then (3) becomes $X^{″}+{\lambda }^{2}X=0$ and $T^{\prime }+{\lambda }^{2}kT=0$. Then the general solution of the form $X\left(x\right)=A\cos \lambda x+B\sin \lambda x$, also $T=C{\text{e}}^{-{\lambda }^{2}kt}$ with the boundary conditions $X\left(0\right)=0$ which means $A=0$.

And $X\left(L\right)=0$ becomes $B\sin \lambda L=0\Rightarrow \sin \lambda L=0$.

Which implies $\lambda L=n\pi$.

This leads to ${\lambda }^2=\frac{n^2{\pi }^2}{L^2}$.

Therefore $X\left(x\right)=B\sin \left(\frac{n\pi x}{L}\right)$.

Considering $B=1$, ${\mu }_n=\frac{n^2{\pi }^2}{L^2}$ and $X_n\left(x\right)=\sin \left(\frac{n\pi x}{L}\right)$.

Here ${\mu }_n$ is an Eigen-value of the Sturn-Liouville problem and $X_n\left(x\right)$ is an eigenfunction.

The group of eigenvalues and eigenfunctions for $n\to 1$ to $n\to \infty$ is the entire solution to the Sturn-Liouville problem.

Therefore the solution becomes, $u_n\left(x,t\right)={\text{e}}^{-\frac{n^2{\pi }^{2}kt}{L^2}}\sin \left(\frac{n\pi x}{L}\right)$ that meets the boundary conditions.

Any linear combination of two solutions is, once again, a solution.

So, $u\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{B}_n{\text{e}}^{-\frac{n^2{\pi }^{2}kt}{L^2}}\sin \left(\frac{n\pi x}{L}\right)}$ (4)

that is the solution of the heat equation. Using initial condition the coefficient B_n becomes,

$u\left(x,0\right)=f\left(x\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{B}_n\sin \left(\frac{n\pi x}{L}\right)}$ implies $B_n=\frac{2}{L}{\displaystyle \underset{0}{\overset{L}{\int }}f\left(x\right)}\sin \left(\frac{n\pi x}{L}\right)\text{d}x$.

As a result, Equation (4) is the heat equation’s analytic solution.

1.2. Thermal Diffusivity

The one-dimensional Heat Equation is as follows:

$\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}$

where k is Thermal diffusivity [10]. Thermal diffusivity is different for different materials. Thermal diffusivity is a material attribute that specifies how quickly heat flows through a substance. Thermal diffusivity measures how quickly heat moves through a substance. Materials with high thermal diffusivity quickly regulate their temperature to the temperature around them, as they conduct heat faster. The thermal diffusivity of a substance is calculated by dividing its thermal conductivity by its density and specific heat capacity at constant pressure.

$k=\frac{K}{\rho c}$

where K = Thermal conductivity of the rod’s material (W/m·K).

$\rho$ = Density or the mass of the rod per unit volume (K/m³);

c = the rod’s specific heat capacity (J/kg·K).

The denominator $\rho c$ can be considering the volumetric heat capacity (J/m³·K).

The units of k must be: $\frac{{\left(\text{length}\right)}^{\text{2}}}{\text{time}}$.

Thermal diffusivity is typically measured in units of m²/s, cm²/s, mm²/s.

The thermal diffusivity is calculated using the material’s properties. Table 1 below shows typical values for thermal diffusivity.

Table 1. Value of thermal diffusivity for some materials.

1.3. Numerical Method

The explicit center finite difference technique, implicit finite difference method, and Crank-Nicolson finite difference method are three types of finite difference methods for heat equations. This paper presents the explicit center difference method, a well-known and well-understood numerical method that is required for the heat equation [11].

The finite difference approach for solving the heat equation is described next. Consider a rectangle $\left[0,L\right]\times \left[0,T\right]$ into a finite number of nodes $\left(x_i,t_n\right)$. The notation $u_i^n=u\left(x_i,t_n\right)$ to represent the numerical solution on a finite rod $x\in \left[0,L\right]$ and $t\in \left[0,T\right]$.

1.4. Finite Difference Method for the Heat Equation

The explicit centered difference scheme, which is obtained by a first-order forward difference in time and a second-order centered difference in space, is the simplest numerical discretization of the heat equation [12].

The forward difference in time: $\frac{\partial u}{\partial t}\left(x_i^n\right)\approx \frac{u_i^{n+1}-u_i^n}{\Delta t}$.

The centered difference in space: $\frac{{\partial }^{2}u}{\partial x^2}\left(x_i^n\right)=\frac{1}{{\left(\Delta x\right)}^2}\left[u_{i-1}^n-2u_i^n+u_{i+1}^n\right]$.

Therefore the explicit centered difference scheme of the 1-D heat equation is

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=\frac{k}{{\left(\Delta x\right)}^2}\left[u_{i-1}^n-2u_i^n+u_{i+1}^n\right]$

$\Rightarrow u_i^{n+1}=u_i^n+r\left[u_{i-1}^n-2u_i^n+u_{i+1}^n\right]$ where $r=\frac{k\Delta t}{{\left(\Delta x\right)}^2}$

$\Rightarrow u_i^{n+1}=\left(1-2r\right)u_i^n+r\left[u_{i-1}^n+u_{i+1}^n\right]$.

Stability Condition

If $0<r<\frac{1}{2}$ then the solution at the new time is a weighted average of the solution at the old time. This implies a discrete maximum principle $\max \left|u_i^{n+1}\right|\le \max \left|u_i^n\right|$ and therefore numerical stability. The explicit centered difference scheme is stable if $r<\frac{1}{2}$ that is $\frac{k\Delta t}{{\left(\Delta x\right)}^2}<\frac{1}{2}$.

2. Numerical Computations and Results

Here a computer program (code) in MATLAB Scientific programming language is developed and implement the explicit centered difference method for heat equation for some initial and Dirichlet boundary conditions. In this paper, seven different materials such as Silver, Gold, Copper, Aluminum, Cast Iron, Granite, and Brick are taken and these seven different materials will provide seven different numerical solutions. For the seven materials the values of the thermal diffusivities that is k are different which are given in Table 1. So, these different thermal diffusivities with initial and boundary conditions provide seven different one-dimensional initial-boundary value problems (IBVP). Consider the same initial and boundary condition for the described IBVP. Also, we compare these seven results and analysis the stabilities (Figure 1).

Figure 1. Initial condition for temperature distribution of heat equation for seven materials.

$u\left(x,0\right)=\{\begin{array}{l}0\text{if}0\le x\le 50\\ 10\text{if}50\le x\le 100\end{array}$

and the boundary conditions:

$u\left(0,t\right)=0$ for all t;

$u\left(100,t\right)=10$ for all t.

Case-I

Here we take total time 5 minutes that is 300 seconds and the length of the rod is 100 cm. Total time step 12,000 and time step size, $\Delta t=\frac{300}{12000}=0.025$.

The total length step 40 and space step size, $\Delta x=\frac{100}{40}=2.5$.

Problem 1: we consider Silver material.

(a)(b)

Figure 2. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Silver material. (b) Variation of temperature distribution along rod with time for Silver material.

Problem 2: we consider Gold material.

(a)(b)

Figure 3. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Gold material. (b) Variation of temperature distribution along rod with time for Gold material.

Problem 3: we consider Copper material.

(a)(b)

Figure 4. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Copper material. (b) Variation of temperature distribution along rod with time for Copper material.

Problem 4: we consider Aluminium material.

(a)(b)

Figure 5. (a) Variation of TEmperature distribution along rod according to length at various time steps with initial condition for Aluminium material. (b) Variation of temperature distribution along rod with time for Aluminium material.

Problem 5: we consider Cast Iron material.

(a)(b)

Figure 6. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Cast Iron material. (b) Variation of temperature distribution along rod with time for Cast Iron material.

Problem 6: we consider Granite material.

(a)(b)

Figure 7. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Granite material. (b) Variation of temperature distribution along rod with time for Granite material.

Problem 7: we consider Brick material.

(a)(b)

Figure 8. (a) Variation of temperature distribution along rod according to length at various time steps with initial condition for Brick material. (b) Variation of Temperature distribution along rod with time for Brick material.

Figure 2(a), Figure 3(a), Figure 4(a), Figure 5(a), Figure 6(a), Figure 7(a) and Figure 8(a) show the variation of temperature distribution along rod according to length at various time steps with the above described initial and boundary conditions for the given seven materials such as Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick respectively.

Figure 2(b), Figure 3(b), Figure 4(b), Figure 5(b), Figure 6(b), Figure 7(b) and Figure 8(b) show the variation of temperature distribution along rod with time for Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick respectively.

In these cases time-step $\Delta t=0.025$ and space step $\Delta x=2.5$ then the values of “r” are 0.00684, 0.00508, 0.00456, 0.00344, 0.00048, 0.000044, and 0.0000152 for Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick respectively, which satisfy the stability condition that is, $r=\frac{k\Delta t}{{\left(\Delta x\right)}^2}<0.5$ (Figure 9).

Figure 9. Comparison of Heat distribution between seven different materials w.r.t. initial stage for $\Delta t=0.025$ after 300 seconds.

Case-II

Here we take total time 5 minutes that is 300 seconds and the length of the rod is 100 cm. Total time step 36,000 and time step size, $\Delta t=\frac{300}{36000}=0.0083$.

The total length step 40 and space step size, $\Delta x=\frac{100}{40}=2.5$

In these cases time-step $\Delta t=0.0083$ and space step $\Delta x=2.5$ then the values of “r” are 0.0022708, 0.0016865, 0.0015139, 0.00114208, 0.00015936, 0.000014608, and 0.00000504 for Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick respectively, which satisfy the stability condition that is, $r=\frac{k\Delta t}{{\left(\Delta x\right)}^2}<0.5$ (Figure 10).

Figure 10. Comparison of Heat distribution between seven different materials w.r.t. initial stage for $\Delta t=0.0083$ after 300 seconds.

Case-III

Here we take total time 20 minutes that is 1200 seconds and the length of the rod is 100 cm. Total time step 12,000 and time step size, $\Delta t=\frac{1200}{12000}=0.1$.

The total length step 40 and space step size, $\Delta x=\frac{100}{40}=2.5$

Figure 11. Comparison of Heat distribution between seven different materials w.r.t. initial stage for $\Delta t=0.1$ after 300 seconds.

In these cases time-step $\Delta t=0.1$ and space step $\Delta x=2.5$ then the values of “r” are 0.02736, 0.02032, 0.01824, 0.01376, 0.00192, 0.000176, and 0.0000608 for Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick respectively, which satisfy the stability condition that is, $r=\frac{k\Delta t}{{\left(\Delta x\right)}^2}<0.5$ (Figure 11).

3. Conclusions

In this paper, one-dimensional heat equation as an Initial Boundary Value Problem (IBVP) for the seven different materials such as Silver, Gold, Copper, Aluminium, Cast Iron, Granite, and Brick with seven different thermal diffusivities are considered. The numerical method of the parabolic partial differential equation such as an explicit centered difference scheme with the stability conditions and advantages has been described. By using the MATLAB codes and at the same time step size, the sample tests were performed with Silver, Gold, Copper, Aluminium, Cast Iron, Granite and Brick materials and seven different results were found. The most efficient results were found with Silver, Gold and Copper materials, the third one was Aluminium whereas, efficiency of Cast Iron was found medium. However, the results with Granite and Brick materials were not at all satisfactory. But when we decrease the time step then we were obtained the same result for all materials. Also for better results with Granite and Brick materials, we increase the total time but we were found a small change.

From the definition of the thermal diffusivity and the comparative results of this analysis, we observed that Gold, Silver and Copper materials act as the best heat conductor compared to Aluminium and Cast Iron. Aluminium was found the second good conductor and Cast Iron was found medium whereas, Granite and Brick materials conduct very less amount of heat.

Explicit scheme is conditionally stable. It depends on time step size, $\Delta t$. From stability conditions, we obtained that for Silver, Gold, Copper, Aluminium, Cast Iron, Granite and Brick materials stability conditions are satisfied.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Marwah, J. and Chopra, M.G. (1992) Transient Heat Transfer in a Slab with Heat Generation. Defence Science Journal, 32, 143-149. https://doi.org/10.14429/dsj.32.6274
[2]	Elias, E.A., Cichota, R., Torrioni, H.H. and van Lier, Q.D.J. (2004) Analytical Soil-Temperature Model: Correction for Temporal Variation of Daily Amplitude, Soil Science Society of America Journal, 68, 784-788. https://doi.org/10.2136/sssaj2004.7840
[3]	de Monte, F. (2000) Transient Heat Conduction in One-Dimensional Composite slab: A ‘Natural’ Analytic Approach. International Journal of Heat and Mass Transfer, 43, 3607-3619. https://doi.org/10.1016/S0017-9310(00)00008-9
[4]	Lu, X., Tervola, P. and Viljanen, M. (2005) An Efficient Analytical Solution to Transient Heat Conduction in a One-Dimensional Hollow Composite Cylinder. Journal of Physics A: Mathematical and General, 38, 10145-10155. https://doi.org/10.1088/0305-4470/38/47/007
[5]	Rozin Khatun, M. and Shajib Ali, Md. (2020) Stability Analysis of Finite Difference Schemes for Heat Equation with Various Thermal Conductivity. International Journl of Scientific & Engineering Reserch, 11, 1275-1282.
[6]	Čiegis, R., Jankevičiūtė, G. and Suboč, O. (2010) Numerical Simulation of the Heat Conduction in Composite Materials. Mathematical Modelling and Analysis, 15, 9-22. https://doi.org/10.3846/1392-6292.2010.15.9-22
[7]	Yang, X.-J. and Gao, F. (2017) A New Technology for Solving Diffusion and Het Equations, Thermal Science, 21, 133-140. https://doi.org/10.2298/TSCI160411246Y
[8]	Maturi. D., Alsulami, N. and Alaidarous, E. (2020) Finite Difference Approximation for Solving Transient Heat Conduction Equation of Copper. Advances in Pure Mathematics, 10, 350-358. https://doi.org/10.4236/apm.2020.105021
[9]	Crank, J. and Nicolson, P. (1947) A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat Conduction Type. Mathematical Proceedings of the Cambridge Philosophical Society, 43, 50-67. https://doi.org/10.1017/S0305004100023197
[10]	Gerald, C.F. and Wheatley, P.O. (2002) Applied Numerical Analysis. Pearson Education, Singapore.
[11]	Dhawan, S. and Kumar, S. (2009) A Numerical Solution One Dimensional Heat Equation Using Cubic B-Spline Basis Functions. International Journal of Research and Reviews in Applied Science, 1, 71-77.
[12]	Leveque, R.J. (1992) Numerical Methods for Conservation Law. Birkhauser Verlag, Germany. https://doi.org/10.1007/978-3-0348-8629-1