Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 4, pp. 837-845, Warsaw 2013

ESTIMATION OF RELAXATION AND THERMALIZATION TIMES
IN MICROSCALE HEAT TRANSFER MODEL

Bohdan Mochnacki
Częstochowa University of Technology, Czestochowa, Poland and

Higher School of Labour Safety Management, Katowice, Poland

e-mail: bohdan.mochnacki@im.pcz.pl

Marek Paruch
Silesian University of Technology, Institute of Computational Mechanics and Engineering, Gliwice, Poland

e-mail: marek.paruch@polsl.pl

The energy equation corresponding to the dual phase lag model (DPLM) results from the
generalized form of the Fourier law, in which the two ‘delay times’ (relaxation and ther-
malization time) are introduced. The DPLM should be used in the case of microscale heat
transfer analysis, in particular when thermal processes are characterized by extremely short
duration (e.g. ultrafast laser pulse), considerable temperature gradients and very small di-
mensions (e.g. thinmetal film). In this paper, the problem of relaxation and thermalization
time identification is discussed, at the same time the heat transfer processes proceeding
in the domain of a thin metal film subjected to a laser beam are analyzed. The solution
presented bases on the application of evolutionary algorithms. The additional information
concerning the transient temperature distribution on a metal film surface is assumed to be
known. At the stage of numerical realization, the finite difference method (FDM) is used.
In the final part of the paper, an example of computations is presented.

Key words: microscale heat transfer, laser heating, dual phase lagmodel, inverse problem

1. Governing equations

The following form of the generalized Fourier law is considered

q(x,t+ τq)=−λ∇T(x,t+τT) (1.1)

where q is the unitary heat flux, λ is the thermal conductivity, ∇T is the temperature gradient,
τq is the relaxation time (the mean time for electrons to change their energy states), τT is the
thermalization time (the mean time required for electrons and lattice to reach equilibrium).

The DPL equation results from considerations concerning the parabolic two-temperature
model (Al-Nimr, 1997; Chen and Beraun, 2001; Majchrzak and Poteralska, 2011; Majchrzak et
al., 2009), in particular the energy equations determining the heat transfer in the electron gas
and the metal lattice are taken into account. In the case of pure metals, the two-temperature
model takes the form

ce(Te)
∂Te
∂t
=∇[λe(Te)∇Te]−G(Te−Tl) cl(Tl)

∂Tl
∂t
=G(Te−Tl) (1.2)

where Te =Te(x,t),Tl =Tl(x,t) are the temperatures of electrons and lattice, ce(Te), cl(Tl) are
the volumetric specific heats (thermal capacities), λe(Te), λl(Tl) are the thermal conductivities,
G is the coupling factor – this parameter characterizes the energy exchange between a phonon
and electrons (Lin and Zhigilei, 2008).



838 B.Mochnacki,M. Paruch

In Fig. 1 (see Majchrzak and Poteralska, 2011) an example of numerical solution obtained
using the two temperature parabolic model is shown. The metal film (Ti) subjected to a laser
pulse has been considered and the heating (cooling) curves correspond to the film surface. Dif-
ferences between Te and Tl are visible, the time of both temperatures equalization corresponds
to the thermalization one.

Fig. 1. Changes of surface temperatures

Using a certain elimination technique, equations (1.2) canbe substitutedbya single equation
containing the second derivative of temperature with respect to time and a higher-order mixed
derivative in both time and space. Let us assume that the volumetric specific heats ce and cl
are constant values (this simplification is acceptable when the changes of temperature during
the heating process are not very big). So, from equation (1.2)1 it results that

Te =Tl+
cl
G

∂Tl
∂t

(1.3)

Putting (1.2)2 into (1.1), one has

ce
(∂Tl
∂t
+
cl
G

∂2Tl
∂t2

)

=∇(λe∇Tl)+
cl
G
∇

[

λe
∂

∂t
(∇Tl)

]

− cl
∂Tl
∂t

(1.4)

this means

(ce+ cl)
[∂Tl
∂t
+
ce
cl
G(ce+ cl)

∂2Tl
∂t2

]

=∇(λe∇Tl)+
cl
G

∂

∂t
[∇λe(∇Tl)] (1.5)

Denoting

τT =
cl
G

τq =
1

G

( 1

ce
+
1

cl

)

−1

(1.6)

finally one obtains

c
[∂T(x,t)

∂t
+ τq
∂2T(x,t)

∂t2

]

=∇[λ∇T(x,t)]+ τT∇
[

λ
∂∇T(x,t)

∂t

]

(1.7)

where T(x,t) = Tl(x,t) is the macroscopic lattice temperature (Ozisik and Tzou, 1994),
c= cl+ce is the effective volumetric specific heat resulting from the serial assembly of electrons
and phonons and λ=λe (Ozisik and Tzou, 1994).



Estimation of relaxation and thermalization times in ... 839

Another approach to theDPLM formulation is also possible and the followingmathematical
considerations were presented inMajchrzak et al. (2009).

Using the Taylor series expansions, the following first-order approximation of equation (1.1)
can be taken into account

q(x,t)+ τq
∂q(x,t)

∂t
=−λ

[

∇T(x,t)+ τT
∂∇T(x,t)

∂t

]

(1.8)

or

−q(x,t)= τq
∂q(x,t)

∂t
+λ∇T(x,t)+ τTλ

∂∇T(x,t)

∂t
(1.9)

This formula should be introduced to the well knownmacroscopic energy equation

c
∂T(x,t)

∂t
=−∇·q(x,t) (1.10)

This means

c
∂T(x,t)

∂t
= τq
∂

∂t
[∇q(x,t)]+∇[λ∇T(x,t)]+ τT∇

[

λ
∂∇T(x,t)

∂t

]

(1.11)

Substituting −∇q by c(T)∂T/∂t, one obtains the same equation as equation (1.7)

c
[∂T(x,t)

∂t
+ τq
∂2T(x,t)

∂t2

]

=∇[λ∇T(x,t)]+ τT∇
[

λ
∂∇T(x,t)

∂t

]

(1.12)

In this paper, the problem of heat diffusion in the presence of volumetric internal heat sources
Q(x,t) is considered. It can be shown that in this case, equation (1.12) must be supplemented
by additional components, in particular

c
[∂T(x,t)

∂t
+ τq
∂2T(x,t)

∂t2

]

=∇[λ∇T(x,t)]+ τT∇
[

λ
∂∇T(x,t)

∂t

]

+Q(x,t)+ τq
∂Q(x,t)

∂t

(1.13)

It results from the fact that the thermal interactions between external heating (laser pulse) and
the domain of the metal film are taken into account by the introduction of an additional term
supplementing theDPLM, in particular a function corresponding to the volumetric internal heat
sources Q(x,t) is considered. This approach is often used (e.g. Al-Nimr, 1997). The formula
determining the capacity of the internal heat sources takes a form (1D problem – Chen and
Beraun (2001), Chen et al. (2004))

Q(x,t)=

√

µ

π

1−R

tpδ
I0exp

[

−
x

δ
−µ
(t−2tp
tp

)2]

(1.14)

where I0 is the laser intensity which is defined as the total energy carried by a laser pulse per
unit cross-section of the laser beam, tp is the characteristic time of the laser pulse, δ is the
characteristic transparent length of irradiated photons called the absorption depth, R is the
surface reflectivity, µ=4ln2. The local and temporary value of Q results from the distance x
between the surface subjected to laser action and the point considered.Using this approach, the
no-flux boundary conditions for x=0 and x=L should be assumed.

In Fig. 2 the domain considered is shown and its geometrical features allows one to treat the
problem as a 1D one.



840 B.Mochnacki,M. Paruch

Fig. 2. Domain considered

2. Numerical solution basing on fdm (a direct problem)

At the stage of numerical modeling, the finite difference method has been used. The version
proposed constitutes a generalization of the FDM variant proposed by Mochnacki and Suchy
(1995). So, the following energy equation is considered

c
[∂T(x,t)

∂t
+ τq
∂2T(x,t)

∂t2

]

=
∂

∂x

[

λ
∂T(x,t)

∂x

]

+ τT
∂

∂t

∂

∂x

[

λ
∂T(x,t)

∂x

]

+Q(x,t)+ τq
∂Q(x,t)

∂t

(2.1)

The differential mesh is created as a Cartesian product of spatial ∆h and time ∆tmeshes. The
time grid is defined as follows

∆t : t0 <t1 < ... < tf−2 <tf−1 <tf < ... < tF <∞ (2.2)

while the spatial mesh is shown in Fig. 3.

Fig. 3. Themesh

It is visible that the ’boundary’ nodes are located at the distance 0.5h from real boun-
daries (this type of discretization assures a very simple and exact approximation of boundary
conditions, Mochnacki and Suchy (1995)).

The FDM approximation of the spatial differential operator can be taken in the form

∂

∂x

(

λ
∂T

∂x

)f

i
=
T
f
i+1−T

f
i

R
f−1
i+1

Ψi+1+
T
f
i−1−T

f
i

R
f−1
i−1

Ψi−1 (2.3)

where Ψi+1 =Ψi−1 =1/h are the mesh shape functions, while

R
f−1
i+1 =

h

2λ
f−1
i

+
h

2λ
f−1
i+1

R
f−1
i−1 =

h

2λ
f−1
i

+
h

2λ
f−1
i−1

(2.4)

are the thermal resistances between the node i and adjoining the nodes i+1, i−1. The index f
in formula (2.3) shows that the implicit differential scheme will be used here, at the same time,



Estimation of relaxation and thermalization times in ... 841

the thermal conductivities are taken for the time tf−1 to obtain the linear form of final FDM
equations. The FDM approximation of equation (2.1) for transition tf−1 → tf is of the form

c
T
f
i −T

f−1
i

∆t
+ cτq

T
f
i −2T

f−1
i +T

f−2
i

(∆t)2
=
T
f
i+1−T

f
i

R
f−1
i+1

Ψi+1+
T
f
i−1−T

f
i

R
f−1
i−1

Ψi−1

+
τT
∆t

(T
f
i+1−T

f
i

R
f−1
i+1

Ψi+1+
T
f
i−1−T

f
i

R
f−1
i−1

Ψi−1
)

−
τT
∆t

(T
f−1
i+1 −T

f−1
i

R
f−1
i+1

Ψi+1+
T
f−1
i−1 −T

f−1
i

R
f−1
i−1

Ψi−1
)

+Q
f
i + τq

(∂Q

∂t

)f

i

(2.5)

and the last formula can be written as follows (i=1,2, . . . ,N))

AiT
f
i−1+BiT

f
i +CiT

f
i+1 =DiT

f−1
i−1 +EiT

f−1
i +FiT

f−1
i+1 +

τq
(∆t)2

T
f−2
i

−
Q
f
i

c
−
τq
c

(∂Q

∂t

)f

i

(2.6)

where

Ai =
Ψi−1

cR
f−1
i−1

(

1+ τT
∆t

)

Ci =
Ψi+1

cR
f−1
i+1

(

1+ τT
∆t

)

Bi =−
1

∆t

(

1+
τq
∆t

)

−Ai−Ci Di =
Ψi−1

cR
f−1
i−1

τT
∆t

Fi =
Ψi+1

cR
f−1
i+1

τT
∆t

Ei =−
1

∆t

(

1+
2τq
∆t

)

−Di−Fi

(2.7)

Finally

AiT
f
i−1+BiT

f
i +CiT

f
i+1 =G

f
i (2.8)

where

G
f
i =DiT

f−1
i−1 +EiT

f−1
i +FiT

f−1
i+1 +

τq
(∆t)2

T
f−2
i −

Q
f
i

c
−
τq
c

(∂Q

∂t

)f

i
(2.9)

The same equations are accepted for the nodes close to boundaries. It is enough to assume that
the thermal resistances in the directions ’to boundary’ are sufficiently big (e.g. 1010), and then
the non-flux condition is taken into account. The starting point of the numerical simulation
process results from the initial conditions, in particular T0i = T

1
i = T0, i = 1,2, . . . ,N. The

system of FDM equations (2.6) has been solved using the very effective Thomas algorithm for a
three-diagonal linear system of algebraic equations (e.g. Mochnacki and Suchy, 1995).

3. Inverse problem

To solve the inverse problem, the least squares criterion is applied

S(τq,τT)=
1

MF

M
∑

i=1

F
∑

f=1

(T
f
i −T

f
di
)2 (3.1)



842 B.Mochnacki,M. Paruch

where T
f
di
and T

f
i =T(xi, t

f) are themeasured and estimated temperatures, respectively, M is
the number of sensors (e.g. Majchrzak et al., 2008; Mochnacki and Majchrzak, 2007). The
minimum of functional (3.1) has been found using the evolutionary algorithms. So the direct
problems have been solved and the results allow one to determine the time dependent surface
temperature (x=0). Because the temperature history resulting from the numerical solution for
the basic input data is very close to the experimental ones quoted in Chen and Beraun (2001),
Tang and Araki (1999), Fig. 4, therefore this undisturbed numerical solution is assumed to be
the base for the identification problem solution (‘measured surface temperature’). So, the laser
parameters determining capacity of the internal source function Q(x,t) and also the thermal
conductivity and volumetric specific heat of gold are known, the parameters τq, τT should be
determined (from the practical standpoint, the experimental estimation of τq, τT is not easy).

Fig. 4. Comparison with experimental data (Chen and Beraun, 2001)

Fig. 5. (a) Temperature profiles, (b) cooling (heating) curves

InFigs. 5a and5b, an example of the direct problem solution is shown.The layer is subjected
to a short-pulse laser irradiation whose parameters are R=0.93 (reflectivity), I0 =13.7J/m

2

(intensity), tp = 0.1ps = 10
−13 s (time of laser pulse), δ = 15.3nm (absorption depth). The

following parameters of the gold thin film are assumed: thermal conductivity λ=317W/(mK),
volumetric specific heat c=2.4897MJ/(m3K), relaxation time τq =8.5ps, thermalization time
τT = 90ps. The initial temperature equals T0 = 20

◦C (see: Majchrzak et al., 2009; Majchrzak
and Poteralska, 2010).

Using the algorithm presented in the previous Section on the assumption that N = 200
and ∆t=0.005ps the transient temperature field has been found. In Fig. 5a, the temperature
profiles are shown, while Figs. 4 and 5b illustrate the courses of heating (cooling) curves at the
points selected from the domain considered.



Estimation of relaxation and thermalization times in ... 843

The identification of ‘delay’ times has been done using the evolutionary algorithms. In Ta-
ble 1, the algorithmparameters are collected. The results obtained are presented inTable 2, and
they are quite satisfactory.

Table 1.Evolutionary algorithm parameters

Number of
generations

Number of
chromosomes

Prob. of Prob. of Prob. of
Prob. of
cloning

uniform nonuniform arithmetic
mutation mutation crossover

100 75 20% 30% 50% 10%

Table 2.Result of computations using the EA

Design variable Exact value Found value Error [%]

τq 8.5 ·10
−12 8.499999 ·10−12 0

τT 90 ·10
−12 89.99999 ·10−12 0

Figures 6-9 show the process of identification using the evolutionary algorithm after 1st,
10th, 50th and 100th generations. Themark × shows the position of the real value.

Fig. 6. The process of identification – evolutionary algorithm – after 1st generation

Fig. 7. The process of identification – evolutionary algorithm – after 10th generations



844 B.Mochnacki,M. Paruch

Fig. 8. The process of identification – evolutionary algorithm – after 50th generations

Fig. 9. The process of identification – evolutionary algorithm – after 100th generations

4. Final remarks

The application of evolutionary algorithms for solutions of the identification problems is (from
the numerical point of view) time-consuming. On the other hand, however, the mathematical
and numerical problems connectedwith adequate algorithm construction seems to be essentially
simpler in comparison with the very popular gradient methods. The exactness of estimation
parameters τq, τT is very good, and even in the case of real experimental data application the
quality of the obtained results is also satisfactory. Thedetails concerning themethods of indirect
measurements of surface temperatures in the case of microscale domains can be found in Tang
and Araki (1999). The future works will be connected with the problems concerning the micro
models of crystallization (e.g. Lelito et al., 2012).

Acknowledgements

Thisworkhasbeen supportedbyprojectNo. 2012/05/B/ST8/01477sponsoredby thePolishNational

Science Centre for Research andDevelopment.

References

1. Al-NimrM.A., 1997,Heat transfermechanisms during short duration laser heating of thinmetal
films, International Journal of Thermophysics, 18, 5, 1257-1268

2. Chen J.K., Beraun J.E., 2001,Numerical study of ultrashort laser pulse interactionswithmetal
films,Numerical Heat Transfer, Part A, 40, 1-20



Estimation of relaxation and thermalization times in ... 845

3. Chen G., Borca-Tasciuc D., Yang R.G., 2004, Nanoscale heat transfer,Encyclopedia of Na-
noscience and Nanotechnology, X, 1-30

4. Lelito J., Żak P.L., Greer A.L., Suchy J.S., Krajewski W.K., Gracz B., Szucki M.,
Shirzadi A.A., 2012, Crystallization model of magnesium primary phase in the AZ91/SiC com-
posite,Composites Part B: Engineering, 43, 8, 3306-3309

5. Lin Z., Zhigilei L.V., 2008,Electron-phononcoupling and electronheat capacity ofmetals under
conditions of strong electron-phonon nonequilibrium,Physical Review B, 77, 075133-1-075133-17

6. Majchrzak E., Mochnacki B., Greer A.L., Suchy J.S., 2009, Numerical modeling of short
pulse laser interactions with multi-layered thin metal films, CMES: Computer Modeling in Engi-
neering and Sciences, 41, 2, 131-146

7. Majchrzak E.,Mochnacki B., Suchy J.S., 2008, Identification of substitute thermal capacity
of solidifying alloy, Journal of Theoretical and Applied Mechanics, 46, 2, 257-268

8. MajchrzakE.,Poteralska J., 2010,Two-temperaturemicroscale heat treansfermodel. Part 1:
Determination of electrons parameters, Scientific Research of the Institute of Mathematics and
Computer Science, Częstochowa, 9, 1, 99-108

9. Majchrzak E., Poteralska J., 2011, Two temperature model of microscopic heat transfer,
Computer Methods in Material Science, 11, 2, 337-342

10. Mochnacki B., Majchrzak E., 2007, Identification of macro and micro parameters in solidifi-
cationmodel,Bulletin of the Polish Academy of Sciences – Technical Sciences, 55, 1, 107-113

11. Mochnacki B., Suchy J.S., 1995, Numerical Methods in Computations of Foundry Processes,
Polish Foundrymen’s Technical Association, Cracow

12. OzisikM.N.,TzouD.Y., 1994,On thewave theory in heat conduction,Journal of Heat Transfer,
116, 526-535

13. Tang D.W, Araki N., 1999,Wavy, wavelike, diffusive thermal responses of finite rigid slabs to
high-speed heating of laser-pulses, International Journal of Heat and Mass Transfer, 42, 855-860

Manuscript received November 27, 2012; accepted for print January 24, 2013