Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

47, 2, pp. 383-396, Warsaw 2009

NUMERICAL SIMULATION OF THERMAL PROCESSES

PROCEEDING IN A MULTI-LAYERED FILM SUBJECTED TO

ULTRAFAST LASER HEATING

Ewa Majchrzak
Silesian University of Technology, Gliwice, Poland

e-mail: ewa.majchrzak@polsl.pl

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

e-mail: moch@imi.pcz.pl

Józef S. Suchy
AGH, University of Science and Technology, Cracow, Poland

e-mail: jsuchy@agh.edu.pl

In the paper, the mathematical model, numerical algorithm and examples of
computations concerning thermal processes proceeding in a multi-layered thin
film subjected to an ultrafast laser pulse are discussed. The equations descri-
bing a course of the analysed process correspond to the dual-phase-lag model
and contain both the relaxation time τq and additionally the thermalization
time τT . At the stage of numerical simulation, the finite differencemethod has
been used. The algorithm is based on an artificial decomposition of the domain
considered, while common thermal interactions between successive layers are
taken into account using conditions of heat flux and temperature continuity at
points corresponding to internal boundaries (1D task has been considered).

Keywords:microscaleheat transfer, thinfilms, laserpulse, numericalmodelling

1. Introduction

Classical Fourier’s equation constitutes a quite good mathematical descrip-
tion of heat conduction processes proceeding in macro domains subjected to
external thermal interactions whose duration is not very short, at the same
time the temperature considered T(x,t) should be essentially bigger than 0K
(Al-Nimr, 1997; Chen et al., 2004; Escobar et al., 2006; Özişik andTzou, 1994;
TammaandZhou, 1998). It iswell known that theFourier law assumes instan-



384 E. Majchrzak et al.

taneous heat propagation and this assumption leads to evident errors when
the time considered is comparable with the relaxation time τq of heat carriers
(average time between successive electron-phonon collisions in the conductors
or semiconductors and phonon-phonon collisions in dielectrics). Additionally,
the Fourier equation is acceptable when the characteristic dimension L of the
domain considered is essentially larger than themean free path Λ of the heat
carriers (the average distance that energy carriers travel between successive
collisions). So, generally speaking, considering the processes proceeding in do-
mains for which L ¬ Λ (e.g. thin films) subjected to ultrafast heating (e.g.
short-pulse laser interaction) othermodels of heat transport phenomenamust
be taken into account. The limitations concerning the Fourier model applica-
tions are discussed, among others, in Escobar et al. (2006), Özişik and Tzou
(1994), Tamma and Zhou (1998).

In the paper, the problem of heat transfer proceeding in a multi-layered
thin film subjected to a short pulse laser heating is considered. It should be
pointed out that the heat transport through thin films is of vital importance
in microtechnology applications (Dai and Nassar, 2001a,b; Smith and Nor-
ris, 2003). Themathematical description of the process discussed is based on
the dual-phase-lag model presented, among others, in Escobar et al. (2006),
Özişik and Tzou (1994), Tamma and Zhou (1998), Tzou (1995). Taking into
account characteristic features of thin film geometry one can assume that the
components of heat flux in macro-directions (e.g. x2, x3) result from the tra-
ditional Fourier law, while in the definition of heat flux in the direction x1
the relaxation time τq and additionally the thermalization time τT (themean
time required for electrons and lattice to reach equilibrium) are introduced
(Dai and Nassar, 2001). So, in this direction the heat flux and temperature
gradient will occur at different times.

Themathematical model presented in the next section concerns a 1D pro-
blemcorrespondingto themicro-direction x1 =x (heat fluxes in thedirections
x2, x3 are neglected). For most short laser pulse interactions with thin films,
the laser spot size is much larger than the film thickness. Therefore, it is re-
asonable to treat the interactions as a one-dimensional heat transfer process
(Chen and Beraun, 2001).

At the stage of numerical realisation, an algorithm based on the finite
difference method is applied, at the same time a certain concept of domain
decomposition is proposed. Temporary temperature fields in successive layers
are calculated separately, while the continuity conditions allow one to find
the temporary solution concerning the entire domain. In the final part of the
paper, examples of computations are shown.



Numerical simulation of thermal processes... 385

2. Heat transport at the microscale

Heat transport equations describing thermalbehaviour of a thinfilm, as shown
in Fig.1, can be written in the form (Dai and Nassar, 2000; Dai and Nassar,
2001; Özişik and Tzou, 1994; Tamma and Zhou, 1998; Tzou, 1995)

c
∂T(x,t)

∂t
=−∇·q(x,t)+Q(x,t) (2.1)

and

q1(x,t+ τq)=−λ
∂T(x,t+ τT)

∂x1
(2.2)

q2(x,t)=−λ
∂T(x,t)

∂x2
q3(x,t) =−λ

∂T(x,t)

∂x3

where x = {x1,x2,x3}, q = [q1,q2,q3]
⊤ is the heat flux, λ is the thermal

conductivity, c is the volumetric specific heat, Q is the capacity of internal
heat sources, τT , τq are the positive constants which are the time lags of the
temperature gradient and heat flux, respectively.

Fig. 1. Thin film

Using the Taylor series expansion, the following first-order approximation
of equation (2.2)1 can be taken into account

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

∂t
=−λ

[∂T(x,t)

∂x1
+ τT

∂

∂t

(∂T(x,t)

∂x1

)]

(2.3)

Equation (2.1), which in the case of 1D problem (x=x1) is of the form

c
∂T(x,t)

∂t
=−
∂q(x,t)

∂x
+Q(x,t) (2.4)



386 E. Majchrzak et al.

where

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

∂t
=−λ

[∂T(x,t)

∂x
+ τT

∂

∂t

(∂T(x,t)

∂x

)]

(2.5)

should be supplemented by adequate boundary and initial conditions.

3. Multi-layered domain

Let us consider amulti-layered thin film of thickness L=L1+L2+ . . .+LM
(as in Fig.2) with the initial temperature distribution T(x,0)=T0, constant
thermal properties of successive layers, ideal thermal contact between the lay-
ers and insulated external boundaries. The front surface x = 0 is irradiated
by a laser pulse whose output intensity equals I(t). According to Tang and
Araki (1999), the conductional heat transfer can bemodeled by equation (2.4)
with internal volumetric heat sources Q(x,t), at the same time for x=0 the
non-flux condition can be assumed. In this paper, the following formula (Kaba
and Dai, 2005; Tzou and Chiu, 2001) has been applied

Q(x,t)=

√

β

π

1−R

tpδ
I0exp

(

−
x

δ
−
√

β
|t−2tp|

tp

)

(3.1)

where I0 is the laser intensity which is defined as the total energy carried by
the 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 reflectivity of the irradiated
surface and β=4ln2 (Chen and Beraun, 2001).
The local and temporary value of Q(x,t) results from the distance x be-

tween the surface subjected to laser action and the point considered. So, the
following system of equations is taken into account

x∈Ωm : cm
∂Tm(x,t)

∂t
=−
∂qm(x,t)

∂x
+Qm(x,t) m=1,2, . . . ,M

(3.2)

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

∂t
=−λm

[∂Tm(x,t)

∂x
+ τTm

∂

∂t

(∂Tm(x,t)

∂x

)]

Theboundary conditions on the contact surfaces between the sub-domains
have the form of continuity ones, which means

x∈Γm :

{

Tm(x,t)=Tm+1(x,t)

qm(x,t)= qm+1(x,t)
m=1,2, . . . ,M−1 (3.3)



Numerical simulation of thermal processes... 387

Fig. 2. Multi-layered domain

The initial conditions are assumed in the following way

t=0 : Tm(x,0)=Tm0
∂Tm(x,t)

∂t

∣

∣

∣

∣

t=0

=0 (3.4)

4. Numerical model

At the stage of numerical computations, the finite differencemethod has been
used, while the final system of equations has been solved using the Thomas
algorithm (Majchrzak andMochnacki, 2004;Mochnacki and Suchy, 1995) (se-
parately for successive layers).
To find a numerical solution to the problem discussed, a staggered grid is

introduced (Dai and Nassar, 2000), as shown in Fig.3. For convenience, we

omit m and denote T
f
i = T(ih,f∆t), where h is the mesh size, ∆t is the

time step, i = 0,2,4, . . . ,N, f = 0,1, . . . ,F, and q
f
j = q(jh,f∆t), where

j=1,3, . . . ,N−1.

Fig. 3. Discretization

Aswasmentioned, the numerical procedure proposed is based on theTho-
mas algorithm for a tridiagonal linear system of equations and decomposition



388 E. Majchrzak et al.

of the domain into M sub-domains corresponding to successive layers. Ad-
ditionally, an adequate procedure of contact temperatures computations is
introduced.

Let us consider an internal point xi ∈Ωm. The finite difference approxi-
mation of equation (3.2)1 can be written as follows (implicit scheme)

ci
T
f
i −T

f−1
i

∆t
=−
q
f
i+1− q

f
i−1

2h
+Qi (4.1)

where the index i corresponds to ’temperature nodes’ (Fig.3) belonging to
the Ωm sub-domain.

Equation (3.2)2 can be transformed to the form

q
f
j+τqj

q
f
j −q

f−1
j

∆t
=−λj

(T
f
j+1−T

f
j−1

2h

)

−
λjτTj
∆t

(T
f
j+1−T

f
j−1

2h
−
T
f−1
j+1 −T

f−1
j−1

2h

)

(4.2)
or

q
f
j =

τqj

∆t
(

1+
τqj
∆t

)q
f−1
j −

λj
(

1+
τTj
∆t

)

2h
(

1+
τqj
∆t

)(T
f
j+1−T

f
j−1)+

(4.3)

+
λjτTj

2h∆t
(

1+
τqj
∆t

)(T
f−1
j+1 −T

f−1
j−1 )

where the index j corresponds to ’heat flux nodes’ (Fig.3) belonging to the
Ωm sub-domain.

The last equation allows one to construct similar formulas for the nodes
i−1, i+1, and then one obtains (τqi = τqi−1 = τqi+1, τTi = τTi−1 = τTi+1,
λi =λi−1 =λi+1, of course)

q
f
i−1−q

f
i+1=

τqi

∆t
(

1+
τqi
∆t

)(q
f−1
i−1 −q

f−1
i+1 )+

λi
(

1+ τTi
∆t

)

2h
(

1+
τqi
∆t

)(T
f
i−2−2T

f
i +T

f
i+2)+

(4.4)

−
λiτTi

2h∆t
(

1+
τqi
∆t

)(T
f−1
i−2 −2T

f−1
i +T

f−1
i+2 )



Numerical simulation of thermal processes... 389

Putting (4.4) into (4.1), one has

ci
T
f
i −T

f−1
i

∆t
=

τqi

2h∆t
(

1+
τqi
∆t

)(q
f−1
i−1 −q

f−1
i+1 )+

+
λi
(

1+ τTi
∆t

)

4h2
(

1+
τqi
∆t

)(T
f
i−2−2T

f
i +T

f
i+2)+ (4.5)

−
λiτTi

4h2∆t
(

1+
τqi
∆t

)(T
f−1
i−2 −2T

f−1
i +T

f−1
i+2 )+Q

f
i

or
AiT

f
i−2− (1+2Ai)T

f
i +AiT

f
i+2 =D

f
i (4.6)

where

Ai =
λi∆t

(

1+ τTi
∆t

)

4h2ci
(

1+
τqi
∆t

) (4.7)

and

D
f
i =BiT

f−1
i−2 − (1+2Bi)T

f−1
i +BiT

f−1
i+2 +Ci(q

f−1
i+1 −q

f−1
i−1 )−

∆t

ci
Q
f
i (4.8)

while

Bi =
λiτTi

4h2ci
(

1+
τqi
∆t

) Ci =
τqi

2hci
(

1+
τqi
∆t

) (4.9)

Let us assume that the contact temperatures T
f
i = T

f
cm at the boundary

points xm,m=1,2, . . . ,M−1 are known. Then the temperature field at the
time tf results from the following systems of equations:
— first layer

T
f
0 −T

f
2 =

τT1
∆t

1+ τT1
∆t

(T
f−1
2 −T

f−1
0 )

AiT
f
i−2− (1+2Ai)T

f
i +AiT

f
i+2 =D

f
i i=2,4, . . . ,N1−2 (4.10)

T
f
N1
=T

f
c1

— internal layers

T
f
Nm−1

=T
f
cm−1

AiT
f
i−2− (1+2Ai)T

f
i +AiT

f
i+2 =D

f
i

i=Nm−1+2,Nm−1+4, . . . ,Nm−2 (4.11)

T
f
Nm
=Tfcm



390 E. Majchrzak et al.

— last layer

T
f
NM−1

=T
f
cM−1

AiT
f
i−2− (1+2Ai)T

f
i +AiT

f
i+2 =D

f
i

i=NM−1+2,NM−1+4, . . . ,N −2 (4.12)

T
f
N−2
−T

f
N
=

τTM
∆t

1+ τTM
∆t

(T
f−1
N
−T

f−1
N−2
)

Finally, the problem of computations of contact temperatures will be expla-
ined.Thecontinuity condition qm(x,t)= qm+1(x,t)= qcm(x,t), formula (3.3),
leads to the equation x∈Γm

τqm
∂qm(x,t)

∂t
+λm

∂Tm(x,t)

∂x
+λmτTm

∂2Tm(x,t)

∂t∂x
=

(4.13)

= τqm+1
∂qm+1(x,t)

∂t
+λm+1

∂Tm+1(x,t)

∂x
+λm+1τTm+1

∂2Tm+1(x,t)

∂t∂x

This formula should be written down using the finite difference convention,
and then

αmT
f
cm =λm

(

1+
τTm
∆t

)

T
f
Nm−2

+λm+1
(

1+
τTm+1
∆t

)

T
f−1
Nm+2

+

+
λmτTm
∆t
(Tf−1cm −T

f−1
Nm−2

)+λmτTm(T
f−1
cm −T

f−1
Nm−2

)+ (4.14)

+λm+1τTm+1(T
f−1
cm −T

f−1
Nm+2

)+
2h

∆t
(τqm+1− τqm)(q

f
cm−q

f−1
cm )

where

αm =λm
(

1+
τTm
∆t

)

+λm+1
(

1+
τTm+1
∆t

)

(4.15)

In the place of qfcm and q
f−1
cm , the arithmetic means of heat flux values at

the points Nm−1, Nm+1 are introduced. The starting point of computations
consists in assumption that T0cm =T

1
cm =T0 and q

0
cm =0.Next, the systemof

equations (4.10), (4.11), (4.12) is solved and the heat fluxes at the odd internal
nodes are found bymeans of equation (4.3). Finally, the contact temperatures
Tfcm are calculated using formula (4.14) and the next loop of computations
can be realised. The method proposed is very quick and effective owing to
application of the Thomas algorithm and decomposition of the domain.



Numerical simulation of thermal processes... 391

5. Results of computations

To test the accuracy and effectiveness of the method proposed, at first the
following task has been solved. The layer of thickness L = 10−4 and ther-
mophysical parameters λ = 1, c = 1, τq = 1/π

2 +100, τT = 1/π
2 +10−6,

Q(x,t)= 0 has been considered. For the data assumed, the problemdescribed
by equations (2.4), (2.5) and boundary-initial conditions in the form

T(0, t)= 0 T(L,t)= 0

T(x,0)= sin(104πx)
∂T(x,t)

∂t

∣

∣

∣

∣

t=0

=−π2 sin(104πx)
(5.1)

has the following analytical solution (Dai and Nassar, 2001)

T(x,t)= exp(−π2t)sin(104πx) (5.2)

So, the domain considered has beendivided in an artificial way into 4parts
of the same thickness, and ideal thermal contact between the sub-domains has
been assumed. Using the algorithm presented in the previous sections on the
assumption that h=5·10−7 and ∆t=0.0001, the transient temperature field
has been found and the results have been compared with the exact solution.
Both solutions are very close as shown in Fig.4.

Fig. 4. Analytical (lines) and numerical (symbols) solutions

The second task is connected with the numerical solution presented inDai
and Nassar (2000) which concerns a two-layer domain (gold layer and chro-



392 E. Majchrzak et al.

mium layer of thicknesses 50nm). In order to test the algorithmdiscussed, the
domain consideredhasbeendivided into 4parts (Ω1 and Ω2 correspond to the
gold sub-domain, Ω3 and Ω4 correspond to the chromium sub-domain). The
layers are subjected to short-pulse laser irradiation (R=0.93, I0 =13.7J/m

2,
tp = 100fs, δ = 15.3nm). Thermophysical parameters of the sub-domains
are the following: λ = 317W/(mK), c = 2.4897MJ/(m3K), τq = 8.5ps
(1ps=10−12s), τT = 90ps (gold), λ = 93W/(mK), c = 3.2148MJ/(m

3K),
τq =0.136ps, τT =7.86ps (chromium). Themesh step: h=1nm, time step:
∆t=0.005ps.

In Fig.5, the temperature profiles (temperature rise above T0 = 27
◦C)

for the instants 0.2ps and 0.25ps are shown. The results of both solutions
are close. The temperatures obtained using the algorithm presented here are
bigger, indeed. It results from the fact that the laser interaction was proba-
bly approximated in a little different way, additionally the approach to the
continuity conditions and the concept of decomposition were different, too.

Fig. 5. Temperature profiles – comparison with solution (symbols) presented in Dai
and Nassar (2000)

The last example concerns the alternating gold-chromium-gold-chromium
layers. Thermophysical parameters of thematerials are the sameas previously,
the laser characteristic is also the same.

In Fig.6 the temperature profiles (temperature rise above T0 =27
◦C) for

1 – 0.4ps, 2 – 0.6ps, 3 – 0.8ps and 4 – 1ps are shown. Figure 7 illustrates
the course of temperature at the surface subjected to laser heating (x = 0)
and the internal surfaces x=L1, x=L1+L2.



Numerical simulation of thermal processes... 393

Fig. 6. Temperature profiles in the multi-layer domain

Fig. 7. Heating (cooling) curves at points selected from the domain Ω

6. Final remarks

The presentedmodel based on the dual-phase-lag approach contains both the
relaxation time τq and additionally the thermalization time τT . In literature
concerning the microscale heat transfer, one can also find models for which
only the relaxation time is taken into account. In this place the well known
Cataneo equation can be mentioned. According to present opinions resulting
mainly from experiments (Özişik and Tzou, 1994; Tank and Araki, 1999), it



394 E. Majchrzak et al.

seems that the assumption concerning a non-zero value of τT gives results
closer to real physical conditions of the microscale heat transfer.
The algorithm presented can be simply generalised for 2D or 3D tasks.

The components determining q2(x,t) and q3(x,t) result then directly from
the classical Fourier law. A numerical solution obtained in this way gives a
possibility to analyse the influence of laser pulse distribution in the direc-
tions x2 and x3 on the course of heating and cooling processes in the domain
considered.
Themodel presented here can beused for analysis of a heat transfer proce-

eding in multi-layered domains being a composition of an optional number of
thin filmswith different parameters. The choice ofmaterials considered in this
paper results, first of all, from the available in literature input data. It seems
thatmore close to real thermal conditions is the 2Dmodel corresponding to an
axially symmetrical domain, and this problem will be a subject of the future
research.

Acknowledgement

The paper is a part of projectMTKD-CT-2006-042468.

References

1. Al-Nimr M.A., 1997, Heat transfer mechanisms during short duration laser
heating of thin metal films, International Journal of Thermophysics, 18, 5,
1257-1268

2. Chen G., Borca-Tasciuc D., Yang R.G., 2004, Nanoscale heat transfer,
In: Encyclopedia of Nanoscience and Nanotechnology, Edited by H.S. Nalwa,
Vol. X

3. Chen J.K., Beraun J.E., 2001, Numerical study of ultrashort laser pulse
interactions with metal films,Numerical Heat Transfer, Part A, 40, 1-20

4. DaiW.,NassarR., 2000,Adomain decompositionmethod for solving three-
dimensional heat transport equations in a double-layered thin filmwithmicro-
scale thickness,Numerical Heat Transfer, Part A, 38, 243-255

5. Dai W., Nassar R., 2001a, A compact finite difference scheme for solving a
one-dimensional heat transport equation at themicroscale, Journal of Compu-
tational and Applied Mathematics, 132, 431-441

6. Dai W., Nassar R., 2001b, A finite difference scheme for solving a three-
dimensional heat transport equation in a thin film with microscale thickness,
International Journal for Numerical Methods in Engineering, 50, 1665-1680



Numerical simulation of thermal processes... 395

7. EscobarR.A.,GhauS.S., JhonM.S.,AmonC.H., 2006,Multi-length and
time scale thermal transport using the lattice Boltzmann method with appli-
cation to electronic cooling, International Journal of Heat and Mass Transfer,
49, 97-107

8. Kaba I.K., Dai W., 2005, A stable three-level finite difference scheme for
solving theparabolic two-stepmodel in a 3Dmicro-sphereheatedbyultrashort-
pulsed lasers,Journal ofComputational andAppliedMathematics,181, 125-147

9. Majchrzak E., Mochnacki B., 2004, Numerical Methods. Theoretical Ba-
ses, Practical Aspects and Algorithms, Publication of the Silesian University of
Technology, Gliwice [in Polish]

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

11. Özişik M.N., Tzou D.Y., 1994, On the wave theory in heat conduction,
Journal of Heat Transfer, 116, 526-535

12. Smith A.N., Norris P.M., 2003, Microscale heat transfer, Chapter 18 in:
Heat Transfer Handbook, JohnWiley and Sons

13. Tamma K.K., Zhou X., 1998, Macroscale and microscale thermal transport
and thermo-mechanical interactions: some noteworthy perspectives, Journal of
Thermal Stresses, 21, 405-449

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

15. Tzou D.Y., 1995, A unified field approach for heat conduction frommicro- to
macro-scales, Journal of Heat Transfer, 117, 8-16

16. Tzou D.Y., Chiu K.S., 2001, Temperature-dependent thermal lagging in ul-
trafast laser heating, International Journal of Heat and Mass Transfer, 44,
1725-1734

Symulacja numeryczna procesów cieplnych zachodzących

w wielowarstwowych mikroobszarach poddanych działaniu ultrakrótkiego

impulsu laserowego

Streszczenie

W pracy przedstawiono model matematyczny, algorytm numeryczny i przykła-
dy symulacji dotyczących przebiegu procesów cieplnych w wielowarstwowymmikro-
obszarze nagrzewanym ultraszybkim impulsem generowanym przez laser. Równanie
opisujące przebieg procesu odpowiada modelowi z dwoma opóźnieniami wynikają-



396 E. Majchrzak et al.

cymi z czasu relaksacji i czasu termalizacji. Na etapie obliczeń numerycznych wy-
korzystanometodę różnic skończonych. Algorytm bazuje na sztucznej dekompozycji
obszaru wielowarstwowego, przy czym wzajemne oddziaływania między warstwami
uwzględniono poprzez założenie ciągłości strumienia ciepła i pola temperatury na po-
wierzchniach kontaktu. Biorąc pod uwagę geometrię obszaru, rozpatrywano zadanie
jednowymiarowe.

Manuscript received September 22, 2008; accepted for print November 7, 2008