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.
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 where is the temperature at the point and time . The mesh points in a plane parallel to the 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
where is the current time index, is the step size for time discretization, and , , and are mesh spacing in , , and directions. The total truncation error for the numerical method is
After grouping like terms from Equation (2), we obtain
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 at as [14] where is an amplification factor, , and , , , , and are constants. The iteration scheme in Equation (4) is stable if [15]. Using Equation (5) in Equation (4) and simplifying, we get
where , , , , , , , and .
Using the relations and , Equation (6) becomes
For the worst case, we must have [16]. This gives
The stability condition yields
Hence, the time step size should satisfy Equation (9) to ensure numerical stability. To determine a suitable value for that gives a stable iteration, one can take a minimum of , , and 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 , , and [17]. It is assumed that the stainless steel has density , specific heat capacity , and thermal conductivity [18]. Thus, the thermal conductivity . The following initial and boundary conditions were considered to illustrate the numerical solution method [14, 19, 20].
Initial condition: where is the Bessel function of order one given by [20] and is the first positive zero of .
Boundary condition:
The exact solution is chosen as
Mesh spacing , , and 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 . The time step size of was used in the computation. The approximate value was taken to determine the initial condition [21]. A MATLAB code was developed to perform the numerical computations and to graph the solutions.
Figure 2 is a 3D plot of temperature distribution in the solid cylinder for the time value at obtained by the numerical method and exact solutions. Figure 3 shows temperature contours at different cross-section of the cylinder from the two methods.
(a)
(b)
Using the numerical method, the maximum temperature value on the cylinder obtained at 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 to any time . 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 of and in cylindrical coordinate expression which affects the accuracy of the finite difference scheme near the origin (). It also takes more computational time when we use small mesh spacings in -direction since small 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 at with the two numerical methods and the exact solution are presented in Table 1 for , , and . The absolute errors of the two numerical methods at and 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.
(a)
(b)
If there are , , and nodes in , , and 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 , , and . 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
where 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.
From Table 2, 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The author declares that there is no conflict of interest regarding the publication of this article.