J. of Thermal Science Vol.9, No.2

Difference Scheme for Hyperbolic Heat Conduction Equation with Pulsed Heating Boundary Li Ji

Zhang Zhengfang

Liu Dengying

Institute o f Engineering T h e r m o p h y s i c s , Chinese A c a d e m y o f Sciences, Beijing 100080, China

In this paper, the authors have analyzed the second-order hyperbolic differential equation with pulsed heating boundary, and tried to solve this kind of equation with numerical method. After analyzing the error of the difference scheme and comparing the numerical results with the theoretical results, the authors present some examples to show how the thermal wave propagates in materials. By analyzing the calculation results, the conditions to observe the thermal wave phenomena in the laboratory are conferred. Finally the heat transfer in a complex combined structure is calculated with this method. The result is quite different from the calculated result from the parabolic heat conduction equation.

Keywords: second-order hyperbolic differential equation, pulsed heating boundary, numerical method

Introduction It is well known that the basic law of heat conduction is the Fourier law. It has the form

q = - k . V T , and one-dimensional heat conduction differential equation is ~ T - a -~ -2T • The above ~t ~x 2 equations are derived from the hypothesis that the velocity to establish the thermal balance is infinitely great. In the modern heat conduction theory heat transmits in materials in a limited velocity. The factors to affect the velocity are the thermal properties of the materials. In order to describe this problem, the scholars brought forward a parameter: thermal delay time z0, z0 is related to the thermal movement of the particles. From analysis the heat conduction equation becomes q ( r , x 0 ) --- - k • V T • Thus from derivation, the onedimensional heat conduction differential equation has 3 2T O 2T the hyperbolic form 0 T + x0 - a 2 3t 3t 2 0x This equation is also called thermal wave equation. From some researchers' experiments, the heat conduction has the non-Fourier characteristics under micro-scale level. Consideration of the solvability of the second-order hyperbolic differential equation is necessary. Some analytical methods are studied by Received 2000

many scholars. But those methods have some drawbacks, e.g., the questions must be very simple, and the expressions of results are complex and not easy to use. Especially heat transfer in a combined structure can not be solved by analytical methods. In many experimental research works, the calculations were carried out only in a single material with the hyperbolic heat conduction differential equation to compare with the experimental results. No matter what kind of the model is, the one step or the two steps model, no matter whether the space equilibrium is considered or not, the governing equation is hyperbolic form when we study the thermal wave propagation in materials. In this paper, the author tried to solve this kind of equation with numerical method. After analyzing the different pulsed heating boundaries and the stability of the implicit difference scheme, the conditions for the convergence of the numerical result to the theoretical result are attained. And some new calculated results are shown, which are useful for the comparison with the experimental data.

Mathematical Model As shown in Fig. 1, a high energy pulsed heating resource is put on a semi-indefinite plate. The mathematical model for the structure shown in Fig.1 under pulsed heating is:

Li Ji et al.

Difference Scheme for Hyperbolic Heat Conduction Equation with Pulsed Heating Boundary

~

153

Xo

Materiall

q

q

,

x

~>

X

,

>

>

\

Fig.1 one-dimensional heat conduction model 1

Material 1

at

OZT o

7

O2T =a--

ax 2

(1)

T(O, x) = To,

Iat, c3T ~ ,=o.,=o = 0

(2)

Material 3

k oT

u(t) = [ - at + b,

"c, < t < t2 "

(7)

for the O--~qare definite. But the rectangle form has 3t

Oq

-~21,=0 =q+'r°

0 < t < "cp

It is easy to deal with this kind of boundary waveform

Boundary condition:

a~"'

q = F-U(t)

(3)

Here, F is the maximum heat flux and U(t) is time function, 0_
_kxOT Oq ~x x=xO = q + "co (x) ~ -

(4)

The difference scheme for the main equation is implicit: T' - T' ~ T ' - 2T '-~ + T '-2 P " + %(x) " P P dt dt 2 = a(x)

Material 2

In equation (3) the time function for the triangle pulse waveform is

Initial condition:

-

\

Fig.2 one-dimensional heat conduction model 2

Main equation:

OT --+"c

\

L ' - 2 T ; +T"

u(t)=

~q

and ~

0,

t<0

1,

O
0,

"cp < t

(8)

is indefinite at t=O and t='cp. "cp is pulse

duration. There are two methods to solve this kind of problem: (1) Method of approximate treatment - The rectangle wave is regarded as a very steep echelon; (2) Mathematical manipulation - To transform the indefinite partial derivative into a solvable definite derivative as used in reference [3]. Because the second method is not convenient for use and is not usable for complex conditions, the first method is used here. 1.20E+13

(5)

dx 2

8.00E+12

From derivation, we have

ap.,T,=ae.,L +a,T, +ap~_,T,_~+a.,_zT.,_2

(6)

The ADI method is used to solve the difference equations.

Method to deal with boundary conditions Here, two forms of pulsed heating resources are shown in Fig.3. In Fig.3 the waveforms are a triangle pulse waveform and a rectangle waveform respectively.

4.00E+12

O.OOE+O O.OOE+O

2.00E-12 4.00E-12 6.00E-12 Time s

Fig.3 Pulse waveforms

154

Journal of Thermal Science, Vol.9, No.2, 2000

The point on the boundary has the conservative form as follows: (~T + ~t

0:T ~)T ~q ~0 ~-~'---k-~-x )j=l = q + ~ o ~)-"~-

,

4000 00

I

(9) 3000 00

,

I

,

I

,

- -

NonFouner

dt=T011000

- -

NonFourier

dt=Z'0/100

....

NonFourier theoretical result

,o Stability of the difference scheme

It is very difficult to analyze the stability of the difference scheme of the second-order hyperbolic differential equation. Little information is found in any textbooks. In general the implicit difference scheme is stable unconditionally and convergent to the theoretical result. But in this paper it is found that the implicit scheme for the second-order hyperbolic equation is convergent, but there is a deviation to the theoretical results. Through the analysis, the size of the mesh and the time step value is payoff. Time step value In Fig.4 the dot line is a theoretical result under rectangle pulsed heating in model 1 shown in Fig.l, which is calculated from the expression (10). After pulsed heating there is a temperature step change. The change is determined by lim T(0, t)= F° ~ ' z ° • Because the rectangle wave is t-~0 k approximately treated as a very steep echelon, the choice of the time step is a critical factor. From the calculation the smaller the time step is, the more accurate the result is, as shown in Fig.4. But the calculation needs more computer time.

x

T(x,t)

Foaqr~o.%/k= n ( t - - )vx ~

e

"~

2000.00

E 1000 00 -

0.00

'

0.00E+0

I

I

4 00E-12

'

Time

-H(t-t o

/

1,60E-11

Fig.4 Results of different time steps

4 0 0 0 O0

t

' ....

I , I theoretical result

~

I

numerical solutim dx=lE-lOm /

3 0 0 0 O0

~

- -

,

dt=To/lO0

Fourier null edcal sotut=on

~ 2000.00~ Q..

I---

1000.00 -

....... I t s x2 Io((2--~0)-'~xo) +

I1 ......... . . / . u . =

[~o joe

x'._

,oL~t~-~-,o) - 4--~Zo,a~

0.00

(10)

t 2 t-- o 2 X [e_(l_tO)/2to "10(,/(-:---) --7--) + V 2:'ro am'to

v li~-'o ---)×1 [z o

'

1.20E-11

s

I

0.00E+0

400E-12

I

[

8.00E-12

1.20E-11

Time

x

I

8 00E-12

e .......

io(_l(u)~_ ~ 2zo

x 2 )& 4aZ~

here H ( t ) is Heaviside step functiontll. Space step value The grid division is very influential to the calculation results. In the difference scheme of the second-order hyperbolic equation (6), ap,t_2 is negative and

this is against the requests of

the reliability for the numerical results of the difference scheme-[4]. So it is requested that ap,t_l must be bigger

a "L"o than ap,,_2, here-d--~ - > d t 2 . From calculations the above analysis is correct. The result is shown in Fig.5. From comparison among the calculation results shown in Fig.4 and Fig.5, the result is approaching to

1.60E-11

s

Fig.5 Results of different space steps

the theoretical result more exactly when the space step is changed from 1E-9 m to 1E-10 m. The reason is the control volume method is used in the grid division. But there is no such problem in solving the parabolic heat conduction equation for there is no t e r m ap,t_2Tp,t_ 2 in the implicit difference scheme. As same as the time step, it needs more computer time to decrease the space step. Results a n d Analysis

Model 1 (shown in Fig.l): Heat transmits in material(platinum) with a triangle and a rectangle pulsed heating boundary respectively, Fm~=IE13 W I m 2, Zp=6.35E-12 s, z0=6.35E-12 s.

Li Ji et al.

Difference Scheme for Hyperbolic Heat Conduction Equation with Pulsed Heating Boundary

The numerical results of the model 1 with two pulsed heating boundary are shown in Fig.6 and Fig.7 respectively. In Fig.6(a) the temperature of the three neighbor points begin to change at different time, but the curves of the temperature show the boundary heating information - triangle form. From Fig.6(b) it is clear how the thermal wave transmits in the material. In Fig.7 the pulse waveform is rectangle. The temperature variation takes nearly rectangular form. Model 2 (shown in Fig.2): Heat transmits in a combined structure (1-water, 2-platinum, 3-quartz glass) with an inner rectangle pulsed heating, Fm~=IE13W/m 2, ~p=6.35E-12s, ~o,l=lE-6s, Zo,2=6.36E-12s, ~o,3=lE-10s (~0,1 and "Co,3 a r e estimated);

4000 O0 ~

'

I

,

I

155

,

I

I

,

2000 00

E

(D F-

,

rectangle pulse heating

3000 O0 I

P e

I

Point 1 Point 2

4

1000 O0

0 O0 O00E+O

4.00E-t2 8 0 0 E - 1 2

1.20E-11

Time

1,80E-11 2.00E-11

s

i

2500 00

(a) Temperature variation of Three neighbored points 2000.00 3000.00

.o

i

I

p

[

,

I

,

4

Rectangle pulse

150000

10*dt 100*dt 1000.00

P

2000.00

500"dt

F--

t

-

1000*tit

-

1500*dt 500.00

E 0.00 0.00E+0

p

I 1.00E-11

I 200E-11

J

I--

i 3 00E-11

100000

4.00E-11

Time s

(a) Temperature variation of three neighbored points

1

i;/'~

I I

!

0.130-

'

O.OOE+O

I

4.00E-8

'

I

8.00E-8

'

I

1.20E-7

'

I

i

1.60E-7

2,00E-7

Length m I

2500 O0

i

I

,

(b) Thermal wave propagation in material Fig.7 Model 1 under rectangle pulsed heating

2000 O0

Thermal wave transmission

P .e

8_

E

t=-1000*dt t= 1500*dt

1500.00

t=-500*dt t=-100*dt 1000 00

t=10*dt

I-500 00

0,00 0 00E+0

400E-8

800E-8

1 20E-7

1.60E-7

Length m (b) Thermal wave propagation in material Fig.6 Model 1 under triangle pulsed heating

2.00E-7

The results are shown in Fig.8. From the results it can be found that heat transfer in a combined structure is in a more complex mode than in a single material. The reason may be the thermal waves in different materials affect each other. From Fig.8(b) the depth of the thermal wave movement is different in different materials. In quartz glass the thermal wave can be transfered further than in water at the same time while the magnitude of the temperature variation are almost the same. Heat transfer in a combined structure with an inner pulsed heating boundary governed by the parabolic heat conduction equation is also calculated. It is found that under same pulsed heating boundaries the magnitude of the temperature rise governed by the parabolic equation

156

Journal of Thermal Science, Vol.9, No.2, 2000

is much higher than that governed by the hyperbolic equation. This result is different from the calculated result in a single material as shown in Fig.5. The reason may also be the thermal waves in different materials affect each other and decrease the temperature magnitude. This result is important for the analysis of the experimental results. There are many models for the non-Fourier heat conduction research, but this result is firstly presented in this paper.

3000 O0

~

I

I

~

h

heat conduction can be described by the Fourier Law quite accurately t2].

Discussions Through analyzing the hyperbolic heat conduction equation and the third kind of heating boundary, the non-Fourier heat conduction behavior in materials is mainly affected by the properties of the materials and the turning rate of the heat flux. So if we want to observe the wave movement in heat conduction, the large thermal delay time

in glass

In w a t e r

__

- K ~T ~x

bq

-

l

q

~.

E

should be considered

firstly. And then a proper boundary condition is required. From the boundary condition, it can be derived:

in platium 2000. O0

To

+

r o ~-. -

q

In the above equation the left

term is the dimensionless heat flux in the material. In the right hand the second term is the relative magnitude of the thermal wave to the imposed heat flux. If

1000 O0

11)

0q

r o 0 t > or =1, the thermal wave would be observed in 0 O0

q [ 4.00E-11

O.OOE+O

'

I 8.00E-11

Time

'

I 1 20E-10

i

1.60E-10

the experiment. Otherwise the heat conduction will abide by the Fourier Law.

s

Conclusions (a) Temperature variation of Three neighbored points

2000.00

,

,

J

I

i

XPt=4*dx

I

t

Ox=IE-IO t=-1500*d×

1800.00

lO00*dx

.o

300*dx lO0*d×

120Q00

~"

E

800.00

}400.00

0.00

I -4.00E-9

I 000E+O

I 4.00E-9

I

800E-9

Length m (b) Thermal wave propagation in different materials

Fig.8 Model 2 under rectangle pulsed heating The calculations with amplitude and long pulse The calculation results of equations are same. So in

a rectangle pulse of small duration Zp are carried out. the two kinds of governing conventional conditions the

Numerical method for the two-order hyperbolic heat conduction equation is discussed. The parameters, which affect the thermal wave propagation, are discussed. The numerical results are compared between the hyperbolic heat conduction equation and parabolic heat conduction equation. The main results obtained are as follows: 1)As far as the second-order hyperbolic equation is concerned, the implicit scheme is not enough for the convergence and stability. The reliability of the results must be analyzed. 2) The boundary conditions are quite important for the emergency of the thermal wave; the high turn rate of the heat flux added to the boundary surface is encouraged. 3) Under the micro-scale level conditions compared to the thermal delay time

TO

and the depth of the

thermal wave movement, the thermal wave phenomena could be observed; otherwise, it is difficult to observe such phenomena. 4) In a complex combined structure, the thermal wave propagates in a different way compared with that in a single material. The reason is the thermal waves affect each other. This result is firstly calculated and useful to analyze the experimental data.

Li Ji et al.

Difference Scheme for Hyperbolic Heat Conduction Equation with Pulsed Heating Boundary

Acknowledgements

This article is accomplished under the sponsor of Grant No.KJ951-B 1-704, CAS.

[4]

Reference

[5]

[1]

[6]

[2]

[3]

Xiu Jingzhong, "Non-Fourier effect in transient heat conduction induced by a pulsed heating," Proceedings of China Thermophysics Congress, (1997). Li Ji; Zhang Zhengfang, "Transient Temperature of Liquid on Micro Metal Layer Heated by Pulsed Laser," J.of Thermal Science, 18, No.2, (1999). Z h a n g Zhe; Liu D e n g y i n g , " H y p e r b o l i c Heat Propagation in a Spherical Solid Medium under'

[7]

157

Extremely High Heating Rates," 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, USA,(1998). Tao Wenquan, "Numerical Heat Transfer (in Chinese)," Xi'an Jiaotong University Press, (1995). Lu Jinpu; Guan Ye, "Numerical Method of Differential Equation (in Chinese)," Tsinghua Press, (1989). T.Q. Qiu; C. L. Tien, "Femtosecond Laser Heating of Multi-layer Metal - I. Analysis," Int. J. Heat Mass Transfer, 39, No.17(1994). T.Q.Qiu; T. Juhasz; C. Suarez; W. E. Bron and C. L. Tien, "Femtosecond Laser Heating of Multi-layer Metal - II. Experiments," Int. J. Heat Mass Transfer, 39, No.17, (1994).