Finite Difference Method Applied in Two-Dimensional Heat Conduction Problem in the Permanent Regime in Rectangular Coordinates

Vânia Gonçalves de Brito dos Santos¹, Paula Tamires Gomes dos Anjos²

¹Departamento de Ciências Exatas e da Terra, Universidade do Estado da Bahia (UNEB), Salvador, Brazil.
²Centro Educacional Expoente (CEDEX), Brasília, Brazil.

DOI: 10.4236/apm.2022.129038                

Abstract

The solution of many conduction heat transfer problems is found by two-dimensional simplification using the analytical method since different points have different initial temperatures. The temperature at each point of a given element can be analyzed through the Heat Equation that, in some cases, converges to analytical solutions without precision and is far from the real. However, with the application of the Finite Difference Method (FDM), it is possible to solve it numerically in a relatively fast way, providing satisfactory results for the most varied boundary conditions and diverse geometries, characteristics of heat transfer problems by conduction. This study solved two problems inside a plate with and without heat generation involved in temperature distribution. Algorithms were built with the aid of the Matlab programming language, and applied to obtain a numerical solution using the FDM numerical method. The computational and analytical solutions were then compared. Under certain conditions of the parameters involved in the phenomenon of each problem, the numerical method was very efficient for presenting errors less than or equal to 0.003 and 0.03, respectively, for cases without and with heat generation.

Keywords

Heat Transfer, Conduction, Finite Difference Method, Numerical Method, Matlab

1. Introduction

Heat conduction is a physical process in which heat propagates from a temperature source along a section of area or volume, producing thermal gradients. This phenomenon can be classified as a steady state when the heat transmitted in a system is constant and only the temperature varies at each point in the system, or transient, the heat flow in a system is considered transient in which the temperature at various points in the system varies with time [1].

This study aimed to analyze this process of heat conduction or transmission through a problem situation of a flat square plate, with a heating temperature source on one of its edges, determining the temperature field that develops under the conditions imposed on its borders and model such a problem situation.

The study of heat flow using numerical models is a powerful tool because of their advantages over experimental methods and also interesting in view of their low cost and accuracy. There are several numerical methods that can be used to analyze this problem, such as the finite difference method—FDM, which solves partial differential equations by discretizing the continuous physical domain in a finite discrete mesh.

2. Theoretical Background

2.1. Heat Transfer

Heat is a thermal movement that facilitates the exchange of energy from one body to another due to the difference in temperatures in space. Through the energy balance of a material, based on the law of conservation of energy, we can study heat transfer. Depending on the magnitude of the heat transfer in each direction and the determination of the problem solution, heat transmission can be classified into one-direction, two-direction and three-direction.

In the applications presented here, the two-dimensional (2D) mesh conduction mechanisms are addressed, considering an internal region in a steady state in rectangular coordinates using the finite difference method. Through the temperature distribution in a medium, it is possible to determine the energy flow $\stackrel{\dot{}}{Q}$ in two cases, without heat generation resulting in a zero-energy variation as in Equation (1), and with heat generation resulting in the sum of the generated energies in Equation (2):

${\stackrel{\dot{}}{Q}}_{\text{inflow}}-{\stackrel{\dot{}}{Q}}_{\text{outflow}}=0$ (1)

${\stackrel{\dot{}}{Q}}_{\text{inflow}}+{\stackrel{\dot{}}{Q}}_{\text{generated}}-{\stackrel{\dot{}}{Q}}_{\text{outflow}}=0$. (2)

2.2. Heat Conduction in Permanent Regime

Aiming at balancing temperature differences, the thermal conduction mechanism can be understood as a gradual transmission of thermal agitations in different ways [2]. In solids, it occurs through the propagation of vibration of the most energetic molecules.

Considering a rectangular plate subjected to boundary conditions, we can quantify the heat flow from the region of higher to lower temperature, applying the first law of thermodynamics (conservation of energy) together with Fourier’s Law [3]. So, we have the equation:

${\stackrel{\dot{}}{Q}}_{\text{Conduction}}=-kA\frac{\partial T}{\partial x},$ (3)

where ( $\stackrel{\dot{}}{Q}$ ) is given in (W), and A in (m²) is the surface area over which transfer takes place, $\partial T/\partial x$ is the temperature gradient in the x-direction, is given in k/m. The k (W/m∙˚C) material property depends on the material studied, and is known as thermal conductivity.

Considering the amount of heat absorbed in the change of temperature, equivalent to the heat that enters and leaves a material of heat propagation in solid bodies, the heat flow in Equation (3) can be submitted in the expression of Equation (4), this Fourier introduces the ρ density of the substance and the c specific heat. The latter is essential to understand the process of heat propagation as a function of a given time t in three coordinates (x, y, z):

$\frac{\partial T}{\partial t}=\frac{k}{c\rho }\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$. (4)

2.3. Finite Difference Method (FDM)

Analytical solutions of differential equations exist for a limited number of problems with boundary conditions and diverse geometries. The most realistic solutions must be numerically determined at discrete points in the domain, called mesh points. These points are equally spaced; that is, the mesh is regular. The FDM consists of discretizing the internal surface of a material and replacing the derivatives present in the differential equation with approximations using only the numerical values of the function.

In order to obtain the temperature at the mesh point (m, n), Figure 1 illustrates an example of a two-dimensional mesh, however it depends on the temperatures of the mesh points around it. Displacements in the horizontal direction are identified by m, and those in the vertical direction by n. The closer the points

Figure 1. Domain $D\left(x,y\right)\subset ℝ$ —finite difference mesh. (Source: p. 4, [4]).

(Δx and Δy), the greater the accuracy of the method. In this case, the substitutions of the first and second order derivatives, by the finite difference method, adapted from [5], are presented in Equation (11) and Equation (14) corresponding to the value of each point of the temperatures of the material rectangular and an internal environment.

The Finite Difference technique seeks to write the differential operators in their discrete form, that is, as a function of point values of the solution. For explanation, the derivative of a function $u\left(x\right)$ at a point $x_i$ is given by:

${\frac{\partial T}{\partial x}|}_{x=x_i}=\underset{h\to 0}{\lim }\frac{u\left(x_i+h\right)-u\left(x_i\right)}{h},$ (5)

where $h=\Delta x$. From this limit, taking $h\ne 0$ small (not too small to avoid catastrophic cancellation) is expected to obtain a reasonable approximation for the derivative [6]. Approximately, using an increment h, however finite, we can denote the progressive finite difference, being:

${u^{\prime }}_i\approx \frac{u_{i+1}-u_i}{h}.$ (6)

2.4. Variable Separation Method

Fourier’s Law of Equation (3) can also be applied to the energy balance of a two-dimensional material, assuming that there is only heat transfer by conduction under boundary conditions, giving rise to the general heat conduction according to Equation (7). The term $\stackrel{\dot{}}{q}$ represents the rate of energy generated by an infinitesimal volume element, while $\alpha$ is the thermal diffusivity of the material, t is the time of the material diffusivity in the transient term and k is the conductivity of the material studied [4]:

$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{\stackrel{\dot{}}{q}}{k}=\frac{1}{\alpha }\frac{\partial T}{\partial t}.$ (7)

Considering an internal region of the material, with no convection effect, in a steady state and zero heat generation, it is possible to formulate the following heat equation:

$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}=0.$ (8)

With substitutions of the second order derivative by the finite difference method, we obtain:

$\frac{{\partial }^{2}T}{\partial x^2}=\frac{T_{m+1,n}+T_{m-1,n}-2T_{m,n}}{\Delta x^2}$ (9)

$\frac{{\partial }^{2}T}{\partial y^2}=\frac{T_{m,n+1}+T_{m,n-1}-2T_{m,n}}{\Delta y^2}$ (10)

As we determined that the values Δx and Δy are equal, and substituting in Equation (8), we have the calculation of the temperature at a point without heat generation in an internal environment in the steady state:

$T_{m,n}=\frac{T_{m+1,n}+T_{m-1,n}+T_{m,n+1}+T_{m,n-1}}{4}.$ (11)

Analyzing a discretized material with heat generation inside the piece in a steady state, since the real direction of the heat flow (into the mesh) is often unknown; it becomes convenient to formulate the energy balance considering that all the heat flow is inside the point.

Since we do not know the heat transfer that takes place in a control volume around the point (m, n), all heat conductions at the adjacent points are the input energies ( ${\displaystyle \sum {\stackrel{\dot{}}{q}}_{\text{conduction}}}$ ). As there is no convection effect, the generated energies ( ${\stackrel{\dot{}}{E}}_{\text{generated}}$ ) are equal at all points. The proper form for the energy conservation equation is then:

${\stackrel{\dot{}}{E}}_{\text{inflow}}+{\stackrel{\dot{}}{E}}_{\text{generated}}=0.$ (12)

For two-dimensional conditions with the generation, the heat exchange is exerted by conduction between (m, n) and its four adjacent points and considering the unit depth (volume), their points around it are summarized as:

${\displaystyle {\sum }_{i=1}^4\stackrel{\dot{}}{q}cond}+\stackrel{\dot{}}{q}.\left(\Delta x.\Delta y.1\right).$ (13)

Considering the values adopted for Δx and Δy, and substituting the sum of the energies generated at points adjacent to (m, n), we have the calculation of the temperature at a point with heat generation and an internal environment in the steady state.

$T_{m,n}=T_{m+1,n}+T_{m-1,n}+T_{m,n+1}+T_{m,n-1}+\frac{q\Delta x^2}{k}.$ (14)

As stated by [7], the essence of numerical methods is in the discretization of the continuum. It is this discretization that makes the problem finite and thus makes its solution possible through computers.

The representation of the finite difference method in practice is to replace the derivatives by the incremental ratio converging to the value of the derivative when the increment tends to zero and thus, we say that the problem has been discretized. In this study, we deal with the totally discrete domain problem $\left(x,t\right)$, that is, both space (x) and time (t) are discretized. Then, we can define the step mesh $\Delta x,\Delta t$ associated with $\left(x_0,t_0\right)$ as the set of points $\left(x_j,t_n\right)=\left(j\Delta x,n\Delta t\right)$, for all $j=0,1,\cdots ,j+1$ and $n=0,1,\cdots ,n+1$.

The essential mathematical tool in calculating approximations for derivatives is the Taylor Series [8]. The heat conduction energy given by Equation (3) requires approximations for the first derivative in time and second derivative in space. Therefore, assuming that there is the conductivity of material k without total loss and applying the progressive finite difference operators, we have:

$\frac{u_j^{n+1}-u_j^n}{\Delta t}-\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{\Delta x^2}=0.$ (15)

So, from that and taking $\theta =\frac{\Delta t}{\Delta x^2}$, we have:

$u_j^{n+1}=\left(1-2\theta \right)u_j^n+\theta \left(u_{j+1}^n+u_{j-1}^n\right),\forall n\in N.$ (16)

3. Methodology—Problem Analysis

Two heat conduction problems are presented, considering a steady state, which was analyzed by two-dimensional simplifications at the internal point of a mesh: with and without heat generation. They were subjected to a computational procedure, through the Matlab software, version 2017b. Initially, an algorithm was constructed, seeking a solution to the problem of the heat equation only using FDM on a square plate in order to know the mathematical behavior of this heat distribution.

3.1. Numerical Simulation Results

When creating an algorithm that considers only finite difference derivative calculations using heat conduction without energy generation (11), with dimensions $A=20\text{m}\times 20\text{m}$, and considering $k=1\left(\text{W}/\text{m}\cdot ˚\text{C}\right)$, the temperature reached in the middle of the plate was 37.5˚C The heat distribution in Figure 2 shows a quadratic behavior in temperature distribution.

3.2. Computational Modeling

For better stability in the behavior of the program via numerical analysis, it was found necessary to build a quadratic function to model this problem of temperature distribution on the plate. Then, a function was used to represent these data, with which it was possible to determine a temperature matrix at the points of the mesh of this plate, obtaining the analytical model. Then, a second algorithm was built using FDM to evaluate the temperatures of the mesh on the plate, obtaining data for the numerical model. Then, a comparison was made between the two models: analytical and numerical, evaluating the parameters responsible for their approximation and distance.

Figure 2. Temperature distribution on the 2D plate by FDM. (Source: prepared by the authors).

It is worth noting that the continuous heat source used to simulate heat distribution in this problem was an electric heater applied to the upper edge of the plate [L, t], where L is the distance corresponding to the total length of the plate given in meters, x is the measure of the abscissa traveled to the studied point and t is the variable which coincides with the ordinate at the studied point, this resulted in a function that fits the model of temperature variation, in which the highest temperature reached is located in the middle of the plate. In this way, a quadratic function with an exponential coefficient was created:

$\mu \left(x,t\right)=\left(\frac{-T}{L}\right)\ast x\ast \left(x-L\right)\ast \exp \left(-0.545\ast \left(\frac{t}{2}\right)\right)$. (17)

The model is presented on a uniform flat geometry 2D, square with an area of 400 m², exposed to different temperatures in its boundary region, which was divided into 100 parts horizontally, characterizing a discretization in space [0, L] and in 20 parts vertically representing the discretization for the time domain [0, t], as shown in Figure 3.

The Cartesian product of discretization gives the domain [0, L] × [0, t]. To reach the highest temperature in the middle of the plate, we considered the upper edge a temperature ranging from 85˚C to 500˚C, in the lower edge, a temperature of 25˚C, and on the left and right edges, a temperature of 25˚C. With this, we obtained a mesh of coordinate points represented by $\left(x_j,t_n\right)=\left(j\Delta x,n\Delta t\right)$, where:

$\Delta x=\frac{L}{N+1}\text{and}\Delta t=\frac{T}{j+1}.$ (18)

The algorithm script built follows the flowchart in Figure 4, which calculates the temperatures of the mesh nodes referring to the first column of the matrix $u\left(x,t\right)$, using the model equation given by Equation (17) where $t=10$ (ordinate of the mesh midpoint), and from the second column, the FDM is applied, which solves the problem using Equation (16). And the temperature distributions with and without heat generation are obtained, which can be seen in Figure 5 and Figure 6.

Figure 3. Domain [0, L] × [0, t] discretized: [100 × 0.20] × [1 × 20]. (Source: prepared by the authors).

Figure 4. Flowchart of the script used to solve the numerical case without heat generation by the FDM. (Source: prepared by the authors).

Figure 5. Graphic with numerical solution by the FDM with temperature source without heat generation: initial condition T= 100˚C. (Source: prepared by the authors).

Figure 6. Graphic with numerical solution by the FDM with temperature source without heat generation: initial condition T = 100˚C. (Source: prepared by the authors).

4. Results and Discussion

4.1. Without Heat Generation

For comparison between the analytical model given by Equation (17) and the computational model given by Equation (16), we considered $k=1\left(\text{W}/\text{m}\cdot ˚\text{C}\right)$, the values of temperatures reached at the midpoint of the plate are listed in Table 1.

The results obtained by the numerical method were very close to the calculated analytical values. Thus, expressions obtained by the FDM through the energy balance can be considered valid for numerical analysis of temperatures. Figure 7 and Figure 8 show that the numerical model used is validated by the good results obtained.

4.2. With Heat Generation

A simulation of temperature distribution on the plate was also performed, considering heat generation, as well as in the situation without heat generation, we also used the hypothetical situation of a plate whose material is a type of concrete, with $k=1.8\left(\text{W}/\text{m}\cdot ˚\text{C}\right)$. In order to adapt to the analytical model based on Equation (17), the addition of the term of energy generated as the heat was introduced, following the same idea presented in heat conduction with energy generation as in Equation (14). The term is given by $\left(q\cdot \Delta x\right)/k$, where $q=\Delta t=T_{\text{reachedatthemidpoint}}-T_{\text{initial}}$ and $\Delta x=0.2$.

The results are presented in Table 2 with the temperature values obtained in

Figure 7. Comparison between analytical and numerical data without heat generation. (Source: prepared by the authors.

Figure 8. Relative error of numerical data without heat generation. (Source: prepared by the authors).

Table 1. Comparasion between FDM temperature without heat generation.

the problem simulation. There was a greater increase in values of temperatures reached analytically and there was a gap between the results obtained by the two methods, as the source temperature increased (Figure 9).

This caused a larger error (Figure 10) between the values obtained by the analytical and numerical methods. So, the numerical method did not achieve high accuracy in the heat distribution process in this case, but this can be justified by the insertion of the heat generation term, since this increase was not linear, it depended on the variation between the initial temperatures and the one reached at the midpoint. There was a resistance of the algorithm to closely follow high values of temperature, but even so, results with an acceptable error were obtained.

Figure 9. Comparison between analytical and numerical results with heat generation. (Source: prepared by the authors).

Table 2. Comparison between FDM temperature with heat Generation.

4.3. Relative Error Evaluation

According to Figure 8 and Figure 10, using the Euclidean norm, the relative errors were estimated, for the cases with and without heat generation, indicating values smaller than or equal to 0.003, and 0.03, respectively, evidencing a better accuracy for the case without heat generation (Figure 11). Probably, in the simulation with heat generation, besides the non-linearity of the term added in the equation of the analytical model, the lower accuracy is justified by the resistance of the algorithm to follow higher temperatures.

Figure 10. Relative error of numerical data with heat generation. (Source: prepared by the authors).

Figure 11. Relative error between analytical and numerical data without and with heat generation. (Source: prepared by the authors).

5. Final Considerations

The initial objective of this study was to use the finite difference method (FDM) to solve the problem of heat conduction in a plate under steady state conditions. It was realized that it was possible to develop a model. Numerical tests were then performed, by two-dimensional simplifications, at the internal point of a mesh with and without heat generation. The results show proximity between the results in the analytical and numerical calculations, determining the temperature value reached the midpoint of the plate. Even for the case in which heat generation in the temperature distribution process was considered, the result was still satisfactory. The errors found using the Euclidean norm indicate values smaller than or equal to 0.003, in the simulations without heat generation, and errors smaller than or equal to 0.03, in simulations with heat generation.

The first study done by an algorithm that only evaluated the temperatures by the FDM on the plate, indicated that the temperature distribution behavior was similar to that of a quadratic function, this information was fundamental for the construction of computational modeling. The function of Equation (17) used to calculate the first column of the discretized numerical data of plate temperatures, in the models, is quadratic following what is already known in the literature, but something novel was used, an exponential coefficient that depends on the variable t, which was taken as the midpoint of ordinates in the mesh, it served to dampen the sudden rise of temperatures in the algorithm process, avoiding divergence.

With the results obtained, it can be stated that the FDM is a practical and efficient choice for thermal analysis since it has good precision and accuracy, especially in the case where there is no heat generation. In addition, its background and formulation are understandable, which facilitates the use of the computational resource, as the distributions of temperatures and thermal gradients in graphic form were evident when generated by the developed algorithm, facilitating the problem analysis and understanding.

Acknowledgements

The authors acknowledge the contributions of the Departamento de Ciências Exatas e da Terra (DCET-I), of the Universidade do Estado da Bahia (UNEB), for encouraging the research work developed in this manuscript.

Conflicts of Interest

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

References

[1]	Incropera, F.P., Dewitt, D.P., Bergman, T.L. and Lavine, A.S. (2008) Fundamentos da Transferência do Calor e de Massa. 6th Edition, Livros Técnicos e Científicos Editora S.A., Rio de Janeiro.
[2]	Almeida, G.V., Coelho, N.A. and Pedroso, L.J. (2016) Distribuição de temperatura em placas em regime transiente: comparação entre solução analítica e numérica. XXXVII Iberian-Latin American Congress on Computational Methods in Engineering, Brasília, 6-9 November 2016, 2-4.
[3]	Butkov, E. (1988) Física Matemática. Livros Técnicos e Científicos Editora Ltda, Rio de Janeiro.
[4]	Neto, G.C.S. and Valda, L.H.C. (2009) O método das diferenças finitas aplicado a problemas de transmissão de calor em regime transiente. CILAMCE 2009 Iberian-Latin-American Congress on Computational Methods in Engineering, Búzios, 8-11 November 2009, 4-7.
[5]	Almeida, G.V., Coelho, N.A. and Alkmim, N. (2017) Comparative Analysis of a Transient Heat Flow and Thermal Stresses by Analytical and Numerical Methods. VII International Conference on Computational Methods for Coupled Problems in Science and Engineering, Rhodes Island, 12-14 June 2017, 1204-1206.
[6]	Sauter, E., Justo, D.A.R., de Azevedo, F.S., Guidi, L.F. and Konzen, P.H.A. (2020) Cálculo Numérico, Um Livro Colaborativo Versão Scilab. https://www.ufrgs.br/reamat/CalculoNumerico/livro-sci/livro-sci.pdf
[7]	Cunha, C. (1993) Métodos Numéricos para as Engenharias e Ciências Aplicadas. UNICAMP, São Paulo.
[8]	Franco, N.B. (2006) Cálculo Numérico. Person Prentice Hall, São Paulo.