# Non-Fourier Heat Conduction Analysis with Temperature-Dependent

Authors: H. Rahideh; P. Malekzadeh; M. R. Golbahar Haghighi
Journal: ISRN Mechanical Engineering
 (2011)

---

## 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.


---

## Body

## 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(k) is assumed to be a linear function of temperature [32](1)k(T)=k0[1+λ(T-T0)].
The other material properties of the domain such as density (ρ) and specific heat capacity (Cp) are assumed to be constant.The hyperbolic constitutive equation governing the transient heat transfer may be written as(2)τ∂q∂t+q=-k∇T,
where τ is the relaxation time, which is related to the thermal wave speed and thermal diffusivity as τ=α/C2. Moreover, the energy equation can be written as(3)ρCp∂T∂t+∇⋅q=g(x,y,t).
Using (1) and (2), (3) takes the following form: (4)(∂k∂x∂T∂x+∂k∂y∂T∂y)+k(∂2T∂x2+∂2T∂y2)+g(x,y,t)+τ∂g∂t=τρCp∂2T∂t2+ρCp∂T∂t.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 ofNx and Ny discrete grid points in the x- and y-direction, respectively. Then, at a given grid point (xi,yj), the first- and second-order derivatives of an arbitrary function f(x,y,t) can be approximated as (5)∂f(x,y,t)∂x|(xi,yj)=∑m=1NxAimxf(xm,yj,t)=∑m=1NxAimxfmj(t),∂2f(x,y,t)∂x2|(xi,yj)=∑m=1NxBimxfmj(t),∂f(x,y,t)∂y|(xi,yj)=∑n=1NyAjnyfin(t),∂2f(x,y,t)∂y2|(xi,yj)=∑n=1NyBjnyfin(t),
where i=1,2,…,Nx and j=1,2,…,Ny.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 (ζ=x,y, and t) are thus determined, respectively, as (6)Aijζ={1LζM(ζi)(ζi-ζj)M(ζj)i≠j,-∑k=1,k≠iNζAikζi=j,i,j=1,2,…,Nζ,Bijζ={2[AiiζAijζ-Aijζζi-ζj]i≠j,-∑k=1,k≠iNζBikζi=j,i,j=1,2,…,Nζ,
where,
(7)M(ζi)=∏j=1,i≠jNζ(ζi-ζj).In numerical calculation, Chebyshev-Gauss-Lobatto quadrature points are used to generate the grid points, that is,(8)ζi=Lζ2{1-cos⁡⁡[(i-1)πNζ-1]}i=1,2,…,Nζ,ζ=x,y,t.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 (∂T/∂t) 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 k(=2,3,...,NT) in the governing equation as(9)dTijdt|t=tk=∑m=2NtAkmtTijm+Ak1tTij1,d2Tijdt2|t=tk=∑n=2Nt∑m=2NtAkmtAmntTijn,+Ak1tdTijdt|t=t1+(∑m=2NtAkmtAm1t)Tij1,
where Tij1 and dTij/dt|t=t1 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 (10)(∂k∂x)ijk∑m=1NxAimxTmjk+(∂k∂y)ijk∑n=1NyAjnyTink+kijk(∑m=1NxBimxTmjk+∑n=1NyBjnyTink)+gijk+τ(∂g∂t)ijk=τρCp[∑n=2Nt∑m=2NtAkmtAmntTijn+Ak1tdTijdt|t=t1+(∑n=2NtAkntAn1t)Tij1]+ρCp(∑m=2NtAkmtTijm+Ak1tTij1).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.(11)Boundaryconditions:T(0,t)=0,T(Lx,t)=0,Initialconditions:{Type(1):T(x,0)=0,∂T(x,0)∂t=1,Type(2):T(x,0)=sin(πxLx),∂T(x,0)∂t=0.
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 t*=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.The time history of the nondimensional temperature for both types of the initial conditions atξ=0.5.
(a)
Type (1)(b)
Type (2)The nondimensional temperature distribution for both types of the initial conditions along theξ-axis at t*=0.5.
(a)
Type (1)(b)
Type (2)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: (12)T(0,t)=β(t),T(Lx,t)=0,T(x,0)=0,∂T(x,0)∂t=0,
where,(13)β(t)={12+34(2tt0-1)-14(2tt0-1)3if0≤t≤t0,1ift≥t0.
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 k=1, and t0=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 t*=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].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).
(a)
The time history of the nondimensional temperature atξ=0.5(b)
Temperature along theξ-axis for different methods at t*=0.5


### 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:(14)-k∂T(0,y,t)∂x=qoβ(t),T(Lx,y,t)=0,T(x,0,t)=0,T(x,Ly,t)=0,T(x,y,0)=0,∂T(x,y,0)∂t=0.
Otherwise specified, the following values are used in the numerical computations: (15)Fo=1,Ve=0.5,r=1,To=1,λ=0.1.As a first example, the convergence of the method against the number of DQ grid points in the spatial and temporal domains(Nx,Ny,NT,andNIT) 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 t*=1 are presented in Figures 4(a) and 4(b), respectively. The fast rate of convergence of the method is quite evident.The convergence behavior of IDQM against the number of DQ grid points.
(a)
The time history of the nondimensional temperature atη=0.5(b)
The nondimensional temperature along the axisη=0.5 and t*=1After 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).Figures8(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.The effects of heat supply at the boundary on the temperature distribution of the domain.
(a)
The time history of the nondimensional temperature at (ξ=0, η=0.5)(b)
η=0.5The 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 (r) 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.The effect of temperature dependence of the thermal conductivity on the nondimensional temperature distribution alongη=0.5 and ξ=0 axes.
(a)
η=0.5(b)
ξ=0The influence of the aspect ratio(r) on the temperature distribution of the domain.
(a)
η=0.5(b)
ξ=0


## 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.(11)Boundaryconditions:T(0,t)=0,T(Lx,t)=0,Initialconditions:{Type(1):T(x,0)=0,∂T(x,0)∂t=1,Type(2):T(x,0)=sin(πxLx),∂T(x,0)∂t=0.
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 t*=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.The time history of the nondimensional temperature for both types of the initial conditions atξ=0.5.
(a)
Type (1)(b)
Type (2)The nondimensional temperature distribution for both types of the initial conditions along theξ-axis at t*=0.5.
(a)
Type (1)(b)
Type (2)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: (12)T(0,t)=β(t),T(Lx,t)=0,T(x,0)=0,∂T(x,0)∂t=0,
where,(13)β(t)={12+34(2tt0-1)-14(2tt0-1)3if0≤t≤t0,1ift≥t0.
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 k=1, and t0=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 t*=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].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).
(a)
The time history of the nondimensional temperature atξ=0.5(b)
Temperature along theξ-axis for different methods at t*=0.5


## 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:(14)-k∂T(0,y,t)∂x=qoβ(t),T(Lx,y,t)=0,T(x,0,t)=0,T(x,Ly,t)=0,T(x,y,0)=0,∂T(x,y,0)∂t=0.
Otherwise specified, the following values are used in the numerical computations: (15)Fo=1,Ve=0.5,r=1,To=1,λ=0.1.As a first example, the convergence of the method against the number of DQ grid points in the spatial and temporal domains(Nx,Ny,NT,andNIT) 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 t*=1 are presented in Figures 4(a) and 4(b), respectively. The fast rate of convergence of the method is quite evident.The convergence behavior of IDQM against the number of DQ grid points.
(a)
The time history of the nondimensional temperature atη=0.5(b)
The nondimensional temperature along the axisη=0.5 and t*=1After 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).Figures8(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.The effects of heat supply at the boundary on the temperature distribution of the domain.
(a)
The time history of the nondimensional temperature at (ξ=0, η=0.5)(b)
η=0.5The 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 (r) 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.The effect of temperature dependence of the thermal conductivity on the nondimensional temperature distribution alongη=0.5 and ξ=0 axes.
(a)
η=0.5(b)
ξ=0The influence of the aspect ratio(r) on the temperature distribution of the domain.
(a)
η=0.5(b)
ξ=0


## 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.



---

*Source: 321605-2011-04-26.xml*