Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

56, 2, pp. 393-402, Warsaw 2018
DOI: 10.15632/jtam-pl.56.2.393

IMPLICIT SCHEME OF THE FINITE DIFFERENCE METHOD FOR THE

SECOND-ORDER DUAL PHASE LAG EQUATION

EWA MAJCHRZAK

Silesian University of Technology, Gliwice, Poland

e-mail: ewa.majchrzak@polsl.pl

Bohdan Mochnacki

University of Occupational Safety Management, Katowice, Poland

e-mail: bmochnacki@wszop.edu.pl

The second-orderdual phase lag equation(DPLE)asamathematicalmodel of themicroscale
heat transfer is considered. It is known that the starting point determining the final form of
this equation is the generalized Fourier law in which two positive constants (the relaxation
and thermalization times) appear. Depending on the order of the generalized Fourier law
expansion into theTaylor series, different formsof theDPLEcanbe obtained.Asanexample
of the problemdescribed by the second-orderDPLE equation, thermal processes proceeding
in the domain of a thin metal film subjected to a laser pulse are considered. The numerical
algorithm is based on an implicit scheme of the finite difference method. At the stage of
numerical modeling, the first, second and mixed order of the dual phase lag equation are
considered. In the final part of the paper, examples of different solutions are presented and
conclusions are formulated.

Keywords:microscale heat transfer, dual phase lagmodel, implicit schemeof finite difference
method

1. Introduction

The Fourier heat conduction model is based on the assumption of instantaneous propagation
of the thermal wave in the domain considered. Intuitively, this approach seems to be incorrect,
but it has worked for solving a number of macroscopic heat conduction problems. However, it
turned out that for certain non-typical materials with a complex internal structure, the Fourier
model is insufficient (Roetzel et al., 2003). Even more, deviations from the real course of the
process can be seen in the case of microscale heat transfer.

It is obvious that accumulating enough energy to transfer to the nearest neighborhoodwould
take time in theprocess ofheat transfer (Zhang, 2007). So, the lag timeof theheatflux in relation
to the temperature gradient referred to as “a relaxation time” was introduced by Cattaneo
(1948) and Vernotte (1958), and the appropriate energy equation (a hyperbolic PDE) became
known as the Cattaneo-Vernotte equation. In the recent years, the heat conduction model in
which two delay times appear has become more and more popular. This model is called the
dual-phase lag one (Zhang, 2007; Tzou, 2015). The starting point for considerations is the
generalized formof theFourier law, e.g. (Faghri et al., 2010; Smith andNorris, 2003). Depending
on the number of terms in the Taylor series expansion of this law, different forms of the dual
phase lag equation (DPLE) can be obtained (see Section 2). The lag times appearing in DPLE
are called the relaxation time and the thermalization time. Some simple tasks described by
this equation (supplemented by appropriate boundary and initial conditions) can be solved
analytically, e.g. (Ciesielski, 2017a;Tang andAraki, 1999;Askarizadeh et al., 2017;Mohammadi-



394 E.Majchrzak, B.Mochnacki

-Fakhar andMomeni-Masuleh, 2010). However, most of the practical problems have been solved
using numerical methods. Examples of such solutions in the field of themicroscale heat transfer
may be the papers (Majchrzak and Mochnacki, 2014; Ciesielski, 2017b; Dai and Nassar, 2000;
Mochnacki and Paruch, 2013; Chen and Beraun, 2001) concerning the first-order DPLE.
The similar problems have been considered for non-homogeneous (multilayered) domains. In

this place, the papers (Majchrzak et al., 2009; Qiu et al., 1994; Al-Nimr et al., 2004;Wang et al.,
2006, 2008) can be (as the examples) mentioned. The correct form of the boundary conditions
between subdomains (here, the macroscopic boundary conditions are often used, which is a
significant simplification) can be found in (Ho et al., 2003) while the detailed mathematical
considerations were shown in (Majchrzak and Kałuża, 2017). In turn, in the paper (Majchrzak
and Mochnacki, 2016), the problem of stability condition (explicit scheme of the FDM) was
analyzed.
The numerical solutions concerning the second-order DPLE (based on the finite difference

method) are the subject of works prepared by Castro et al. (2016) and Deng et al. (2017). The
similar problems are discussed in the paper presented, but the wider class of equations and the
other numerical algorithm are taken into account.
The applications of DPLE in the scope of bioheat transfer will not be discussed here.

2. Dual-phase lag model

The following well known thermal diffusion equation is considered

c
∂T(X,t)

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

where c is a volumetric specific heat, q is a heat flux vector, Q is a capacity of the internal
volumetric heat source,X, t denote the geometrical co-ordinates and time.
The relationship between the heat flux q and the temperature gradient ∇T is given in the

form of the generalized Fourier law (Zhang, 2007; Smith and Norris, 2003), namely

q(X,t+ τq)=−λ∇T(X,t+ τT) (2.2)

where λ is thermal conductivity, τq and τT are the relaxation time and thermalization time,
respectively. The relaxation time τq is themean time for electrons to change their energy states,
while the thermalization time τT is the mean time required for electrons and lattice to reach
equilibrium.
Using theTaylor series expansions, the following second-order approximation of formula (2.2)

can be taken into account

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

∂t
+
τ2q
2

∂2q(X,t)

∂t2
=−λ

[

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

∂t
+
τ2T
2

∂2∇T(X,t)

∂t2

]

(2.3)

which means

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

∂t
+
τ2q
2

∂2q(X,t)

∂t2
+λ∇T(X,t)+λτT

∂∇T(X,t)

∂t
+λ
τ2T
2

∂2∇T(X,t)

∂t2
(2.4)

From equation (2.4) it results that

−∇·q(X,t) = τq
∂[∇·q(X,t)]

∂t
+
τ2q
2

∂2[∇·q(X,t)]

∂t2
+∇[λ∇T(X,t)]

+ τT
∂{∇[λ∇T(X,t)]}

∂t
+
τ2T
2

∂2{∇[λ∇T(X,t)]}

∂t2

(2.5)



Implicit scheme of the finite difference method for... 395

The last dependence is introduced in to equation (2.1), and then

c
∂T(X,t)

∂t
= τq
∂[∇·q(X,t)]

∂t
+
τ2q
2

∂2[∇·q(X,t)]

∂t2
+∇[λ∇T(X,t)]

+ τT
∂{∇[λ∇T(X,t)]}

∂t
+
τ2T
2

∂2{∇[λ∇T(X,t)]}

∂t2
+Q(X,t)

(2.6)

Equation (2.1) can also be written as

∇·q(X,t) =−c
∂T(X,t)

∂t
+Q(X,t) (2.7)

Putting equation (2.7) into (2.6), one obtains

c
∂T(X,t)

∂t
= τq
∂

∂t

[

−c
∂T(X,t)

∂t
+Q(X,t)

]

+
τ2q
2

∂2

∂t2

[

−c
∂T(X,t)

∂t
+Q(X,t)

]

+∇[λ∇T(X,t)]+ τT
∂{∇[λ∇T(X,t)]}

∂t
+
τ2T
2

∂2{∇[λ∇T(X,t)]}

∂t2
+Q(X,t)

(2.8)

Assuming the constant value of the volumetric specific heat c, one has

c
[∂T(X,t)

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

∂t2
+
τ2q
2

∂3T(X,t)

∂t3

]

=∇[λ∇T(X,t)]+ τT
∂{∇[λ∇T(X,t)]}

∂t

+
τ2T
2

∂2{∇[λ∇T(X,t)]}

∂t2
+Q(X,t)+ τq

∂Q(X,t)

∂t
+
τ2q
2

∂2Q(X,t)

∂t2

(2.9)

Additionally, for λ= const the last equation takes form

c
[∂T(X,t)

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

∂t2
+
τ2q
2

∂3T(X,t)

∂t3

]

=λ∇2T(X,t)+λτT
∂[∇2T(X,t)]

∂t

+λ
τ2T
2

∂2[∇2T(X,t)]

∂t2
+Q(X,t)+ τq

∂Q(X,t)

∂t
+
τ2q
2

∂2Q(X,t)

∂t2

(2.10)

As previously mentioned, dual phase lag equation (2.10) is often simplified by omitting appro-
priate components. For example, in several works (e.g. Tzou, 1995) the second order Taylor
expression of heat flux and the first order Taylor expression of the temperature gradient are
applied to take into account the phase lagging behavior. Ignoring the inner heat source (as in
Tzou, 1995), the governing equation of temperature based on the DPLmodel is the following

c
[∂T(X,t)

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

∂t2
+
τ2q
2

∂3T(X,t)

∂t3

]

=λ∇2T(X,t)+λτT
∂[∇2T(X,t)]

∂t
(2.11)

It is also possible to consider the energy equation in the form (assuming thatQ(X,t)= 0)

c
[∂T(X,t)

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

∂t2

]

=λ∇2T(X,t)+λτT
∂[∇2T(X,t)]

∂t
+λ
τ2T
2

∂2[∇2T(X,t)]

∂t2
(2.12)

The most popular DPLE results from the assumption that the first-order approximation of
formula (2.2) is used, and then (e.g. Tang andAraki, 1999; Al-Nimr et al., 2004; Majchrzak and
Mochnacki, 2014)

c
[∂T(X,t)

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

∂t2

]

=λ∇2T(X,t)+λτT
∂[∇2T(X,t)]

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

∂Q(X,t)

∂t
(2.13)

One can see that for τT =0, DPLE (2.13) takes form of the Cattaneo-Vernotte equation, while
for τq =0 and τT =0 the well knownmacroscopic Fourier equation is obtained.



396 E.Majchrzak, B.Mochnacki

Taking into account the numerical examples presented in the final part of the paper, a
modified form of the Neumann boundary condition must still be formulated, namely

qb(X,t)+ τq
∂qb(X,t)

∂t
+
τ2q
2

∂2qb(X,t)

∂t2

=−λ
[

n ·∇T(X,t)+ τT
∂[n ·∇T(X,t)]

∂t
+
τ2T
2

∂2[n ·∇T(X,t)]

∂t2

]

(2.14)

where n ·∇T(X,t) denotes normal derivative and qb(X,t) is the known boundary heat flux. In
the case of simplified forms of theDPLE, the appropriate components in condition (2.14) should
be neglected.

3. Formulation of the problem

Thermal processes proceeding in a thin metal film subjected to laser pulse are considered. A
1D problem is analyzed (heat transfer in the direction perpendicular to the layer is taken into
account). The front surface x = 0 is irradiated by a laser pulse and according to (Tang and
Araki, 1999; Kaba andDai, 2005), the conductional heat transfer in the domain considered can
be modeled using the DPLE in which the volumetric heat source Q(x,t) is introduced. At the
same time, forx=0andx=L, thenon-fluxconditions shouldbe assumed.The laser irradiation
is described by the following source term

Q(x,t)=

√

β

π

1−R

tpδ
I0exp

[

−
x

δ
−β
(t−2tp)

2

t2p

]

(3.1)

where I0 is the laser intensity, tp is the characteristic time of the laser pulse, δ is the optical
penetration depth,R is the reflectivity of the irradiated surface, and β=4ln2.
In the most general case, the following DPLE is considered::

— for 0<x<L

∂T(x,t)

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

∂t2
+wq
τ2q
2

∂3T(x,t)

∂t3
= a
∂2T(x,t)

∂x2
+aτT

∂3T(x,t)

∂t∂x2

+wTa
τ2T
2

∂4T(x,t)

∂t2∂x2
+
1

c
Q(x,t)+

τq
c

∂Q(x,t)

∂t
+wq
τ2q
2c

∂2Q(x,t)

∂t2

(3.2)

where a=λ/c is the diffusion coefficient,wT andwq are bivalent parameters. HerewT =1 and
wq =1. For the “simplified” forms of DPLE, they are equal to (0,1), (1,0) and (0,0).
As previously mentioned, qb(0, t) = qb(L,t) = 0 and the appropriate boundary conditions

are of the form (Eq. (2.14)):
— for x=0

∂T(x,t)

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

∂t∂x
+wT

τ2T
2

∂3T(x,t)

∂t2∂x
=0 (3.3)

— for x=L

∂T(x,t)

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

∂t∂x
+wT

τ2T
2

∂3T(x,t)

∂t2∂x
=0 (3.4)

The initial condition is also given for t=0

T(x,0)=Tp
∂T(x,t)

∂t

∣

∣

∣

∣

∣

t=0

=u(x)
∂2T(x,t)

∂t2

∣

∣

∣

∣

∣

t=0

= v(x) (3.5)

where Tp is the initial temperature, while u(x) and v(x) are known functions.



Implicit scheme of the finite difference method for... 397

4. Numerical algorithm

The algorithm presented below is based on the implicit scheme of the finite difference method
(FDM).

Let T
f
i = T(xi,f∆t), where ∆t is the time step, xi = ih (h is the geometrical mesh step)

and f =0,1, . . . ,F. Taking into account initial conditions (3.5), on the assumption that u(x)=
v(x) = 0, one has T0i = T

1
i = T

2
i = Tp. For the transition t

f−1 → tf (f ­ 3), the approximate
form of equation (3.2) resulting from the introduction of adequate differential quotients is as
follows

T
f
i −T

f−1
i

∆t
+ τq
T
f
i −2T

f−1
i +T

f−2
i

(∆t)2
+wq
τ2q
2

T
f
i −3T

f−1
i +3T

f−2
i −T

f−3
i

(∆t)3

= a
T
f
i−1−2T

f
i +T

f
i+1

h2
+
aτT
∆t

(T
f
i−1−2T

f
i +T

f
i+1

h2
−
T
f−1
i−1 −2T

f−1
i +T

f−1
i+1

h2

)

(4.1)

+wT
aτ2T
2(∆t)2

(T
f
i−1−2T

f
i +T

f
i+1

h2
−2
T
f−1
i−1 −2T

f−1
i +T

f−1
i+1

h2
+
T
f−2
i−1 −2T

f−2
i +T

f−2
i+1

h2

)

+
1

c
Q
f
i +
τq
c

(∂Q

∂t

)f

i
+wq
τ2q
2c

(∂2Q

∂t2

)f

i

After mathematical transformations, one has

−
a[2(∆t)2+2τT∆t+wTτ

2
T ]

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

[2(∆t)2+2τq∆t+wqτ
2
q

2(∆t)3

+
2a[2(∆t)2+2τT∆t+wTτ

2
T ]

2h2(∆t)2

]

T
f
i −
a[2(∆t)2+2τT∆t+wTτ

2
T ]

2h2(∆t)2
T
f
i+1

=
2(∆t)2+4τq∆t+3wqτ

2
q

2(∆t)3
T
f−1
i −

2τq∆t+3wqτ
2
q

2(∆t)3
T
f−2
i

+
wqτ

2
q

2(∆t)3
T
f−3
i −

aτT(∆t+wTτT)

h2(∆t)2
(T
f−1
i−1 −2T

f−1
i +T

f−1
i+1 )

+
awTτ

2
T

2h2(∆t)2
(T
f−2
i−1 −2T

f−2
i +T

f−2
i+1 )+

1

c
Q
f
i +
τq
c

(∂Q

∂t

)f

i
+wq
τ2q
2c

(∂2Q

∂t2

)f

i

(4.2)

Denoting

A=−
a[2(∆t)2+2τT∆t+wTτ

2
T ]

2h2(∆t)2
B=
2(∆t)2+2τq∆t+wqτ

2
q

2(∆t)3
−2A

C
f
i =
2(∆t)2+4τq∆t+3wqτ

2
q

2(∆t)3
T
f−1
i −

2τq∆t+3wqτ
2
q

2(∆t)3
T
f−2
i +

wqτ
2
q

2(∆t)3
T
f−3
i

−
aτT(∆t+wTτT)

h2(∆t)2
(T
f−1
i−1 −2T

f−1
i +T

f−1
i+1 )+

awTτ
2
T

2h2(∆t)2
(T
f−2
i−1 −2T

f−2
i +T

f−2
i+1 )

+
1

c
Q
f
i +
τq
c

(∂Q

∂t

)f

i
+wq
τ2q
2c

(∂2Q

∂t2

)f

i

(4.3)

one obtains

AT
f
i−1+BT

f
i +AT

f
i+1 =C

f
i (4.4)

The FDM equation resulting from the boundary condition for x=0 is of the form

T
f
1 −T

f
0

h
+
τT
∆t

(T
f
1 −T

f
0

h
−
T
f−1
1 −T

f−1
0

h

)

+
wTτ

2
T

2(∆t)2

(T
f
1 −T

f
0

h
−2
T
f−1
1 −T

f−1
0

h
+
T
f−2
1 −T

f−2
0

h

)

=0

(4.5)



398 E.Majchrzak, B.Mochnacki

or

− [2(∆t)2+2τT∆t+wTτ
2
T ]T
f
0 +[2(∆t)

2+2τT∆t+wTτ
2
T ]T
f
1

=(2τT∆t+2wTτ
2
T)(T

f−1
1 −T

f−1
0 )−wTτ

2
T(T

f−2
1 −T

f−2
0 )

(4.6)

Let us denote

D=2(∆t)2+2τT∆t+wTτ
2
T E=2τT∆t+2wTτ

2
T (4.7)

then

−DT
f
0 +DT

f
1 =E(T

f−1
1 −T

f−1
0 )−wTτ

2
T(T

f−2
1 −T

f−2
0 ) (4.8)

In a similar way, for x=L, one has

−DT
f
n−1+DT

f
n =E(T

f−1
n −T

f−1
n−1)−wTτ

2
T(T

f−2
n −T

f−2
n−1) (4.9)

So, the final form of the system of equations corresponding to the transition tf−1 → tf (f ­ 3)
is the following

−DT
f
0 +DT

f
1 =E(T

f−1
1 −T

f−1
0 )−wTτ

2
T(T

f−2
1 −T

f−2
0 )

AT
f
i−1+BT

f
i +AT

f
i+1 =C

f
i i=1,2, . . . ,n−1

−DT
f
n−1+DT

f
n =E(T

f−1
n −T

f−1
n−1)−wTτ

2
T(T

f−2
n −T

f−2
n−1)

(4.10)

So, the transition from tf−1 to tf (f ­ 3) requires solving of the system of equations with a
three-bandmain matrix which is the fastest solved using the Thomas algorithm.

5. Examples of computations

Thin metal films (L = 100nm) made of chromium, nickel and gold have been considered.
The surface x = 0 of the domain is subjected to the laser pulse. The parameters determi-
ning the capacity of the internal heat source (Eq. (3.1)) are equal to I0 =13.7J/m

2, tp =0.1ps,
δ = 15.3nm, R = 0.93. The initial temperature of the domain equals Tp = 300K, while the
initial values of functions are u(x) = 0, v(x) = 0. Differential mesh parameters are n = 1000,
∆t=0.0001ps.
At the stage of numerical computations, constant values of thermophysical parameters have

been assumed (mainly due to lack of other data in the literature) – see Table 1.

Table 1.Thermophysical parameters (Tzou, 2015)

Chromium Gold Nickel

c [MJ/(m3K)] 3.21484 2.4897 4

λ [W/(mK)] 93 315 90.8

τq [ps] 0.136 8.5 0.82

τT [ps] 7.86 90 10

Computations have been performed in versions corresponding towT =0,wq =0 (first-order
DPLE), wT = 1, wq = 1 (second-order DPLE), wT = 0, wq = 1 and wT = 1, wq = 0 (mixed
orderDPLE).Additionally, for comparative purposes, numerical solutions of the classical Fourier
problems have been also found. The results are presented in the form of heating/cooling curves
at the irradiated surface. The set of solutions for the chromium layer is shown in Fig. 1. For
the other materials (Figs. 2 and 3), the solutions corresponding to the Fourier model, wT =0,
wq = 0 and wT = 1, wq = 1 are distinguished. The discussion of the results obtained will be
carried out in the next Section.



Implicit scheme of the finite difference method for... 399

Fig. 1. Temperature history at the irradiated surface for different models (chromium)

Fig. 2. Temperature history at the irradiated surface for different models (gold)

6. Conclusions

Different (in the sense of the order) models using the dual phase lag equation give different
results. Here, one can see some regularities. In relation to the model based on the second-order
DPLE, the solution resulting from the first-order equation is clearly overstated. This is the case
for all the materials in question. The fact that the Fourier model gives a solution over DPLE
has been repeatedly confirmed in numerous papers. This is a natural consequence of the delay



400 E.Majchrzak, B.Mochnacki

Fig. 3. Temperature history at the irradiated surface for different models (nickel)

times introduced. In the case of mixed models, the omission of the component containing τ2T
(Eq. (2.11)) leads to results close to the solution of the first-order DPLE – see Fig. 1. On the
other hand, the omission of the component containing τ2q (Eq. (2.12)) gives a solution similar to
the solution of the second-order DPLE.The same trend is observed for the remainingmaterials.
This results from themuch larger (in the case ofmetals) value of the thermalization time versus
the relaxation one. Therefore, more components of the Taylor series should be included on the
right hand side of the generalized Fourier law. Summing up, the problems connected with the
modeling of thermal processes in metal microdomains should be solved using the second-order
dual phase lag equation. If the delay times vary less, then the solution based on the first-order
model is sufficiently accurate.

Acknowledgement

The paper and research were financed within the project 2015/19/B/ST8/01101 sponsored by The

National Science Centre (Poland).

References

1. Al-Nimr M.A., Naji M., Abdallah R.I., 2004, Thermal behavior of a multi-layered thin slab
carryingperiodic signalsunder the effectof thedual-phase-lagheat conductionmodel, International
Journal of Thermophysics, 25, 3, 949-966

2. AskarizadehH., Baniasadi E., AhmadikiaH., 2017,Equilibrium and non-eqilibrium thermo-
dynamic analysis of high-order dual-phase lag heat conduction, International Journal of Heat and
Mass Transfer, 104, 301-309

3. Castro M.A., Rodriguez F., Cabrera J., Martin, J.A., 2016, A compact difference sche-
me for numerical solution of second order dual-phase-lagging models of microscale heat transfer,
Journal of Computational and Applied Mathematics, 201, 432-440

4. Cattaneo M.C., 1948, Sulla conduzione de calor, Atti de Seminario Matematico e Fisico Della
Universita di Modena, 3, 3, 3-21



Implicit scheme of the finite difference method for... 401

5. Chen J. K., Beraun J. E., 2001, Numerical study of ultrashort laser pulse interactions with
metal films,Numerical Heat Tranfer A – Applications, 40, 1-20

6. Ciesielski M., 2017a, Analytical solution of the dual phase lag equation describing the laser
heating of thin metal film, Journal of Applied Mathematics and Computational Mechanics, 16, 1,
33-40

7. Ciesielski M., 2017b, Application of the alternating direction implicit method for numerical
solution of the dual phase lag equation, Journal of Theoretical and Applied Mechanics, 55, 3,
839-852

8. Dai W., Nassar R., 2000, A compact finite difference scheme for solving a three-dimensional
heat transport equation in a thin film, Numerical Methods for Partial Differential Equations, 16,
441-458

9. Deng D., Jiang Y., Liang D. L., 2017, High-order finite difference methods for a second-order
dual-phase-laggingmodels ofmicroscaleheat transfer,AppliedMathematics andComputation,309,
31-48

10. Faghri A., ZhangY.,Howell J., 2010,AdvancedHeat andMass Transfer, GlobalDigital Press

11. Ho A.J.R., Kuo Ch.P., Jiaung W.S., 2003, Study of heat transfer in multilayered structure
within the framework of dual-phase-lag heat conduction model using lattice Boltzmann method,
International Journal of Heat and Mass Transfer, 46, 55-69

12. Kaba I., Dai W., 2005, A stable three-level finite difference scheme for solving the parabolic
two-stepmodel in a 3Dmicro-sphere heated by ultrashort-pulsed lasers, Journal of Computational
and Applied Mathematics, 181, 125-147

13. Majchrzak E., Kałuża G., 2017, Analysis of thermal processes occurring in the heated multi-
layeredmetal films using the dual-phase lag model,Archives of Mechanics, 69, 4-5, 275-287

14. Majchrzak E., Mochnacki B., 2014, Sensitivity analysis of transient temperature field in mi-
crodomains with respect to dual-phase-lag-model parameters, International Journal for Multiscale
Computational Engineering, 12, 1, 65-77

15. Majchrzak E., Mochnacki B., 2016, Dual-phase lag equation. Stability conditions of a nume-
rical algorithm based on the explicit scheme of the finite difference method, Journal of Applied
Mathematics and Computational Mechanics, 15, 3, 89-96

16. Majchrzak E., Mochnacki B., Suchy J.S., 2009, Numerical simulation of thermal processes
proceeding in a multi-layered film subjected to ultrafast laser heating, Journal of Theoretical and
Applied Mechanics, 47, 2, 383-396

17. Mochnacki B., Paruch M., 2013, Estimation of relaxation and thermalization times in micro-
scale heat transfer model, Journal of Theoretical and Applied Mechanics, 51, 4, 837-845

18. Mohammadi-Fakhar V., Momeni-Masuleh S. H., 2010, An approximate analytic solution of
the heat conduction equation at nanoscale, Physics Letters A, 374, 595-604

19. Qiu T.Q., Juhasz T., Suarez C., Bron N.E., Tien C.L., 1994, Femtoscond laser heating
of multi-layers metals. II Experiments, International Journal of Heat and Mass Transfer, 37, 17,
2799-2808

20. Smith A. N., Norris P. M., 2003, Microscale Heat Transfer, Chapter 18, [In:] Heat Transfer
Handbook, JohnWilley & Sons

21. RoetzelW., PutraN.,Das S.K., 2003,Experiment and analysis for non-Fourier conduction in
materials with non-homogeneous inner structure, International Journal of Thermal Sciences, 42,
541-552

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

23. TzouD.Y., 1995,Aunifiedfieldapproach forheat conduction frommacro- tomicro-scales,Journal
of Heat Transfer, ASME, 117, 1, 8-16



402 E.Majchrzak, B.Mochnacki

24. TzouD.Y., 2015,Macro toMicroscale Heat Transfer: The Lagging Behavior, JohnWiley&Sons,
Ltd.

25. Vernotte P., 1958, Les paradoxes de la theorie continue de l’equation de la chaleur, Comptes
Rendus de l’Academie des Sciences, 246, 3154-3155

26. WangH.,DaiW.,MelnikR., 2006,Afinite differencemethod for studying thermaldeformation
in a double-layered thin film exposed to ultrashort pulsed lasers, International Journal of Thermal
Sciences, 45, 1179-1196

27. WangH.,DaiW.,HewavitharanaL.G., 2008,Afinite differencemethod for studying thermal
deformation in a double-layered thin film with imperfect interfacial contact exposed to ultrashort
pulsed lasers, International Journal of Thermal Sciences, 47, 7-24

28. Zhang Z.M., 2007,Nano/Microscale Heat Transfer, McGraw-Hill

Manuscript received November 11, 2017; accepted for print January 9, 2018