- Email: [email protected]
Numerical simulation for non-linear thermal wave
Applied Mathematics and Computation 175 (2006) 1385–1399
www.elsevier.com/locate/amc
Numerical simulation for non-linear thermal wave Kuo-Chi Liu Department of Mechanical Engineering, Far East College, 49 Chung Hua Road, Hsin-Shih, Tainan Prefecture 744, Taiwan
Abstract The present work applies the hybrid method of the Laplace transform technique and a modified discretization scheme to analyze the non-linear hyperbolic heat conduction problems in a semi-infinite domain with either linearly or exponentially temperaturedependent thermal conductivity. The Laplace transform method is used to remove the time-dependent terms in the governing differential equation and the boundary conditions, and then the transformed equations are discretized by a modified discretization scheme. To show the rationality and reliability of the present results, a comparison between the present numerical results and those in the literature is made. Results show the behavior of hyperbolic heat conduction is sensitive to the function form of thermal conductivity. The speed of heat propagation changes for the variation of thermal conductivity. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Control volume; Laplace transform method; Modified discretization scheme; Nonlinear
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.08.033
1386
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1. Introduction The finite propagation speed of heat has been demonstrated for extremely short times, for temperatures near absolute zero and for high heat flux [1–3]. Under the circumstances, the heat propagation is a wave-like. Hence, the hyperbolic heat conduction equation was proposed to describe the phenomena of heat transfer at a finite speed. This equation has been a satisfactory extension of the classical heat conduction equation. Due to the mathematics difficulty, the hyperbolic heat conduction problem is always discussed under the assumption of constant thermal properties. For the practical interest, the temperature-dependent effect of thermal properties on heat transfer should be taken into account. The literatures [4–17] studying the hyperbolic heat conduction problem are numerous, but the most have paid attention on the linear problems. Glass et al. [14] used the MacCormackÕs predictor–corrector scheme to solve the HHC problems in a semi-infinite slab with the temperature-dependent thermal conductivity. The results show that the temperature-dependent thermal conductivity strongly affects the temperature distribution. Kar et al. [15] assumed the thermal diffusivity constant and applied the Kirchhoff transformation to linearize the non-linear HHC equation. Chen and Lin [16] have ever proposed the hybrid application of the Laplace transform method and control volume method in conjunction with the TaylorÕs series approximation and the suitable hyperbolic shape functions to solve non-linear HHC problems. Pulvirenti et al. [17] used the MacCormackÕs predictor–corrector scheme to investigate the HHC problem in an infinite solid cylinder with temperature-dependent thermal properties. The previous works [14–17] unanimously assumed that the thermal properties linearly depend on temperature. In practice, it is possible that the thermal properties are not a linear function of temperature. In order to investigate the temperature-dependent effect of thermal conductivity on the temperature distribution of hyperbolic heat conduction, this work considers the thermal conductivity is a linear function or an exponential function of temperature. It can be observed from the previous literatures [4–17] that the major difficulty encountered in the numerical solution of such problems is numerical oscillations in the vicinity of sharp discontinuities. When the thermal conductivity is an exponential function of temperature, the present problem becomes the high non-linear problem. Thus, solving the present problem is a hard task. The present work divides the whole space domain into several sub-space domains, and the thermal conductivity in each sub-space domain is assumed to be constant. Then, the hybrid application of the Laplace transform technique and a modified discretization scheme in conjunction with the suitable hyperbolic shape functions is applied to solve the present problem. When the thermal conductivity inversely changes to temperature, the results in the literatures [14,16] show that the sharp thermal wavefront can be diminished and the speed
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1387
of thermal wave propagation is almost same as that for the constant thermal conductivity. It is an interesting phenomenon and is studied again in the present work.
2. Mathematical formulation In a local energy balance, the energy conservation equation without heat source or sink in the medium can be given as qC P
oT ðx; tÞ oqðx; tÞ þ ¼ 0; ot ox
ð1Þ
where the density q and the specific heat CP are assumed to be constant. In order to accommodate the finite propagation speed of the observed thermal waves, Cattaneo [18] and Vernotte [19] suggested a modified heat flux model in the form s
oqðx; tÞ oT ðx; tÞ þ qðx; tÞ ¼ kðT Þ ; ot ox
ð2Þ
where q is the heat flux, T the temperature, t the time and x is the space coordinate. The thermal conductivity k(T) is assumed to vary with temperature. s denotes relaxation time. In the present study, it is assumed that the thermal conductivity k(T) can be expressed as k(T) = k0k*(T), where k*(T) is a arbitrary function of temperature and k0 is the reference thermal conductivity. Thus, the relaxation time s can be defined as s ¼ qCkP0 c2 , where c0 denotes the reference propagation speed of the 0 thermal wave. Elimination of q between Eqs. (1) and (2) leads to the HHC equation as o2 T oT o oT ¼ aðT Þ s 2 þ ; ð3Þ ot ot ox ox where a is the thermal diffusivity and is defined as a = k(T)/(qCP). For convenience of numerical analysis, the following dimensionless parameters are introduced: g¼
x 1=2
2ðsk 0 =qC P Þ
;
n¼
t ; 2s
Q¼
q ðqC P k 0 T 20 =sÞ1=2
;
and h ¼
T ; T0
ð4Þ
where T0 is the reference temperature. Substitution of Eq. (4) into Eqs. (2) and (3) gives the following dimensionless forms as oQ oh þ 2Q ¼ k ðhÞ on og
ð5Þ
1388
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
and o2 h oh o oh k ðhÞ þ2 ¼ . on og og on2
ð6Þ
In all examples of this work, the dimensionless initial conditions are assumed as hðg; 0Þ ¼ 0;
ohðg; 0Þ ¼ 0; on
and Qðg; nÞ ¼ 0.
ð7Þ
Various types of boundary conditions will be discussed in the following examples.
3. Numerical analysis The present work divides the whole space domain into several sub-space domains, as shown in Fig. 1. The thermal properties in each sub-space domain are assumed to be constant. Thus Eqs. (5) and (6) in the jth sub-space domain can be written as oQj ohj þ 2Qj ¼ k j on og
for gi 6 g 6 giþ1 ;
j¼i
ð8Þ
for gi 6 g 6 giþ1 ;
j ¼ i;
ð9Þ
and 2 o 2 hj ohj o hj ¼ k þ 2 j on og2 on2
where k j is a function with respect to the thermal conductivity of the jth sub h þh space domain and is defined as k j ¼ k i;j 2 iþ1;j . During the operation for computing the values of the dimensionless temperature h, k j is replaced with h þh hiþ1;j are the initial guess values or the previous iterk i;j 2 iþ1;j where hi;j and ation values of hi,j = hj(gi) and hi+1,j = hj(gi+1).
1
1
j-1
i-1
2
j
i
i+1
x Fig. 1. Geometry and coordinates.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1389
The Laplace transform method is employed to map the transient Eqs. (8) and (9) into steady ones. The Laplace transform of a function /(n) with respect to n is defined as follows: Z 1 ~ /ðsÞ ¼ /ðnÞesn dn; ð10Þ 0
where s is the Laplace transform parameter. Eqs. (8) and (9) in the jth subspace domain can be written in the transform domain as ~j ¼ Q
k j d~ hj s þ 2 dg
for gi 6 g 6 giþ1 ;
j¼i
ð11Þ
and hj d2 ~ hj ¼ 0 k2j ~ dg2
for gi 6 g 6 giþ1 ;
j ¼ i;
ð12Þ
where k2j ¼ ðs2 þ 2sÞ=k j :
ð13Þ
Eq. (12) subjected to the boundary conditions ~ hj ðgi Þ ¼ ~ hi;j
and
~ hj ðgiþ1 Þ ¼ ~ hiþ1;j .
ð14Þ
has the analytical solution in the interval [gi, gi+1] as ~ hj ðgÞ ¼
1 ½sinh kj ðgiþ1 gÞ~ hi;j þ sinh kj ðg gi Þ~hiþ1;j . sinh kj ‘
ð15Þ
Similarly, the analytical solution of Eq. (12) in the interval [gi1, gi] is ~ hj1 ðgÞ ¼
1 ½sinh kj1 ðgi gÞ~ hi1;j1 þ sinh kj1 ðg gi1 Þ~hi;j1 ; sinh kj1 ‘ ð16Þ
where ‘ denotes the length of sub-space domain or the distance between two neighboring nodes. In substance, the heat flux and temperature within the whole space domain are continuous, thus the following conditions can be required ~ hj1 ðgi Þ ¼ ~ hj ðgi Þ
ð17Þ
and k j1
d~ hj1 ðgi Þ d~ hj ðgi Þ ¼ k j . dg dg
ð18Þ
Substituting Eqs. (15)–(17) into Eq. (18) can produce the following discretized form for the present problem at the ith node as
1390
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
hi1 þ C i ~hi þ C iþ1 ~ hiþ1 ¼ 0; C i1 ~
i ¼ 2; 3; . . . ; m 1;
ð19Þ
where m is the number of nodes and the coefficients Ci1, Ci, and Ci+1 are given as k j1 kj1 ; sinh kj1 ‘ k j1 kj1 cosh kj1 ‘ k j kj cosh kj ‘ Ci ¼ sinh kj1 ‘ sinh kj ‘
C i1 ¼
ð20aÞ ð20bÞ
and k j kj . sinh kj ‘
C iþ1 ¼
ð20cÞ
Rearrangement of Eq. (19) in conjunction with the discretized form of the boundary conditions yields the following matrix equation as ½Cf~ hg ¼ fF g;
ð21Þ
~ is a column vector repwhere [C] is a matrix with the complex number s, fhg resenting the unknown dimensionless nodal temperatures in the Laplace transform domain, and {F} is a column vector representing the forcing term. The 0 initial values hi are first guessed to calculate [C] and {F}. Thereafter, the application of the Gaussian elimination algorithm and the numerical inversion of the Laplce transform [20] to Eq. (21) can yield the nodal temperatures in the physical domain. The above computation is performed repeatedly until n
nþ1
j hi hi
j < e for i ¼ 1; 2; 3; . . . ; m
and n ¼ 0; 1; . . .
ð22Þ
where the superscript n denotes the nth iteration. At present the tolerance value e is equal to 0.001.
4. Results and discussion This work regards the function form of k*(h) as k*(h) = 1 + bh or k*(h) = exp(bh), where b is the temperature coefficient. To evidence the efficiency of the present numerical scheme for such problems, the results given by Glass et al. [14] and Chen and Lin [16] are illustrated. The present numerical 0 results are performed with ‘ = 0.02 and hi ¼ 1, i = 1, 2, . . . , m. 4.1. Constant boundary temperature The boundary conditions of the first example are considered as hð0; nÞ ¼ 1
ð23aÞ
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1391
and hð1; nÞ ¼ 0.
ð23bÞ
The Laplace transforms of Eqs. (23a) and (23b) with respect to n are can be written as ~ hð0Þ ¼ 1=s
ð24aÞ
~ hð1Þ ¼ 0.
ð24bÞ
and
The present numerical scheme assumes the thermal conductivity constant in each sub-space domain. For conveniently computing the present results, the thermal in the jth sub-space domain, k j , is estimated with conductivity þ h h k i;j 2 iþ1;j . Thus the distance between two neighboring nodes ‘ and the initial 0 guesses of the dimensionless nodal temperatures hi , i = 1, 2, . . . ,m, should be carefully chosen for the present numerical scheme. They may affect the stability and convergence of the present results and the speed of convergence. Fig. 2 displays the present results of h(g, 1) for various k values under k*(h) = 1 + 0.25h 0 and hi ¼ 1 are coincident. This result shows the present results are convergent. 0 0 0 On other hand, hi ¼ 0, hi ¼ 0:5, and hi ¼ 1 are respectively used to estimate the values of h(g, 1) for ‘ = 0.02 and k*(h) = 1 + 0.25h. The differences between
1
θ(η, ξ)
0.8
0.6
0.4
0.2
0 0
k ∗(θ)=1+0.25θ ξ=1.0 =0.01 =0.02 =0.04
0.25
0.5
η
0.75
1
1.25
0 Fig. 2. Effect of ‘ on the values of h (g, 1) for hi ¼ 1 and k*(h) = 1 + 0.25h.
1392
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
the computation results are below 103. It implies that the effect of the initial guesses on the present numerical results is not very significant. The present results are reliable. Fig. 3 shows the dimensionless temperature distribution for the case of k*(h) = 1 + bh at n = 1 and various b values. Glass et al. [14] and Chen and Lin [16] had ever given the results for b = 0.25, 0, and 0.25. Thus a comparison between the present numerical results and the results given by Glass et al. [14] and Chen and Lin [16] is made. It can be observed from Fig. 3 that the present results agree with the results of Glass et al. [14] and Chen and Lin [16], except around the wave front. In accordance with the definition of the thermal wave speed,pffiffiffiffiffiffiffiffiffiffi the ffi thermal wave speed can be described as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ kðT Þ=sqC P ¼ c0 k ðhÞ [14], where c0 is the reference thermal wave speed. Due to the boundary conditions, the temperature distribution in the sample system should be positive in magnitude for t > 0 at the b values, 0.25, 0, and 0.25. Thus, the heat propagation speed will increase with the b value increasing in accordance with the definition of the propagation speed of thermal wave. The thermal wave speed for b = 0.25 should be slower than that for b = 0 at any specific time. This lower velocity is due to the decrease of the diffusion coefficient relative to the linear case [21]. As a result, the penetration distance of thermal wave for b = 0.25 should be different with that for b = 0 at n = 1. However, the results of Glass et al. [14] and Chen and Lin 1.2 k *(θ) = 1+βθ ξ=1.0 Present results Results of Refs.(14)and(16)
1
θ(η, ξ)
0.8 β=2 β=0.25
0.6
β=1
β=0 0.4 β=-0.25 0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
η Fig. 3. Temperature distributions at n = 1 induced by the constant boundary temperature for k*(h) = 1 + bh.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1393
[16] display the thermal wave penetrates the same distance at n = 1 for b = 0.25 and b = 0. It implies that the present numerical results should be rational. The above results further evidence the reliability of the present numerical results. In addition, the degree of non-linearity of the present problem increases with the absolute value of b increasing. Glass et al. [14] only presented the numerical results for b = ±0.25. Chen and Lin [16] further presented the numerical results for b = ±0.5. However, the present numerical scheme can obtain the results of b = 2. This phenomenon shows the present numerical scheme is efficient for such problems. Fig. 4 presents the dimensionless temperature distributions at n = 1 for k*(h) = exp(bh) and various b values. The degree of non-linearity of the present problem increases as the thermal conductivity exponentially depends on temperature. The stability of the present numerical scheme is affected by the higher non-linearity, so that a mildly numerical oscillation appears at the wave front, as shown in Fig. 4. Due to the temperature-dependent effect of thermal conductivity, the temperature distribution curve tends upward for the positive b value and tends downward for the negative b value, as shown in Figs. 3 and 4. In addition, increasing the value of b can promote the ability of heat transfer in the sample slab. The values of the temperature distribution in the sample slab increase with the value of b increasing.
1.2
k *(θ)=exp(βθ ) ξ=1.0
1
0.8
θ(η, ξ)
β=1 0.6
β=0.25 0.4
β=0.15
β=-0.25 0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
η Fig. 4. Temperature distributions at n = 1 induced by the constant boundary temperature for k*(h) = exp(bh).
1394
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399 1.2 k *(θ)=exp(θ/2) k *(θ)=1+θ/2
1
θ(η, ξ)
0.8
0.6
ξ=0.1 ξ=0.4
0.4
ξ=0.6
0.2
0 0
0.2
0.4
η
0.6
0.8
1
Fig. 5. Temperature distributions at various n values induced by the constant boundary temperature for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h).
Fig. 5 shows a comparison between the temperature distributions for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at various early time values. It can be observed from Fig. 5 that the locations of the jump discontinuities for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at a specific dimensionless time are different. It further displays the propagation speed of thermal wave is dependent of the functional form of k*(h). At the same time, there is a deviation between the temperature distributions for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at a specific dimensionless time. The values of the temperature distributions for k*(h) = exp(0.5h) are bigger than those for k*(h) = 1 + 0.5h. The rapidly varying k*(h) has a significant effect on the profile of h(g, n). In other words, the temperature distribution is sensitive to the functional form of k*(h). 4.2. Pulse boundary heating A pulse source heats on the boundary surface g = 0, thus the boundary conditions are described as Qð0; nÞ ¼ uðnÞ uðn DnÞ
ð25aÞ
hð1; nÞ ¼ 0;
ð25bÞ
and
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1395
where u is a step function and is defined as 1 if Dn P n; u½n Dn ¼ 0 if Dn < n.
ð26Þ
The Laplace transforms of Eqs. (25a) and (25b) with respect to n are can be written as d~ hð0Þ 1 ¼ ð1 þ 2s Þ½1 expðDnsÞ dg k1
ð27Þ
hð1; nÞ ¼ 0.
ð28Þ
and Figs. 6 and 7 demonstrate the temperature distributions at n = 0.5 for k*(h) = 1 + bh and k*(h) = exp(bh), respectively. The major difficulty encountered in the numerical solution of the linear HHC problem is numerical oscillation in the vicinity of sharp discontinuity [4–13]. The pulse disturbance, as described by Eq. (25a), causes two sharp discontinuities in both edges of the pulse wave, so that the possibility of occurring numerical oscillation increases. It can be found from Figs. 6 and 7 that no obviously numerical oscillation appears. This phenomenon implies that the present method can suppress the numerical oscillation and solve the strongly non-linear HHC problems. 1 k *(θ) = 1+βθ ξ=0.5
0.8
θ(η, ξ)
0.6
0.4 β=-0.25 0.2
0
β=0 β=0.25
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 6. Temperature distributions at n = 0.5 induced by a pulse boundary heating of Dn = 0.1 for k*(h) = 1 + bh.
1396
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399 1 k *(θ) = exp(βθ) ξ=0.5
0.8
θ(η, ξ)
0.6
0.4
β=-0.25 0.2
β=0
0
β=0.25
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 7. Temperature distributions at n = 0.5 induced by a pulse boundary heating of Dn = 0.1 for k*(h) = exp(bh).
Because that the thermal conductivity apparently changes in the region of thermal wave for b = 0.25 and 0.25, the reflection phenomenon is created at the left edge of the thermal wave. The thermal conductivity in the region of the thermal wave increases for k*(h) = 1 + 0.25h and k*(h) = exp(0.25h). The higher thermal conductivity makes the heat transfer ability increasing [12]. The more energy moving through the left edge of the thermal wave into the region of the thermal wave. Thus it is observed that the reflected wave is downward for b = 0.25. Conversely, the reflected wave is upward for b = 0.25. Fig. 8 displays the temperature distributions induced by a heating pulse of Dn = 0.1 for k*(h) = 1 and k*(h) = exp(0.5h) at various n values. It is clear that the reflection phenomenon at the left edge of the thermal wave does not happen for k*(h) = 1. It also can be found from Fig. 8 that the reflection phenomenon gradually decays with the n value increasing for k*(h) = exp(0.5h), because the strength of pulse disturbance gradually decays with the n value increasing so that the value of k*(h) = exp(0.5h) approaches to 1. As h > 0, the value of k*(h) = exp(0.5h) is greater than k*(h) = 1. Thus the wave speed for k*(h) = exp(0.5h) is greater than that for k*(h) = 1, and the width of the pulse disturbance expands. The effect of the heating duration Dn on the temperature distribution is studied. The temperature distributions plotted in Fig. 9 are for n = 0.5 and k*(h) = 1 + 0.5h. It is found from Fig. 9 that the leading edges of the pulse for various values of Dn, 0.02, 0.04, and 0.06, are coincident. It implies that
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1397
1.2 k*(θ) = exp(0.5θ) k*(θ) =1 Δξ=0.1
ξ=0.1 ξ=0.4
θ(η, ξ)
0.8
ξ=0.8
0.4
0 0
0.2
0.4
η
0.6
0.8
1
Fig. 8. Temperature distributions at various n values induced by a pulse boundary heating of Dn = 0.1 for k*(h) = 1 and k*(h) = exp(0.5h).
0.8 k *(θ) =1+0.5θ ξ=0.5
0.6 Δξ=0.06 Δξ=0.04
0.4
θ(η, ξ)
Δξ=0.02 0.2
0
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 9. Effect of the heating duration Dn on the temperature distribution at n = 0.5 for k*(h) = 1 + 0.5h.
1398
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
the wave speed is independent of the heating duration Dn. Due to an increase in energy heated on the boundary surface g = 0, the values of the temperature distribution are on the rise for the value of Dn increasing. 5. Conclusion A numerical scheme is developed for the non-linear HHC problem. Numerical results for the temperature distributions induced by a pulse boundary disturbance and by a constant boundary temperature are presented. The degree of non-linearity of the present problems increases with the absolute value of b increasing. The present numerical scheme can obtain the results of the high non-linearity cases without the serious numerical oscillation. This shows the efficiency of the present numerical scheme for such problems. The dependent-temperature thermal conductivity affects the temperature distribution and the thermal wave speed. The behavior of hyperbolic heat conduction is sensitive to the functional form of k*(h). References [1] D.C. Kelly, Diffusion: a relativistic appraisal, Am. J. Phys. 36 (1968) 585–591. [2] A. Vedavarz, K. Mitra, S. Kumar, Hyperbolic temperature profile for laser surface interactions, J. Appl. Phys. 76 (1994) 5014–5021. [3] A. Vedavarz, S. Kumar, M.K. Moallemi, Significance of non-Fourier heat waves in conduction, J. Heat Transfer 116 (1994) 221–223. [4] K.J. Baumeister, T.D. Hamill, Hyperbolic heat-conduction equation—a solution for the semiinfinite body problem, ASME J. Heat Transfer 93 (1971) 126–127. [5] M.J. Mourer, H.A. Thompson, Non-Fourier effects at high heat flux, ASME J. Heat Transfer 95 (1973) 284–286. [6] D.C. Wiggert, Analysis of early-time heat conduction by method of characteristics, J. Heat Transfer 99 (1977) 35–40. [7] G.F. Carey, M. Tsai, Hyperbolic heat transfer with reflection, Numer. Heat Transfer 5 (1982) 309–327. [8] C. Bai, A.S. Lavine, Hyperbolic heat conduction in a superconducting film, ASME. Thermal Eng. Proc. 49 (1991) 87–92. [9] H.T. Chen, J.Y. Lin, Numerical analysis for hyperbolic heat conduction, Int. J. Heat Mass Transfer 36 (1993) 2891–2898. [10] J.Y. Lin, H.T. Chen, Numerical solution of hyperbolic heat conduction in cylindrical and spherical systems, Appl. Math. Model. 18 (1994) 384–390. [11] D.W. Tang, N. Araki, Wavy, wavelike, diffusive thermal responses of finite rigid slabs to highspeed heating of laser-pulses, Int. J. Heat Mass Transfer 42 (1999) 855–860. [12] H.T. Chen, K.C. Liu, Study of hyperbolic heat conduction problem in the film and substrate composite with the interface resistance, Jan. J. Appl. Phys. 41 (2002) 6267–6275. [13] K.C. Liu, H.T. Chen, Numerical analysis for the hyperbolic heat conduction problem under a pulsed surface disturbance, J. Appl. Math. Comput. 159 (2004) 887–901. ¨ zisik, D.S. McRae, Hyperbolic heat conduction with temperature[14] D.E. Glass, M.N. O dependent thermal conductivity, J. Appl. Phys. 59 (1986) 1861–1865.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1399
[15] A. Kar, C.L. Chan, J. Mazumder, Comparative studies on nonlinear hyperbolic and parabolic heat conduction for various boundary conditions: analytic and numerical solutions, J. Heat Transfer 114 (1992) 14–20. [16] J.Y. Lin, H.T. Chen, Study of hyperbolic heat conduction with temperature-dependent thermal properties, J. Heat Transfer 116 (1994) 750–753. [17] B. Pulvirenti, A. Barletta, E. Zanchini, Finite-difference solution of hyperbolic conduction with temperature-dependent properties, Numer. Heat Transfer A 34 (1998) 169–183. [18] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, Comptes Rendus 247 (1958) 431–433. [19] P. Vernotte, Les paradoxes de la the´orie continue de lÕequation de la chaleur, Comptes Rendus 246 (1958) 3154–3155. [20] G. Honig, U. Hirdes, A method for the numerical inversion of Laplace transforms, J. Comp. Appl. Math. 10 (1984) 113–132. [21] V. Me´ndez, Nonlinear hyperbolic heat conduction, J. Non-equilibrium Thermodyn. 22 (1997) 217–232.
www.elsevier.com/locate/amc
Numerical simulation for non-linear thermal wave Kuo-Chi Liu Department of Mechanical Engineering, Far East College, 49 Chung Hua Road, Hsin-Shih, Tainan Prefecture 744, Taiwan
Abstract The present work applies the hybrid method of the Laplace transform technique and a modified discretization scheme to analyze the non-linear hyperbolic heat conduction problems in a semi-infinite domain with either linearly or exponentially temperaturedependent thermal conductivity. The Laplace transform method is used to remove the time-dependent terms in the governing differential equation and the boundary conditions, and then the transformed equations are discretized by a modified discretization scheme. To show the rationality and reliability of the present results, a comparison between the present numerical results and those in the literature is made. Results show the behavior of hyperbolic heat conduction is sensitive to the function form of thermal conductivity. The speed of heat propagation changes for the variation of thermal conductivity. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Control volume; Laplace transform method; Modified discretization scheme; Nonlinear
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.08.033
1386
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1. Introduction The finite propagation speed of heat has been demonstrated for extremely short times, for temperatures near absolute zero and for high heat flux [1–3]. Under the circumstances, the heat propagation is a wave-like. Hence, the hyperbolic heat conduction equation was proposed to describe the phenomena of heat transfer at a finite speed. This equation has been a satisfactory extension of the classical heat conduction equation. Due to the mathematics difficulty, the hyperbolic heat conduction problem is always discussed under the assumption of constant thermal properties. For the practical interest, the temperature-dependent effect of thermal properties on heat transfer should be taken into account. The literatures [4–17] studying the hyperbolic heat conduction problem are numerous, but the most have paid attention on the linear problems. Glass et al. [14] used the MacCormackÕs predictor–corrector scheme to solve the HHC problems in a semi-infinite slab with the temperature-dependent thermal conductivity. The results show that the temperature-dependent thermal conductivity strongly affects the temperature distribution. Kar et al. [15] assumed the thermal diffusivity constant and applied the Kirchhoff transformation to linearize the non-linear HHC equation. Chen and Lin [16] have ever proposed the hybrid application of the Laplace transform method and control volume method in conjunction with the TaylorÕs series approximation and the suitable hyperbolic shape functions to solve non-linear HHC problems. Pulvirenti et al. [17] used the MacCormackÕs predictor–corrector scheme to investigate the HHC problem in an infinite solid cylinder with temperature-dependent thermal properties. The previous works [14–17] unanimously assumed that the thermal properties linearly depend on temperature. In practice, it is possible that the thermal properties are not a linear function of temperature. In order to investigate the temperature-dependent effect of thermal conductivity on the temperature distribution of hyperbolic heat conduction, this work considers the thermal conductivity is a linear function or an exponential function of temperature. It can be observed from the previous literatures [4–17] that the major difficulty encountered in the numerical solution of such problems is numerical oscillations in the vicinity of sharp discontinuities. When the thermal conductivity is an exponential function of temperature, the present problem becomes the high non-linear problem. Thus, solving the present problem is a hard task. The present work divides the whole space domain into several sub-space domains, and the thermal conductivity in each sub-space domain is assumed to be constant. Then, the hybrid application of the Laplace transform technique and a modified discretization scheme in conjunction with the suitable hyperbolic shape functions is applied to solve the present problem. When the thermal conductivity inversely changes to temperature, the results in the literatures [14,16] show that the sharp thermal wavefront can be diminished and the speed
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1387
of thermal wave propagation is almost same as that for the constant thermal conductivity. It is an interesting phenomenon and is studied again in the present work.
2. Mathematical formulation In a local energy balance, the energy conservation equation without heat source or sink in the medium can be given as qC P
oT ðx; tÞ oqðx; tÞ þ ¼ 0; ot ox
ð1Þ
where the density q and the specific heat CP are assumed to be constant. In order to accommodate the finite propagation speed of the observed thermal waves, Cattaneo [18] and Vernotte [19] suggested a modified heat flux model in the form s
oqðx; tÞ oT ðx; tÞ þ qðx; tÞ ¼ kðT Þ ; ot ox
ð2Þ
where q is the heat flux, T the temperature, t the time and x is the space coordinate. The thermal conductivity k(T) is assumed to vary with temperature. s denotes relaxation time. In the present study, it is assumed that the thermal conductivity k(T) can be expressed as k(T) = k0k*(T), where k*(T) is a arbitrary function of temperature and k0 is the reference thermal conductivity. Thus, the relaxation time s can be defined as s ¼ qCkP0 c2 , where c0 denotes the reference propagation speed of the 0 thermal wave. Elimination of q between Eqs. (1) and (2) leads to the HHC equation as o2 T oT o oT ¼ aðT Þ s 2 þ ; ð3Þ ot ot ox ox where a is the thermal diffusivity and is defined as a = k(T)/(qCP). For convenience of numerical analysis, the following dimensionless parameters are introduced: g¼
x 1=2
2ðsk 0 =qC P Þ
;
n¼
t ; 2s
Q¼
q ðqC P k 0 T 20 =sÞ1=2
;
and h ¼
T ; T0
ð4Þ
where T0 is the reference temperature. Substitution of Eq. (4) into Eqs. (2) and (3) gives the following dimensionless forms as oQ oh þ 2Q ¼ k ðhÞ on og
ð5Þ
1388
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
and o2 h oh o oh k ðhÞ þ2 ¼ . on og og on2
ð6Þ
In all examples of this work, the dimensionless initial conditions are assumed as hðg; 0Þ ¼ 0;
ohðg; 0Þ ¼ 0; on
and Qðg; nÞ ¼ 0.
ð7Þ
Various types of boundary conditions will be discussed in the following examples.
3. Numerical analysis The present work divides the whole space domain into several sub-space domains, as shown in Fig. 1. The thermal properties in each sub-space domain are assumed to be constant. Thus Eqs. (5) and (6) in the jth sub-space domain can be written as oQj ohj þ 2Qj ¼ k j on og
for gi 6 g 6 giþ1 ;
j¼i
ð8Þ
for gi 6 g 6 giþ1 ;
j ¼ i;
ð9Þ
and 2 o 2 hj ohj o hj ¼ k þ 2 j on og2 on2
where k j is a function with respect to the thermal conductivity of the jth sub h þh space domain and is defined as k j ¼ k i;j 2 iþ1;j . During the operation for computing the values of the dimensionless temperature h, k j is replaced with h þh hiþ1;j are the initial guess values or the previous iterk i;j 2 iþ1;j where hi;j and ation values of hi,j = hj(gi) and hi+1,j = hj(gi+1).
1
1
j-1
i-1
2
j
i
i+1
x Fig. 1. Geometry and coordinates.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1389
The Laplace transform method is employed to map the transient Eqs. (8) and (9) into steady ones. The Laplace transform of a function /(n) with respect to n is defined as follows: Z 1 ~ /ðsÞ ¼ /ðnÞesn dn; ð10Þ 0
where s is the Laplace transform parameter. Eqs. (8) and (9) in the jth subspace domain can be written in the transform domain as ~j ¼ Q
k j d~ hj s þ 2 dg
for gi 6 g 6 giþ1 ;
j¼i
ð11Þ
and hj d2 ~ hj ¼ 0 k2j ~ dg2
for gi 6 g 6 giþ1 ;
j ¼ i;
ð12Þ
where k2j ¼ ðs2 þ 2sÞ=k j :
ð13Þ
Eq. (12) subjected to the boundary conditions ~ hj ðgi Þ ¼ ~ hi;j
and
~ hj ðgiþ1 Þ ¼ ~ hiþ1;j .
ð14Þ
has the analytical solution in the interval [gi, gi+1] as ~ hj ðgÞ ¼
1 ½sinh kj ðgiþ1 gÞ~ hi;j þ sinh kj ðg gi Þ~hiþ1;j . sinh kj ‘
ð15Þ
Similarly, the analytical solution of Eq. (12) in the interval [gi1, gi] is ~ hj1 ðgÞ ¼
1 ½sinh kj1 ðgi gÞ~ hi1;j1 þ sinh kj1 ðg gi1 Þ~hi;j1 ; sinh kj1 ‘ ð16Þ
where ‘ denotes the length of sub-space domain or the distance between two neighboring nodes. In substance, the heat flux and temperature within the whole space domain are continuous, thus the following conditions can be required ~ hj1 ðgi Þ ¼ ~ hj ðgi Þ
ð17Þ
and k j1
d~ hj1 ðgi Þ d~ hj ðgi Þ ¼ k j . dg dg
ð18Þ
Substituting Eqs. (15)–(17) into Eq. (18) can produce the following discretized form for the present problem at the ith node as
1390
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
hi1 þ C i ~hi þ C iþ1 ~ hiþ1 ¼ 0; C i1 ~
i ¼ 2; 3; . . . ; m 1;
ð19Þ
where m is the number of nodes and the coefficients Ci1, Ci, and Ci+1 are given as k j1 kj1 ; sinh kj1 ‘ k j1 kj1 cosh kj1 ‘ k j kj cosh kj ‘ Ci ¼ sinh kj1 ‘ sinh kj ‘
C i1 ¼
ð20aÞ ð20bÞ
and k j kj . sinh kj ‘
C iþ1 ¼
ð20cÞ
Rearrangement of Eq. (19) in conjunction with the discretized form of the boundary conditions yields the following matrix equation as ½Cf~ hg ¼ fF g;
ð21Þ
~ is a column vector repwhere [C] is a matrix with the complex number s, fhg resenting the unknown dimensionless nodal temperatures in the Laplace transform domain, and {F} is a column vector representing the forcing term. The 0 initial values hi are first guessed to calculate [C] and {F}. Thereafter, the application of the Gaussian elimination algorithm and the numerical inversion of the Laplce transform [20] to Eq. (21) can yield the nodal temperatures in the physical domain. The above computation is performed repeatedly until n
nþ1
j hi hi
j < e for i ¼ 1; 2; 3; . . . ; m
and n ¼ 0; 1; . . .
ð22Þ
where the superscript n denotes the nth iteration. At present the tolerance value e is equal to 0.001.
4. Results and discussion This work regards the function form of k*(h) as k*(h) = 1 + bh or k*(h) = exp(bh), where b is the temperature coefficient. To evidence the efficiency of the present numerical scheme for such problems, the results given by Glass et al. [14] and Chen and Lin [16] are illustrated. The present numerical 0 results are performed with ‘ = 0.02 and hi ¼ 1, i = 1, 2, . . . , m. 4.1. Constant boundary temperature The boundary conditions of the first example are considered as hð0; nÞ ¼ 1
ð23aÞ
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1391
and hð1; nÞ ¼ 0.
ð23bÞ
The Laplace transforms of Eqs. (23a) and (23b) with respect to n are can be written as ~ hð0Þ ¼ 1=s
ð24aÞ
~ hð1Þ ¼ 0.
ð24bÞ
and
The present numerical scheme assumes the thermal conductivity constant in each sub-space domain. For conveniently computing the present results, the thermal in the jth sub-space domain, k j , is estimated with conductivity þ h h k i;j 2 iþ1;j . Thus the distance between two neighboring nodes ‘ and the initial 0 guesses of the dimensionless nodal temperatures hi , i = 1, 2, . . . ,m, should be carefully chosen for the present numerical scheme. They may affect the stability and convergence of the present results and the speed of convergence. Fig. 2 displays the present results of h(g, 1) for various k values under k*(h) = 1 + 0.25h 0 and hi ¼ 1 are coincident. This result shows the present results are convergent. 0 0 0 On other hand, hi ¼ 0, hi ¼ 0:5, and hi ¼ 1 are respectively used to estimate the values of h(g, 1) for ‘ = 0.02 and k*(h) = 1 + 0.25h. The differences between
1
θ(η, ξ)
0.8
0.6
0.4
0.2
0 0
k ∗(θ)=1+0.25θ ξ=1.0 =0.01 =0.02 =0.04
0.25
0.5
η
0.75
1
1.25
0 Fig. 2. Effect of ‘ on the values of h (g, 1) for hi ¼ 1 and k*(h) = 1 + 0.25h.
1392
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
the computation results are below 103. It implies that the effect of the initial guesses on the present numerical results is not very significant. The present results are reliable. Fig. 3 shows the dimensionless temperature distribution for the case of k*(h) = 1 + bh at n = 1 and various b values. Glass et al. [14] and Chen and Lin [16] had ever given the results for b = 0.25, 0, and 0.25. Thus a comparison between the present numerical results and the results given by Glass et al. [14] and Chen and Lin [16] is made. It can be observed from Fig. 3 that the present results agree with the results of Glass et al. [14] and Chen and Lin [16], except around the wave front. In accordance with the definition of the thermal wave speed,pffiffiffiffiffiffiffiffiffiffi the ffi thermal wave speed can be described as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ kðT Þ=sqC P ¼ c0 k ðhÞ [14], where c0 is the reference thermal wave speed. Due to the boundary conditions, the temperature distribution in the sample system should be positive in magnitude for t > 0 at the b values, 0.25, 0, and 0.25. Thus, the heat propagation speed will increase with the b value increasing in accordance with the definition of the propagation speed of thermal wave. The thermal wave speed for b = 0.25 should be slower than that for b = 0 at any specific time. This lower velocity is due to the decrease of the diffusion coefficient relative to the linear case [21]. As a result, the penetration distance of thermal wave for b = 0.25 should be different with that for b = 0 at n = 1. However, the results of Glass et al. [14] and Chen and Lin 1.2 k *(θ) = 1+βθ ξ=1.0 Present results Results of Refs.(14)and(16)
1
θ(η, ξ)
0.8 β=2 β=0.25
0.6
β=1
β=0 0.4 β=-0.25 0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
η Fig. 3. Temperature distributions at n = 1 induced by the constant boundary temperature for k*(h) = 1 + bh.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1393
[16] display the thermal wave penetrates the same distance at n = 1 for b = 0.25 and b = 0. It implies that the present numerical results should be rational. The above results further evidence the reliability of the present numerical results. In addition, the degree of non-linearity of the present problem increases with the absolute value of b increasing. Glass et al. [14] only presented the numerical results for b = ±0.25. Chen and Lin [16] further presented the numerical results for b = ±0.5. However, the present numerical scheme can obtain the results of b = 2. This phenomenon shows the present numerical scheme is efficient for such problems. Fig. 4 presents the dimensionless temperature distributions at n = 1 for k*(h) = exp(bh) and various b values. The degree of non-linearity of the present problem increases as the thermal conductivity exponentially depends on temperature. The stability of the present numerical scheme is affected by the higher non-linearity, so that a mildly numerical oscillation appears at the wave front, as shown in Fig. 4. Due to the temperature-dependent effect of thermal conductivity, the temperature distribution curve tends upward for the positive b value and tends downward for the negative b value, as shown in Figs. 3 and 4. In addition, increasing the value of b can promote the ability of heat transfer in the sample slab. The values of the temperature distribution in the sample slab increase with the value of b increasing.
1.2
k *(θ)=exp(βθ ) ξ=1.0
1
0.8
θ(η, ξ)
β=1 0.6
β=0.25 0.4
β=0.15
β=-0.25 0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
η Fig. 4. Temperature distributions at n = 1 induced by the constant boundary temperature for k*(h) = exp(bh).
1394
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399 1.2 k *(θ)=exp(θ/2) k *(θ)=1+θ/2
1
θ(η, ξ)
0.8
0.6
ξ=0.1 ξ=0.4
0.4
ξ=0.6
0.2
0 0
0.2
0.4
η
0.6
0.8
1
Fig. 5. Temperature distributions at various n values induced by the constant boundary temperature for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h).
Fig. 5 shows a comparison between the temperature distributions for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at various early time values. It can be observed from Fig. 5 that the locations of the jump discontinuities for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at a specific dimensionless time are different. It further displays the propagation speed of thermal wave is dependent of the functional form of k*(h). At the same time, there is a deviation between the temperature distributions for k*(h) = 1 + 0.5h and k*(h) = exp(0.5h) at a specific dimensionless time. The values of the temperature distributions for k*(h) = exp(0.5h) are bigger than those for k*(h) = 1 + 0.5h. The rapidly varying k*(h) has a significant effect on the profile of h(g, n). In other words, the temperature distribution is sensitive to the functional form of k*(h). 4.2. Pulse boundary heating A pulse source heats on the boundary surface g = 0, thus the boundary conditions are described as Qð0; nÞ ¼ uðnÞ uðn DnÞ
ð25aÞ
hð1; nÞ ¼ 0;
ð25bÞ
and
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1395
where u is a step function and is defined as 1 if Dn P n; u½n Dn ¼ 0 if Dn < n.
ð26Þ
The Laplace transforms of Eqs. (25a) and (25b) with respect to n are can be written as d~ hð0Þ 1 ¼ ð1 þ 2s Þ½1 expðDnsÞ dg k1
ð27Þ
hð1; nÞ ¼ 0.
ð28Þ
and Figs. 6 and 7 demonstrate the temperature distributions at n = 0.5 for k*(h) = 1 + bh and k*(h) = exp(bh), respectively. The major difficulty encountered in the numerical solution of the linear HHC problem is numerical oscillation in the vicinity of sharp discontinuity [4–13]. The pulse disturbance, as described by Eq. (25a), causes two sharp discontinuities in both edges of the pulse wave, so that the possibility of occurring numerical oscillation increases. It can be found from Figs. 6 and 7 that no obviously numerical oscillation appears. This phenomenon implies that the present method can suppress the numerical oscillation and solve the strongly non-linear HHC problems. 1 k *(θ) = 1+βθ ξ=0.5
0.8
θ(η, ξ)
0.6
0.4 β=-0.25 0.2
0
β=0 β=0.25
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 6. Temperature distributions at n = 0.5 induced by a pulse boundary heating of Dn = 0.1 for k*(h) = 1 + bh.
1396
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399 1 k *(θ) = exp(βθ) ξ=0.5
0.8
θ(η, ξ)
0.6
0.4
β=-0.25 0.2
β=0
0
β=0.25
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 7. Temperature distributions at n = 0.5 induced by a pulse boundary heating of Dn = 0.1 for k*(h) = exp(bh).
Because that the thermal conductivity apparently changes in the region of thermal wave for b = 0.25 and 0.25, the reflection phenomenon is created at the left edge of the thermal wave. The thermal conductivity in the region of the thermal wave increases for k*(h) = 1 + 0.25h and k*(h) = exp(0.25h). The higher thermal conductivity makes the heat transfer ability increasing [12]. The more energy moving through the left edge of the thermal wave into the region of the thermal wave. Thus it is observed that the reflected wave is downward for b = 0.25. Conversely, the reflected wave is upward for b = 0.25. Fig. 8 displays the temperature distributions induced by a heating pulse of Dn = 0.1 for k*(h) = 1 and k*(h) = exp(0.5h) at various n values. It is clear that the reflection phenomenon at the left edge of the thermal wave does not happen for k*(h) = 1. It also can be found from Fig. 8 that the reflection phenomenon gradually decays with the n value increasing for k*(h) = exp(0.5h), because the strength of pulse disturbance gradually decays with the n value increasing so that the value of k*(h) = exp(0.5h) approaches to 1. As h > 0, the value of k*(h) = exp(0.5h) is greater than k*(h) = 1. Thus the wave speed for k*(h) = exp(0.5h) is greater than that for k*(h) = 1, and the width of the pulse disturbance expands. The effect of the heating duration Dn on the temperature distribution is studied. The temperature distributions plotted in Fig. 9 are for n = 0.5 and k*(h) = 1 + 0.5h. It is found from Fig. 9 that the leading edges of the pulse for various values of Dn, 0.02, 0.04, and 0.06, are coincident. It implies that
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1397
1.2 k*(θ) = exp(0.5θ) k*(θ) =1 Δξ=0.1
ξ=0.1 ξ=0.4
θ(η, ξ)
0.8
ξ=0.8
0.4
0 0
0.2
0.4
η
0.6
0.8
1
Fig. 8. Temperature distributions at various n values induced by a pulse boundary heating of Dn = 0.1 for k*(h) = 1 and k*(h) = exp(0.5h).
0.8 k *(θ) =1+0.5θ ξ=0.5
0.6 Δξ=0.06 Δξ=0.04
0.4
θ(η, ξ)
Δξ=0.02 0.2
0
-0.2 0
0.2
0.4
η
0.6
0.8
1
Fig. 9. Effect of the heating duration Dn on the temperature distribution at n = 0.5 for k*(h) = 1 + 0.5h.
1398
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
the wave speed is independent of the heating duration Dn. Due to an increase in energy heated on the boundary surface g = 0, the values of the temperature distribution are on the rise for the value of Dn increasing. 5. Conclusion A numerical scheme is developed for the non-linear HHC problem. Numerical results for the temperature distributions induced by a pulse boundary disturbance and by a constant boundary temperature are presented. The degree of non-linearity of the present problems increases with the absolute value of b increasing. The present numerical scheme can obtain the results of the high non-linearity cases without the serious numerical oscillation. This shows the efficiency of the present numerical scheme for such problems. The dependent-temperature thermal conductivity affects the temperature distribution and the thermal wave speed. The behavior of hyperbolic heat conduction is sensitive to the functional form of k*(h). References [1] D.C. Kelly, Diffusion: a relativistic appraisal, Am. J. Phys. 36 (1968) 585–591. [2] A. Vedavarz, K. Mitra, S. Kumar, Hyperbolic temperature profile for laser surface interactions, J. Appl. Phys. 76 (1994) 5014–5021. [3] A. Vedavarz, S. Kumar, M.K. Moallemi, Significance of non-Fourier heat waves in conduction, J. Heat Transfer 116 (1994) 221–223. [4] K.J. Baumeister, T.D. Hamill, Hyperbolic heat-conduction equation—a solution for the semiinfinite body problem, ASME J. Heat Transfer 93 (1971) 126–127. [5] M.J. Mourer, H.A. Thompson, Non-Fourier effects at high heat flux, ASME J. Heat Transfer 95 (1973) 284–286. [6] D.C. Wiggert, Analysis of early-time heat conduction by method of characteristics, J. Heat Transfer 99 (1977) 35–40. [7] G.F. Carey, M. Tsai, Hyperbolic heat transfer with reflection, Numer. Heat Transfer 5 (1982) 309–327. [8] C. Bai, A.S. Lavine, Hyperbolic heat conduction in a superconducting film, ASME. Thermal Eng. Proc. 49 (1991) 87–92. [9] H.T. Chen, J.Y. Lin, Numerical analysis for hyperbolic heat conduction, Int. J. Heat Mass Transfer 36 (1993) 2891–2898. [10] J.Y. Lin, H.T. Chen, Numerical solution of hyperbolic heat conduction in cylindrical and spherical systems, Appl. Math. Model. 18 (1994) 384–390. [11] D.W. Tang, N. Araki, Wavy, wavelike, diffusive thermal responses of finite rigid slabs to highspeed heating of laser-pulses, Int. J. Heat Mass Transfer 42 (1999) 855–860. [12] H.T. Chen, K.C. Liu, Study of hyperbolic heat conduction problem in the film and substrate composite with the interface resistance, Jan. J. Appl. Phys. 41 (2002) 6267–6275. [13] K.C. Liu, H.T. Chen, Numerical analysis for the hyperbolic heat conduction problem under a pulsed surface disturbance, J. Appl. Math. Comput. 159 (2004) 887–901. ¨ zisik, D.S. McRae, Hyperbolic heat conduction with temperature[14] D.E. Glass, M.N. O dependent thermal conductivity, J. Appl. Phys. 59 (1986) 1861–1865.
K.-C. Liu / Appl. Math. Comput. 175 (2006) 1385–1399
1399
[15] A. Kar, C.L. Chan, J. Mazumder, Comparative studies on nonlinear hyperbolic and parabolic heat conduction for various boundary conditions: analytic and numerical solutions, J. Heat Transfer 114 (1992) 14–20. [16] J.Y. Lin, H.T. Chen, Study of hyperbolic heat conduction with temperature-dependent thermal properties, J. Heat Transfer 116 (1994) 750–753. [17] B. Pulvirenti, A. Barletta, E. Zanchini, Finite-difference solution of hyperbolic conduction with temperature-dependent properties, Numer. Heat Transfer A 34 (1998) 169–183. [18] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, Comptes Rendus 247 (1958) 431–433. [19] P. Vernotte, Les paradoxes de la the´orie continue de lÕequation de la chaleur, Comptes Rendus 246 (1958) 3154–3155. [20] G. Honig, U. Hirdes, A method for the numerical inversion of Laplace transforms, J. Comp. Appl. Math. 10 (1984) 113–132. [21] V. Me´ndez, Nonlinear hyperbolic heat conduction, J. Non-equilibrium Thermodyn. 22 (1997) 217–232.