# Numerical Solution of Three-Dimensional Transient Heat Conduction Equation in Cylindrical Coordinates

Authors: Endalew Getnet Tsega
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Heat equation is a partial differential equation used to describe the temperature distribution in a heat-conducting body. The implementation of a numerical solution method for heat equation can vary with the geometry of the body. In this study, a three-dimensional transient heat conduction equation was solved by approximating second-order spatial derivatives by five-point central differences in cylindrical coordinates. The stability condition of the numerical method was discussed. A MATLAB code was developed to implement the numerical method. An example was provided in order to demonstrate the method. The numerical solution by the method was in a good agreement with the exact solution for the example considered. The accuracy of the five-point central difference method was compared with that of the three-point central difference method in solving the heat equation in cylindrical coordinates. The solutions obtained by the numerical method in cylindrical coordinates were displayed in the Cartesian coordinate system graphically. The method requires relatively very small time steps for a given mesh spacing to avoid computational instability. The result of this study can provide insights to use appropriate coordinates and more accurate computational methods in solving physical problems described by partial differential equations.


---

## Body

## 1. Introduction

In science and engineering, partial differential equations are used to express how some quantity varies with position and time [1]. Heat equation is one of such equations used to describe the variation of temperature in a body. Many diffusion processes can also be described with this equation, for instance, fluid flow in porous media, molecular diffusion in biological tissues, and water uptake by plant roots [2]. Solution of the heat equation is required to understand and analyse the heat transfer problems [3]. In many real situations, the complexity of the considered problem makes it difficult to obtain the analytical solutions. In this case, numerical methods such as finite difference, finite element, and finite volume methods are used to obtain approximate solutions [4–7].Finite difference discretizations in polar or cylindrical coordinates are more convenient than in Cartesian coordinates to solve boundary value problems involving circular shapes because they avoid the use of complicated differentiation formulae near the curved boundaries [8]. Finite difference discretizations in polar or cylindrical coordinates have been used by authors to solve partial differential equations. Mori and Romão [9] used the finite difference method to perform numerical simulation of 2D convection-diffusion in cylindrical coordinates. Iyengar and Manohar [10] used the fourth-order difference method for the solution of Poisson’s equation in cylindrical coordinates. They extended the method to solve heat equation in two-dimensional with polar coordinates and three-dimensional with cylindrical coordinates. Shiferaw and Mittal [11] solved three dimensional Poisson’s equation with the finite difference method in cylindrical coordinates. Salehi and Granpayeh [12] presented a finite difference method for solving the two-dimensional Schrödinger equation in polar coordinates.One purpose of this paper is to present finite difference discretization of transient three-dimensional heat equation in cylindrical coordinates and to obtain more accurate solution by a higher-order finite difference method with a computer program. The other purpose is to display the solution results graphically in the Cartesian coordinate system for the better visualization.


## 2. Finite Difference Discretization of Heat Equation

The transient three-dimensional heat equation in cylindrical coordinates is(1)∂T∂t=α∂2T∂r2+1r∂T∂r+1r2∂2T∂2θ+∂2T∂z2,where Tr,θ,z,t is the temperature at the point r,θ,z and time t. The mesh points in a plane parallel to the r−θ plane are defined by the intersection points of the circles and the straight lines as shown in Figure 1. Approximating the time derivative by two-point forward difference (first-order accurate), second-order spatial derivatives by five-point central differences (fourth-order accurate), and first-order spatial derivative by two-point central difference (second-order accurate) [8, 13], we get
(2)Ti,j,kn+1−Ti,j,knΔt=α−Ti−2,j,kn+16Ti−1,j,kn−30Ti,j,kn+16Ti+1,j,kn−Ti+2,j,kn12Δr2+1riTi+1,j,kn−Ti−1,j,kn2Δr+1ri2−Ti,j−2,kn+16Ti,j−1,kn−30Ti,j,kn+16Ti,j+1,kn−Ti,j+2,kn12θ2+−Ti,j,k−2n+16Ti,j,k−1n−30Ti,j,kn+16Ti,j,k+1n−Ti,j,k+2n12Δz2,Figure 1
Sample cylindrical mesh for a solid circular cylinder.wheren is the current time index, ∆t is the step size for time discretization, and ∆r, ∆θ, and ∆z are mesh spacing in r, θ, and z directions. The total truncation error for the numerical method is
(3)Rri,θj,zk,tn=∂Tri,θj,zk,tn∂tΔt2−α∂3Tri,θj,zk,tn∂r3Δr26+α∂6Tri,θj,zk,tn∂r6Δr490+αri2∂6Tri,θj,zk,tn∂θ6Δθ490+α∂6Tri,θj,zk,tn∂z6Δz490+⋯=OΔt+Δr2+Δr4+Δθ4+Δz4.After grouping like terms from Equation (2), we obtain
(4)Ti,j,kn+1=−αΔt12Δr2Ti−2,j,kn+4αΔt3Δr2−αΔt2riΔrTi−1,j,kn+1−2.5αΔtΔr2−2.5αΔtriΔθ2−2.5αΔtΔz2Ti,j,kn+4αΔt3Δr2+αΔt2riΔrTi+1,j,kn−αΔt12Δr2Ti+2,j,kn−αΔt12riΔθ2Ti,j−2,kn+4αΔt3riΔθ2Ti,j−1,kn+4αΔt3riΔθ2Ti,j+1,kn−αΔt12riΔθ2Ti,j+2,kn−αΔt12Δz2Ti,j,k−2n+4αΔt3Δz2Ti,j,k−1n+4αΔt3Δz2Ti,j,k+1n−αΔt12Δz2Ti,j,k+2n.The iterative scheme, Equation (4), is used to compute the temperature at the grids by applying initial and boundary conditions. In this paper, this scheme is named as the “five-point central difference method.” Similarly, if the second-order derivatives in Equation (1) are discretized by three-point central differences, the corresponding scheme is named as the “three-point central difference method.”


## 3. Stability Analysis

The Von Neumann stability analysis is used to check the stability condition of the numerical method. According to this stability analysis method, we replace the temperature atpΔr,qΔθ,kΔz,nΔt as [14]
(5)Tnp,q,k=Aeiα1pΔreiα2qΔθeiα3kΔzeβnΔt=Aeiα1pΔreiα2qΔθeiα3kΔzξβ,where ξ=enΔt is an amplification factor, i=−1, and A, α1, α2, α3, and β are constants. The iteration scheme in Equation (4) is stable if ξ≤1 [15]. Using Equation (5) in Equation (4) and simplifying, we get
(6)ξ=1−ΔtC1e−2iα1Δr−C2e−iα1Δr+C3−C4eiα1Δr+C1e2iα1Δr+C5e−2iα2Δθ−C6e−iα2Δθ−C6eiα2Δθ+C5e2iα2Δθ+C7e−2iα3Δz−C8e−iα3Δz−C8eiα3Δz+C7e2iα3Δz,whereC1=α/12Δr2, C2=4α/3Δr2−α/2riΔr, C3=2.5α/Δr2+2.5α/riΔθ2+2.5α/Δz2, C4=4α/3Δr2+α/2riΔr, C5=α/12riΔθ2, C6=4α/3riΔθ2, C7=α/12Δz2, and C8=4α/3Δz2.Using the relationseiθ=cosθ+isinθ and e−iθ=cosθ−isinθ, Equation (6) becomes
(7)ξ=1−ΔtC12cos2α1Δr−C2cosα1Δr−sinα1Δr+C3−C4cosα1Δr+sinα1Δr+2C5cos2α2Δθ−2C6cosα2Δθ+2C7cos2α3Δz−2C8cosα3Δz.For the worst case, we must haveα1Δr=π=α2Δθ=α3Δz [16]. This gives
(8)ξ=1−Δt2C1+C2+C3+C4+2C5+2C6+2C7+2C8.The stability conditionξ≤1 yields
(9)Δt≤22C1+C2+C3+C4+2C5+2C6+2C7+2C8.Hence, the time step sizeΔt should satisfy Equation (9) to ensure numerical stability. To determine a suitable value for Δt that gives a stable iteration, one can take a minimum of ri, i=1,2,3,⋯,n, and ri≠0 in Equation (9).


## 4. A Test Problem and Solutions

Consider three-dimensional transient heat conduction which satisfies Equation (1) in a solid cylinder made up of stainless steel having dimensions 0≤r≤1, 0≤z≤2, and 0≤θ≤2π [17]. It is assumed that the stainless steel has density ρ=8000kg/m3, specific heat capacity c=475J/kg·K, and thermal conductivity c=15W/m·K [18]. Thus, the thermal conductivity α=k/ρc=3.9474×10−6m/s. The following initial and boundary conditions were considered to illustrate the numerical solution method [14, 19, 20].Initial condition:
(10)Tr,θ,z,0=Fr,θ,z=J1s11Rr2cosθ+2sinθcosπLze−αs11/R2+π/L2,0≤r≤R=1,0≤z≤L=2,0≤θ≤2π,where J1x is the Bessel function of order one given by [20]
(11)J1x=∑k=0∞−1kx/22k+1k!k+1!,and s11 is the first positive zero of J1.Boundary condition:
(12)T0,θ,z,t<∞,Tr,θ,z,t=Tr,θ+2π,z,t,TR,θ,z,t=0,∂T∂zz=0=0,∂T∂zz=L=0.The exact solution is chosen as(13)Tr,θ,z,t=J1s11Rr2cosθ+2sinθcosπLze−αs11/R2+π/L2t,0≤r≤R=1,0≤z≤L=2,0≤θ≤2π,t>0.Mesh spacingΔr=0.01, Δθ=π/32, and Δz=0.02 were taken in the computation with the numerical method. This results in 656,601 grid points in the solid cylinder. The iterative scheme in Equation (4) cannot be used to approximate the temperature at the origin since it involves division by zero at r=0. The time step size of Δt=0.00005 was used in the computation. The approximate value s11=3.83170597020751 was taken to determine the initial condition [21]. A MATLAB code was developed to perform the numerical computations and to graph the solutions.Figure2 is a 3D plot of temperature distribution in the solid cylinder for the time value t=0.5 at z=0,0.5,1.5,2 obtained by the numerical method and exact solutions. Figure 3 shows temperature contours at different cross-section of the cylinder from the two methods.Figure 2
3D plot of temperature distributions fort=0.5.Figure 3
Temperature contours att=0.5 for z between 0 and 2 with a step of (a) 0.5 and (b) 0.02.
(a)(b)Using the numerical method, the maximum temperature value on the cylinder obtained att=1 was 1.644829743701951. The corresponding value in the exact solution was 1.645649737422976. The maximum absolute error of the numerical method for this time was 0.0008199937210258135. This indicates that the solution by the numerical method was in a good agreement with the exact solution. The MATLAB code enables us to simulate heat flow in the cylinder from t=0 to any time T>0. The numerical method requires relatively very small time steps for a given mesh spacing to avoid computational instability [22].For heat equation, one difference between its representation in Cartesian and cylindrical coordinates is the presence of1/r and 1/r2 in cylindrical coordinate expression which affects the accuracy of the finite difference scheme near the origin (r=0). It also takes more computational time when we use small mesh spacings in r-direction since small ∆t is required in order to guarantee stability. In this paper, the center of grid was avoided in the finite-difference computation and the average of the values on the neighbouring nodes was taken to determine the value at the origin [23, 24]. The five-point central difference method in cylindrical coordinates is very sensitive to choose the time step for stability.The accuracy of the five-point central difference method was compared with that of the three-point central difference method in solving heat equation. The solutions of the heat equation att=2 with the two numerical methods and the exact solution are presented in Table 1 for 0≤r≤1, z=0.5, and θ=π/8. The absolute errors of the two numerical methods at t=0.5 and t=2 in the same location are shown in Figure 4. It can be observed that the five-point central difference method is more accurate than the three-point central difference method in solving heat conduction equation in cylindrical coordinates. But the five-point central difference method took more computational time than the three-point central difference method. The errors in the two methods increase as time increase for the selected problem. The numerical method was also applied to solve Equation (1) and gave results which suggested convergence to the exact solution with small time steps.Table 1
A comparison between three-point central difference, five-point central difference, and exact solutions for0≤r≤1, z=0.5, and θ=π/8 at t=2.
rThree-pointFive-pointExact0.000.0000000000000000.0000000000000000.0000000000000000.040.1411673181693720.1411667580645190.1411668842604900.080.2798517225945590.2798514428279400.2798515040834180.120.4136202400600110.4136200538941260.4136200926899100.160.5401349959082540.5401348567133420.5401348837107310.200.6571997364407460.6571996256168690.6571996451150010.240.7628023894535870.7628022977531320.7628023119374490.280.8551539142646790.8551538364542830.8551538466098280.320.9327225953992810.9327225282562340.9327225352219770.360.9942630783952290.9942630198151590.9942630241852580.401.0388395617089141.0388395102579161.0388395124816720.441.0658426824763371.0658426371450321.0658426375808000.481.0749997665542581.0749997266128901.0749997255588720.521.0663782525021741.0663782174100491.0663782151227430.561.0403822420965271.0403822114427701.0403822081482420.600.9977422735694010.9977422470323890.9977422429333490.640.9394985549788060.9394985322976320.9394985275778560.680.8669780309442700.8669780118980220.8669780067251100.720.7817657835324200.7817657679248710.7817657624519230.760.6856713846527560.6856713723009450.6856713666675320.800.5806909204627600.5806909111884470.5806909055211460.840.4689654958184680.4689654894418410.4689654838544840.880.3527370969176300.3527370932527950.3527370878465690.920.2343027415226060.2343027403745490.2343027352380260.960.1159678774854320.1159678786487050.1159678738573901.000.0000000000000000.0000000000000000.000000000000002Figure 4
Errors in finite difference approximations.
(a)(b)If there areMr, Nθ, and Lz nodes in r, θ, and z directions in a cylindrical mesh, the five-point central difference can be used to approximate second-order spatial derivatives in Equation (1), at the nodes for 1≤i≤Mr−2, 1≤j≤Nθ, and 3≤i≤Lz−2. We can use the three-point central difference approximation at all interior nodes of the mesh. In Cartesian grids, more effort is required to apply the five-point central difference method near the boundaries. But the circular nature of cylindrical grids enables us to apply this method easily and implement on a digital computer to obtain solutions.Let(14)EΔr,Δθ,Δz,Δt=max1≤i≤Nr1≤j≤Nθ1≤k≤Nzn=NtTri,θj,zk,tn−Ti,j,kn,whereNt is the number of time steps. The maximum errors of the numerical solution for different mesh spacing and time sizes are shown in Table 2.Table 2
Maximum errors of the numerical solution.
∆r∆θ∆z∆tNtE∆r,∆θ,∆z,∆t0.02π/160.050.00210004.663020E−030.01π/320.0250.00110004.170669E−040.005π/640.010.000510001.565003E−04From Table2, we observe that error decreases as grid spacing is reduced. This shows the convergence of the numerical method for the selected test problem.


## 5. Conclusion

This paper presented the five-point central difference method to solve the three-dimensional transient heat conduction equation in cylindrical coordinates. The numerical method is capable of computing more accurate solution than the three-point central difference method to the heat conduction equation. The five-point central difference method can be applied in cylindrical grids easily because of the circular nature of the grid locations. In the computation, the singularities of heat conduction equation in cylindrical coordinates were avoided by taking average of values on surrounding grids. To acquire computational stability, a relatively very small time step is required for a given mesh spacing. A similar procedure can be applied to derive other higher-order or implicit finite difference schemes in cylindrical coordinates to increase accuracy and stability.



---

*Source: 1993151-2022-01-04.xml*