Non-Fourier Heat Conduction Analysis with Temperature-Dependent Thermal Conductivity

Abstract

As a first endeavor, the one- and two-dimensional heat wave propagation in a medium subjected to different boundary conditions and with temperature-dependent thermal conductivity is studied. Both the spatial as well as the temporal domain is discretized using the differential quadrature method (DQM). This results in superior accuracy with fewer degrees of freedom than conventional finite element method (FEM). To verify this advantage through some comparison studies, a finite element solution ise also obtained. After demonstrating the convergence and accuracy of the method, the effects of different parameters on the temperature distribution of the medium are studied.

1. Introduction

Most of heat conduction or mass diffusion problems are described and analyzed using Fourier's law and Fick's law. It is well known that these classical diffusion theories may break down when one is interested in the transient problems in an extremely short period of time, very high flux, or for very low temperatures. Then, a modified flux model for the transfer processes with a finite speed wave is suggested. The hyperbolic heat conduction equation based on the Cattaneo model for the heat flux incorporates a relaxation mechanism in order to gradually adjust to a change in the temperature gradient. This model has been a satisfactory extension of classical diffusion theory and can yield the hyperbolic diffusion equation within the continuum assumption.

During the past few years, considerable attention has been concerned with using the non-Fourier heat conduction in problems with very low temperatures, extremely short period of time, or very high heat flux; see, for example, [1–5].

There are some interesting analytical approaches for different non-Fourier heat conduction problems [6–11]. However, due to difficulty in applying analytical schemes to investigate problems with complicated geometries, boundary conditions, or variable thermal properties, attention to numerical schemes has been increased.

Carey and Tsai [12] applied central and backward difference schemes to a one-dimensional problem to examine oscillations for the numerical solution of propagating heat wave reflected at a boundary. Glass et al. [13] used McCormack^,s predictor-corrector method to solve the same problem. Tamma and Railkar [14] used finite element method with special shape functions to remove the numerical oscillations. Han-Taw and Jae-Yuh [15] solved one-dimensional hyperbolic heat equation using Laplace-transform and finite volume technique. They extended the presented approach for two-dimensional problems [16]. Yeung and Lam [17] developed a flux-splitting algorithm based on the Godunov numerical scheme for the solution of the one- and two-dimensional non-Fourier heat conduction equation. Hsu and Chu [18] used the central finite difference method to investigate the two- and three-dimensional inverse non-Fourier heat conduction problem in electronic device. Liu and Chen [19] applied a hybrid application of the Laplace transform and control-volume formulation in-conjunction with the hyperbolic shape functions to investigate the hyperbolic diffusion problems. Wu and Li [20] extended the time discontinuous Galerkin finite element method (DGFEM) developed by Li et al. [21] to the heat wave propagation problem. Wu and Li [22] employed a Galerkin finite element method to simulate the heat wave propagation in the medium subjected to different kinds of heat source. Chen [23] combined the Laplace transform, the weighting function scheme, and the hyperbolic shape functions to investigate the effect of the surface curvature of a solid body on hyperbolic heat conduction. Reverberi et al. [24] proposed a numerical study on a nonlinear hyperbolic diffusion equation. The Hartree hybrid method, combining the finite difference techniques with the method of characteristics, was used in the presence of discontinuities between initial and boundary conditions. A sequential method was proposed to estimate boundary condition of the two-dimensional hyperbolic heat conduction problems by Yang [25]. An inverse solution was deduced from a finite difference method, the concept of the future time and a modified Newton-Raphson method. The numerical solution of Dorao [26] was based on the time-space least-squares spectral method. A spectral element method was applied for solving the hyperbolic system treating the heat flux as an independent variable in addition to temperature. T. M. Chen and C. C. Chen [27, 28] investigated the hyperbolic heat conduction problems in the radial-spherical coordinate system by the hybrid Green's function method. Wang [29] used the finite element method with linear least-squares error method to construct a simple noniterative inverse operation procedure for solving non-Fourier heat transfer effect in irregular geometry. The method combines the Laplace transform for the time domain, Green's function for the space domain, and ε-algorithm acceleration method for fast convergence of the series solution. Andarwa and Basirat Tabrizi [30] studied heat and mass transfer coupling under non-Fourier effect. Hsu [31] employed the conventional differential quadrature method to solve the hyperbolic heat conduction problem with constant material properties.

In this work, as a first endeavor, the incremental differential quadrature method (IDQM) is employed for the nonlinear non-Fourier one- and two-dimensional heat-transfer analysis. The initial conditions are directly implemented into the governing differential equations through the IDQM discretization rules without reformulating the IDQM weighting coefficients. The Newton-Raphson algorithm is employed to solve the nonlinear system of algebraic equations. After demonstrating the convergence and accuracy of the method, some parameter studies are performed.

2. Basic Formulations and Solution Procedure

Consider a rectangular domain with arbitrary thermal conditions along its boundaries. Its thermal conductivity (𝑘) is assumed to be a linear function of temperature [32]𝑘(𝑇)=𝑘₀1+𝜆𝑇−𝑇₀.(1) The other material properties of the domain such as density (𝜌) and specific heat capacity (𝐶_𝑝) are assumed to be constant.

The hyperbolic constitutive equation governing the transient heat transfer may be written as 𝜏𝜕𝑞𝜕𝑡+𝑞=−𝑘∇𝑇,(2) where 𝜏 is the relaxation time, which is related to the thermal wave speed and thermal diffusivity as 𝜏=𝛼/𝐶². Moreover, the energy equation can be written as𝜌𝐶_𝑝𝜕𝑇𝜕𝑡+∇⋅𝑞=𝑔(𝑥,𝑦,𝑡).(3) Using (1) and (2), (3) takes the following form: 𝜕𝜕𝑥𝑘𝜕𝜕𝑇𝑥+𝜕𝜕𝑘𝑦𝜕𝜕𝑇𝑦+𝑘𝜕²𝑇𝜕𝑥²+𝜕²𝑇=𝜏𝜌𝐶_𝑝𝜕𝑦²+𝑔(𝑥,𝑦,𝑡)+𝜏𝜕𝜕𝑔𝑡𝜕²𝑇𝜕𝑡²+𝜌𝐶_𝑝𝜕𝑇𝜕𝑡.(4)

Due to nonlinear terms in (4), if it is not impossible to solve it analytically, it is very difficult to obtain such a solution. Hence, the approximate method should be used to solve it. As a first attempt, the differential quadrature method as an efficient and accurate numerical tool [31–36] is employed to solve (4) under arbitrary boundary and initial conditions.

The differential quadrature method (DQM) is an alternative discretization approach for directly solving the governing equations in mathematics and engineering applications. Here, the DQM is used to discretize both the spatial and temporal domains. According to DQ method, the domain is discretized into a set of 𝑁_𝑥 and 𝑁_𝑦 discrete grid points in the 𝑥- and 𝑦-direction, respectively. Then, at a given grid point (𝑥_𝑖,𝑦_𝑗), the first- and second-order derivatives of an arbitrary function 𝑓(𝑥,𝑦,𝑡) can be approximated as 𝜕𝑓(𝑥,𝑦,𝑡)𝜕𝑥|||_{(𝑥𝑖,𝑦𝑗)}=^𝑁𝑥_𝑚₌₁𝐴^𝑥𝑖𝑚𝑓𝑥_𝑚,𝑦_𝑗,𝑡=^𝑁𝑥_𝑚₌₁𝜕𝐴^𝑥𝑖𝑚𝑓_𝑚𝑗(𝑡),₂𝑓(𝑥,𝑦,𝑡)𝜕𝑥²||||_{(𝑥𝑖,𝑦𝑗)}=^𝑁𝑥_𝑚₌₁𝜕𝑓(𝑥,𝑦,𝑡)𝐵^𝑥𝑖𝑚𝑓_𝑚𝑗(𝑡),𝜕𝑦||||_{(𝑥𝑖,𝑦𝑗)}=^𝑁𝑦_𝑛=1𝜕𝐴^𝑦𝑗𝑛𝑓_𝑖𝑛(𝑡),₂𝑓(𝑥,𝑦,𝑡)𝜕𝑦²||||_{(𝑥𝑖,𝑦𝑗)}=^𝑁𝑦_𝑛=1𝐵^𝑦𝑗𝑛𝑓_𝑖𝑛(𝑡),(5) where 𝑖=1,2,…,𝑁_𝑥 and 𝑗=1,2,…,𝑁_𝑦.

In order to determine the weighting coefficients, a set of test functions should be used in (5). For polynomial basis functions DQ, a set of Lagrange polynomials is employed as the test functions. The normalized weighting coefficients for the first- and second-order derivatives in the 𝜁-direction (𝜁=𝑥,𝑦, and 𝑡) are thus determined, respectively, as 𝐴^𝜁𝑖𝑗=⎧⎪⎪⎨⎪⎪⎩1𝐿_𝜁𝑀𝜁_𝑖−𝜁_𝑖−𝜁_𝑗𝑀𝜁_𝑗𝑖≠𝑗,_{𝑁𝜁𝑘=}_{1,𝑘≠𝑖}𝐴^𝜁𝑖𝑘𝑖=𝑗,𝑖,𝑗=1,2,…,𝑁^𝜁𝐵_𝜁𝑖𝑗,=⎧⎪⎪⎨⎪⎪⎩2𝐴^𝜁𝑖𝑖𝐴^𝜁𝑖𝑗−𝐴^𝜁𝑖𝑗−𝜁_𝑖−𝜁_𝑗𝑖≠𝑗,_{𝑁𝜁𝑘=}_{1,𝑘≠𝑖}𝐵^𝜁𝑖𝑘𝑖=𝑗,𝑖,𝑗=1,2,…,𝑁^𝜁,(6) where, 𝑀𝜁_𝑖=^𝑁𝜁_𝑗=_{1,𝑖≠𝑗}𝜁_𝑖−𝜁_𝑗.(7)

In numerical calculation, Chebyshev-Gauss-Lobatto quadrature points are used to generate the grid points, that is, 𝜁_𝑖=𝐿^𝜁21−cos(𝑁𝑖−_𝜁−1)1𝜋𝑖=1,2,…,𝑁^𝜁,𝜁=𝑥,𝑦,𝑡.(8)

Also, to accurately and computationally efficient discretization of the temporal domain, the incremental DQ method (IDQM) is employed [32, 33, 36]. Based on this approach, the total temporal domain is divided into a set of subdomain and in each subdomain, the DQ rule is employed to discretize the temporal derivatives. At the end of each subdomain, the temperature and its first-time derivative (𝜕𝑇/𝜕𝑡) are used as the initial conditions for the next subdomain.

In order to overcome the drawback of the conventional DQM for initial value problem with derivative of order greater or equal to two, which results in multiple initial conditions at the beginning of each temporal subdomain, the proposed method by Malekzadeh et al. [36] is employed. According to this methodology, the initial conditions are directly implemented through the time derivative at the temporal domain grid points 𝑘(=2,3,...,𝑁_𝑇) in the governing equation as𝑑𝑇_𝑖𝑗𝑑𝑡||||_{𝑡=𝑡𝑘}𝑑₂=^𝑁𝑡_𝑚=2𝐴^𝑡𝑘𝑚𝑇_𝑖𝑗𝑚+𝐴^𝑡𝑘1𝑇_𝑖𝑗1,𝑇_𝑖𝑗𝑑𝑡²||||_{𝑡=𝑡𝑘}=^𝑁𝑡_𝑛=2^{𝑁𝑡𝑚=2}+𝐴_𝑡𝑘1𝐴^𝑡𝑘𝑚𝐴^𝑡𝑚𝑛𝑇_𝑖𝑗𝑛,𝑑𝑇_𝑖𝑗𝑑𝑡||||_𝑡=𝑡1+^𝑁𝑡_𝑚=2𝐴^𝑡𝑘𝑚𝐴^𝑡𝑚1𝑇^𝑖𝑗1,(9) where 𝑇_𝑖𝑗1 and 𝑑𝑇_𝑖𝑗/𝑑𝑡|_𝑡=𝑡1 are the initial values of the temperature and its first-order derivative at the beginning of an arbitrary temporal subdomain.

Based on the DQ discretization rules, the DQ discretized form of the governing (4) can be written as 𝜕𝜕𝑥𝑘_𝑖𝑗𝑘^{𝑁𝑥𝑚=1}𝐴^𝑥𝑖𝑚𝑇_𝑚𝑗𝑘+𝜕𝜕𝑘𝑦_𝑖𝑗𝑘^{𝑁𝑦𝑛=1}𝐴^𝑦𝑗𝑛𝑇_𝑖𝑛𝑘+𝑘_𝑖𝑗𝑘^𝑁𝑥_𝑚=1𝐵^𝑥𝑖𝑚𝑇_𝑚𝑗𝑘+^𝑁𝑦_𝑛=1𝐵^𝑦𝑗𝑛𝑇_𝑖𝑛𝑘+𝑔^𝑖𝑗𝑘+𝜏𝜕𝑔𝜕𝑡_𝑖𝑗𝑘=𝜏𝜌𝐶_𝑝^{𝑁𝑡𝑛=2}^{𝑁𝑡𝑚=2}𝐴^𝑡𝑘𝑚𝐴^𝑡𝑚𝑛𝑇_𝑖𝑗𝑛+𝐴^𝑡𝑘1𝑑𝑇_𝑖𝑗𝑑𝑡||||_𝑡=𝑡1+𝜌𝐶_𝑝+^𝑁𝑡_𝑛=2𝐴^𝑡𝑘𝑛𝐴^𝑡𝑛1𝑇^𝑖𝑗1^𝑁𝑡_𝑚=2𝐴^𝑡𝑘𝑚𝑇_𝑖𝑗𝑚+𝐴^𝑡𝑘1𝑇_𝑖𝑗1.(10)

In a similar manner, the boundary conditions can be discretized. An iterative algorithm based on the Newton-Raphson method is employed to solve the resulting nonlinear algebraic system of equations in each temporal subdomain to obtain the temperature distribution at the spatial grid points and at each time step.

3. Numerical Results

3.1. One-Dimensional Problems

In order to demonstrate the accuracy of the method, a non-Fourier one-dimensional transient heat conduction problem subjected to two different types of initial conditions for which the exact solution was obtained by Wang [6] is considered. The boundary and initial conditions are taken as follows.Bou𝑇nd(a0r,y𝑡)c=on0d,itio𝑇ns:𝐿_𝑥Initia⎧⎪⎨⎪⎩lTcyopnedi(ti1o)n:s𝑇:(𝑥,0),=𝑡0,=0𝜕,𝑇(𝑥,0)Type(2):𝑇(𝑥,0)=sin𝜋𝑥𝜕𝑡=1,𝐿_𝑥,𝜕𝑇𝜕(𝑥𝑡,0)=0.(11) Also, in preparation of the numerical results, the values of the relaxation time and thermal diffusivity are assumed to be unit (𝜏=1,𝛼=1). The results for the time history of the temperature at the center of the domain (𝜉=0.5) and also the temperature distribution along the domain at the nondimensional time 𝑡^∗=0.5 are presented in Figures 1 and 2, respectively. Excellent agreement between the results of the presented IDQM and the exact solution [6] is obvious.

(a) Type (1)

(b) Type (2)

Figure 1 

The time history of the nondimensional temperature for both types of the initial conditions at 𝜉=0.5.

(a) Type (1)

(b) Type (2)

Figure 2 

The nondimensional temperature distribution for both types of the initial conditions along the 𝜉-axis at 𝑡∗=0.5.

As another example, a comparison study between the results of the presented method and the finite element method (FEM) for a problem with specified, time-varying temperature at one of its boundary investigated by Dorao [26] is performed. The boundary and initial conditions are taken as follows: 𝑇𝑇(0(,𝑥𝑡,)0=)=𝛽(0𝑡,),𝜕𝑇𝑇(𝑥𝐿,_𝑥0,)𝑡=0,𝜕𝑡=0,(12) where,𝛽(𝑡)=⎧⎪⎨⎪⎩112+342𝑡₀𝑡−1−412𝑡₀𝑡−1³iiff𝑡0≥≤𝑡𝑡₀≤𝑡₀,.(13) Also, the same parameters values as in Dorao [26] are adopted here; that is, the heat capacity 𝛾=1, the relaxation time 𝜏=1, the thermal conductivity 𝑘=1, and 𝑡₀=0.1. By using the FEM, the linear element in conjunction with Galerkin scheme are employed to discretize the spatial and the temporal domain, respectively. The results of the both methods are compared in Figures 3(a) and 3(b) for the time history of the nondimensional temperature at 𝜉=0.5 and the nondimensional temperature along the 𝜁-axis at 𝑡^∗=0.5, respectively. It can be seen that by increasing the FEM number of degrees of freedom (NDOF), its results approaches to those of IDQM. It should also be noticed that the behaviors of the solutions are similar to those of Dorao [26].

(a) The time history of the nondimensional temperature at 𝜉=0.5

(b) Temperature along the 𝜉-axis for different methods at 𝑡∗=0.5

Figure 3 

The convergence of the FEM. Case (1): FEM, (NDOF, CPU time (Second))= (169,0.3); Case (2): FEM, (5329,1); Case (3): FEM, (56943,6); Case (4): IDQM, (8591,1.4).

3.2. Two-Dimensional Problems

In this section, the non-Fourier thermal behaviors of a two-dimensional rectangular domain with temperature-dependent thermal conductivity are investigated using the IDQM. The domain are subjected to the following boundary and initial conditions:−𝑘𝜕𝑇(0,𝑦,𝑡)𝜕𝑥=𝑞^𝑜𝑇(𝑥,0,𝑡)=0,𝛽(𝑡)𝑇,𝑥,𝑇𝐿_𝑦𝐿_𝑥,𝑦,𝑡=0,𝑇(𝑥,𝑦,0)=0,𝜕𝑇(𝑥,𝑦,,0𝑡)=0,𝜕𝑡=0.(14) Otherwise specified, the following values are used in the numerical computations: Fo=1,Ve=0.5,𝑟=1,𝑇_𝑜=1,𝜆=0.1.(15)

As a first example, the convergence of the method against the number of DQ grid points in the spatial and temporal domains (𝑁_𝑥,𝑁_𝑦,𝑁_𝑇,and𝑁_𝐼𝑇) is investigated. For this purpose, the time history of the nondimensional temperature at the point (𝜉=0,𝜂=0.5) and also the nondimensional temperature along the axis 𝜂=0.5 and at the nondimensional time 𝑡^∗=1 are presented in Figures 4(a) and 4(b), respectively. The fast rate of convergence of the method is quite evident.

(a) The time history of the nondimensional temperature at 𝜂=0.5

(b) The nondimensional temperature along the axis 𝜂=0.5 and 𝑡∗=1

Figure 4 

The convergence behavior of IDQM against the number of DQ grid points.

After showing the convergence of the method, the parametric studies are performed.

The time histories of the nondimensional temperature at the different points on the centerline (𝜂=0.5) are shown in Figure 5. The propagation of heat wave can be seen in this figure.

Figure 5 

The time history of the nondimensional temperature at different points along the line 𝜂=0.5.

The time history of the nondimensional temperature at point (𝜉=0,𝜂=0.5) and for the different values of Fourier number (Fo) is shown in Figure 6. It can be seen that as Fourier number or consistently total time increases, the nondimensional temperature approaches to its steady-state value.

Figure 6 

The time history of the nondimensional temperature at the point (𝜉=0, 𝜂=0.5) for different value of Fourier number.

The effect of relaxation time 𝜏, as appeared in Vernotte number (Ve), on the time history of the nondimensional temperature of the point (𝜉=0, 𝜂=0.5) is shown in Figure 7. It is clear that by increasing the Vernotte number, the non-Fourier effect increases. It should be mentioned when Ve=0, the results become equal to those obtained using Fourier law.

Figure 7 

The effect of relaxation time on the time history of the nondimensional temperature at (𝜉=0, 𝜂=0.5).

Figures 8(a) and 8(b) show the time history of the nondimensional temperature at point (𝜁=0, 𝜂=0.5) and nondimensional temperature along the axis 𝜂=0.5 for different value of the nondimensional boundary heat flux Φ, respectively. It can be seen that as the boundary heat flux increases, the temperature and temperature gradient increase.

(a) The time history of the nondimensional temperature at (𝜉=0, 𝜂=0.5)

(b) 𝜂=0.5

Figure 8 

The effects of heat supply at the boundary on the temperature distribution of the domain.

The effect of temperature dependence of the thermal conductivity on the nondimensional temperature distribution of the domain is studied by showing the variation of the nondimensional temperature distribution along 𝜂=0.5 and 𝜁=0 axes in Figures 9(a) and 9(b), respectively. The results are prepared for three different values of the thermal conductivity coefficient 𝜆. Finally, the effect of the aspect ratio (𝑟) on the temperature distribution in the domain is shown in Figure 10. One can see that by increasing the aspect ratio, the temperature at the domain grid point increases.

(a) 𝜂=0.5

(b) 𝜉=0

Figure 9 

The effect of temperature dependence of the thermal conductivity on the nondimensional temperature distribution along 𝜂=0.5 and 𝜉=0 axes.

(a) 𝜂=0.5

(b) 𝜉=0

Figure 10 

The influence of the aspect ratio (𝑟) on the temperature distribution of the domain.

4. Conclusions

As a first attempt, the differential quadrature method (DQM) is employed to study the one- and two-dimensional heat wave propagation in a medium with temperature-dependent thermal conductivity. Both the spatial as well as the temporal domain is discretized using the DQM. The superior accuracy with fewer degrees of freedom than conventional finite element method (FEM) of DQM is shown. The convergence and accuracy of the method is demonstrated. The effects of different parameters such as relaxation time, Fourier number, temperature dependency of thermal conductivity and aspect ratio of the domain on the temperature distribution of the medium are studied. It is shown that these parameters have significant effect on the temperature distribution of the domain.

Nomenclature

𝐴𝜁:	The normalized weighting coefficients of the first-order derivative (𝜁=𝑥,𝑦,𝑡)
𝐵𝜁:	The normalized weighting coefficients of the second-order derivative (𝜁=𝑥,𝑦,𝑡)
𝐶:	Thermal wave speed
𝐶𝑝:	Heat capacity
𝐹𝑜:	Fourier number (=𝛼𝑜𝑡total/𝐿2𝑥)
𝑔:	Internal heat generation source
𝑘:	Thermal conductivity
𝑘𝑜:	Thermal conductivity at reference temperature
𝐿𝑥:	Length of the domain along the 𝑥-direction
𝐿𝑦:	Length of the domain along the 𝑦-direction
𝑁𝐼𝑇:	Number of temporal subdomain
𝑁𝑇:	Number of grid points in of a temporal subdomain
𝑁𝑥:	Number of grid points in the 𝑥-direction
𝑁𝑦:	Number of grid points in the 𝑦-direction
𝑞:	Heat flux
𝑡:	Time
𝑡total:	Total time
𝑡∗:	Nondimensional time
𝑇:	Temperature
𝑇𝑜:	Reference temperature
𝑇∗:	Nondimensional temperature (=𝑇/𝑇𝑜)
𝑉𝑒:	Vernotte number (=√𝛼𝑜𝜏/𝐿𝑥).

Greek Symbols

𝛼:	Thermal diffusivity (𝑘/𝜌𝐶𝑝)
𝛼𝑜:	Thermal diffusivity at reference temperature
Φ:	Nondimensional boundary heat flux (=𝑞𝑜𝐿𝑥/𝑘𝑜𝑇𝑜)
𝜂:	Nondimensional coordinate variable along the 𝑦-axis(=𝑦/𝐿𝑦)
𝜆:	Thermal conductivity coefficient
𝜌:	Density
𝜏:	Relaxation time (=𝛼/𝐶2)
𝜉:	Nondimensional coordinate variable along the 𝑥-axis (=𝑥/𝐿𝑥).