Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation

ARTICLE IN PRESS

JID: APM

[m3Gsc;November 2, 2015;16:16]

Applied Mathematical Modelling 000 (2015) 1–13

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation ´ ∗, E. Majchrzak, L. Turchan M. Jasinski Institute of Computational Mechanics and Engineering, Silesian University of Technology, Konarskiego 18A, 44-100 Gliwice, Poland

a r t i c l e

i n f o

Article history: Received 24 April 2015 Revised 8 October 2015 Accepted 15 October 2015 Available online xxx Keywords: Bioheat transfer Generalized dual-phase lag equation Arrhenius integral Finite difference method

a b s t r a c t In the paper a numerical analysis of thermal processes proceeding in a soft tissue subjected to a laser irradiation is presented. The transient bioheat transfer in the 3D domain considered is described by the generalized dual-phase lag equation. The internal heat source resulting from the laser irradiation based on the solution of the diffusion equation is taken into account. The process of the thermal destruction of the tissue is also analyzed. At the stage of numerical realization the explicit scheme of finite difference method is used. In the final part of the paper the results of computations for different variants of tissue porosity are shown. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Modeling of physical processes proceeding in biological tissues subjected to a laser beam can be divided into three steps, namely the modeling of laser energy deposition, the temperature distribution and the damage in laser irradiated tissue [1]. To describe the light propagation in biological tissues the different mathematical models can be used [2–7]. The most fundamental approach, which is difficult from the mathematical point of view, is the so-called analytical theory which starts from the Maxwell equations. On the other hand, the transport theory concerns directly the transport of light through absorbing and scattering media and in this case the radiative transport equation should be taken into account, e.g. [8–12]. To solve this equation several modifications of the discrete ordinates method (DOM) and statistical Monte Carlo (MC) methods are widely used, e.g. [13,3,14,15]. In some cases it is possible to approximate the light transport using the diffusion equation [16–20]. As the scattering generally dominates over the absorption in soft tissues for wavelengths between 650 and 1300 nm, so in this paper the diffusion approximation is applied. The Pennes equation [21] is one of the earliest bioheat transfer equations which describes the temperature distribution in the living tissues and for the constant thermophysical parameters it has a form:

ρc

∂T = λ∇ 2 T + Q, ∂t

(1)

where c [J/(kg K)] is the specific heat of tissue, ρ [kg/m3 ] is the density, λ [W/(m K)] is the thermal conductivity, T = T(x, t) is the tissue temperature (x denotes spatial coordinates, t is the time). The capacity of internal heat sources Q = Q(x, t) [W/m3 ] can be written in the following form:

Q = wcb (Tb − T ) + Qm + Qex , ∗

(2)

Corresponding author. Tel.: +48 32 237 12 04; fax: +48 32 237 12 82. ´ E-mail address: [email protected] (M. Jasinski).

http://dx.doi.org/10.1016/j.apm.2015.10.025 S0307-904X(15)00683-6/© 2015 Elsevier Inc. All rights reserved.

´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM 2

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

where w [kg/(m3 s)] is the blood perfusion rate, cb is the specific heat of blood, Tb is the arterial blood temperature, Qm is the metabolic heat source and Qex = Qex (x, t ) is the source function associated with the external heating of tissue. The Pennes bioheat equation is based on the classical Fourier’s law of heat conduction which assumes that the thermal disturbance propagates with an infinite speed. So far, it is the most frequently used model to determine the temperature distribution in biological tissues, e.g. [18,22–29]. Heterogeneous structure of biological tissues is taken into account in the other models of bioheat transfer, for example in the hyperbolic thermal wave model (Cattaneo–Vernotte equation) or the dual-phase lag model (DPLM). In the Cattaneo–Vernotte equation [30,31]:



∂ 2T ∂T cρ + τq 2 ∂t ∂t



= λ∇ 2 T + Q + τq

∂Q , ∂t

(3)

the relaxation time τ q [s] is the delay in establishing the heat flux vector, in others words, the temperature gradient always precedes the heat flux vector. In the dual-phase lag equation [32]:



cρ

∂ 2T ∂T + τq 2 ∂t ∂t



= λ∇ 2 T +λτT

∂ 2  ∂Q ∇ T + τq , ∂t ∂t

(4)

both the heat flux and temperature gradient are delayed. The delay of temperature gradient is called the thermalization time

τ T [s].

As can be seen, for τ T = 0 the DPLM is simplified to the Cattaneo–Vernotte equation and for τ T = τ q = 0 to the Pennes equation. In the discussed models the appropriate estimation of the delay times is the significant problem. Values of these parameters are usually determined experimentally [33–35] and they are ranged from 0.01 to 32 s. In recent years, the thermal wave model and the dual-phase lag equation are often used to calculate the temperature distribution in biological tissues, e.g. [11,12,36–41]. The last step in the analysis of tissue heating process is to estimate the degree of its destruction. For this purpose, the so-called Arrhenius integral is most frequently applied, although other models are also used, for example the thermal dose parameter [42,43]. Arrhenius scheme assumes the exponential dependence between the temperature and the degree of tissue destruction, so the following integral is taken into account [26,5]:

(x) = P



tF

0



exp −



E dt, Rg T (x, t )

(5)

where P [1/s] is the pre-exponential factor, E [J/mol] is the activation energy, Rg [J/(mol K)] is the universal gas constant, T (x, t) [K] is the tissue temperature, t denotes time and [0, tF ] is the time interval under consideration. A value of damage integral (x) = 1 corresponds to a 63% probability of cell death at a specific point x, while (x) = 4.6 corresponds to 99% probability of cell death at this point [44]. Besides the basic information on whether the tissue was destructed or not (in accordance with the necrosis criterion, this means (x) ≥ 1), one can deduce the time after which the tissue damage is to occur. The tissue injury integral basically refers only to the irreversible tissue damage, however, there are models which allow to take into account the withdrawal of tissue injury in the case of temporary, small local increasing of temperature [27]. Through the approach presented in [27] it is possible to determine the time after which the thermal lesion is formed. The purpose of this paper is to analyze the phenomena occurring in the laser-treated soft tissue, in which the heat conduction in the tissue is described by the generalized dual-phase lag model [45]. In this model, based on the theory of porous media the phase lag times are expressed in terms of the properties of blood and tissue, interphase convective heat transfer coefficient and the blood perfusion rate. Using this model, for different levels of tissue perfusion the degree of their destruction can be compared. Thus, the paper consists of the following parts. At first, the generalized dual-phase lag equation (GDPLE) is presented and the problem analyzed is formulated. Next, the way of GDPLE solution basing on the explicit scheme of the finite difference method is described and the examples of computations are shown. In the final part of the paper the conclusions are formulated. 2. Generalized dual-phase lag equation Heterogeneous structure of biological tissue can be described in several ways. One of them is to divide the tissue into two regions: the vascular region (blood vessels) and the extravascular region (tissue) [46,47]. Then, the porosity ε is defined as the ratio of blood volume to the total volume and two equations describing the temperature distribution in the tissue Tt and blood vessels Tb are considered, namely

(1 − ε)ρt ct and

 ερb cb

∂ Tt = (1 − ε)λt ∇ 2 Tt + α A(Tb − Tt ) + wcb (Tb − Tt ) + (1 − ε)Qmt + (1 − ε)Qex ∂t

 ∂ Tb + u · ∇ Tb = ελb ∇ 2 Tb + α A(Tt − Tb ) + w cb (Tt − Tb ) + ε Qmb + ε Qex , ∂t

(6)

(7)

´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

3

where α [W/(m2 K)] is the heat transfer coefficient, u [m/s] is the blood velocity, A [m2 /m3 ] is the volumetric transfer area between tissue and blood, c is the specific heat, ρ is the density, λ is the thermal conductivity, t is the time, w is the blood perfusion rate, Qm is the metabolic heat source and Qex is the capacity of internal heat sources associated with the external heating of tissue. The subscripts t and b represent the tissue and blood, respectively. According to the Minkowycz hypothesis [48], before reaching thermal equilibrium, the blood temperature undergoes a transient process described by the relation

ερb cb

∂ Tb = G(Tt − Tb ), ∂t

(8)

where

G = Aα +w cb

(9)

is the parameter called the coupling factor. From Eq. (8) the tissue temperature can be determined

Tt = Tb +

ε ρb cb ∂ Tb , G ∂t

(10)

and then

∂ Tt ∂ Tb ε ρb cb ∂ 2 Tb = + , ∂t ∂t G ∂t2 while

∇ 2 Tt = ∇ 2 Tb +

(11)

ε ρb cb ∂  2  ∇ Tb . G ∂t

(12)

Adding Eqs. (6) and (7), the following equation is obtained

(1 − ε)ρt ct

∂ Tt ∂T +ερb cb b + ερb cb u · ∇ Tb = (1 − ε)λt ∇ 2 Tt + ελb ∇ 2 Tb + (1 − ε)Qmt + ε Qmb + Qex . ∂t ∂t

(13)

Introducing the formulas (10)–(12) into Eq. (13) one has

Ce

∂ Tb ε(1 − ε)ρt ct ρb cb ∂ 2 Tb ε(1 − ε)λt ρb cb ∂ + ερb cb u · ∇ Tb = e ∇ 2 Tb + (∇ 2 Tb ) + ε Qmb + (1 − ε)Qmt + Qex , + ∂t G G ∂t ∂t2 (14)

where

e = ε λb + (1 − ε)λt

(15)

Ce = ε ρb cb + (1 − ε)ρt ct

(16)

and

are the effective thermal conductivity and effective heat capacity, respectively. The phase lags for heat flux and temperature gradient are defined as [47,48]

τq =

ε(1 − ε)ρt ct ρb cb

(17)

GCe

and

τT =

ε(1 − ε)λt ρb cb . G e

(18)

Then Eq. (14) can be written in the form



Ce

∂ 2T ∂ Tb +τq 2b ∂t ∂t



+ ερb cb u · ∇ Tb = e ∇ 2 Tb + e τT

∂ 2  ∇ Tb + ε Qmb + (1 − ε)Qmt + Qex . ∂t

(19)

In the generalized dual-phase lag Eq. (19) the blood temperature is unknown. Note that it is possible, with the additional assumptions, to obtain the equation in which the only unknown is the tissue temperature [45]



Ce

∂ 2 Tt ∂ Tt +τq 2 ∂t ∂t



= e ∇ 2 Tt + e τT

∂ 2  ∇ Tt + G(Tb − Tt ) ∂t

  τqCe ∂ Qmt ∂ Qex ∂ Qmb + (1 − ε) + ε + ε Qmb + (1 − ε)Qmt + Qex + , ∂t ∂t ∂t (1 − ε)ρt ct

(20)

while the blood temperature is determined from the formula (10). To summarize, the two-temperature model (20), supplemented by appropriate boundary and initial conditions allows one to determine the tissue and blood temperatures. In this model of bioheat transfer the phase lag times are expressed in terms of the properties of blood and tissue and the coupling factor between blood and tissue. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM 4

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

Fig. 1. Domain considered.

3. Formulation of the problem The domain of biological tissue (cube of dimensions l×l×l) subjected to the laser action, as shown in Fig. 1, is considered. Tissue and blood temperatures are described by the generalized dual-phase lag Eq. (20) and the ordinary differential Eq. (10), respectively. In Eq. (20) the source function Qex connected with the laser heating is defined as follows [49]

Qex = Qex (x1 , x2 , x3 , t ) = μa φ(x1 , x2 , x3 ) p(t ),

(21)

where μa [1/m] is the absorption coefficient, φ (x1 , x2 , x3 ) [W/m2 ] is the total light fluence rate and p(t) is the function equal to 1 when the laser is on and equal to 0 when the laser is off. The total light fluence rate φ is the sum of collimated part φ c and diffuse part φ d [49]

φ(x1 , x2 , x3 ) = φc (x1 , x2 , x3 ) + φd (x1 , x2 , x3 ).

(22)

The collimated fluence rate is given as [25]

φc (x1 , x2 , x3 ) = I0 exp ( − μt x1 ) exp −s(x22 + x23 ) ,

(23)

where I0 [W/m2 ] is the surface irradiance of laser, s = 2/r2 , where r is the radius of laser beam and μ t [1/m] is the attenuation coefficient defined as [2,13,5]

μt = μa + μs ,

(24)

where μa [1/m] is the absorption coefficient and μs = (1 − g)μs [1/m] is the effective scattering coefficient (μs is the scattering coefficient, g is the anisotropy factor). To determine the diffuse fluence rate the steady-state optical diffusion equation [17,18,7] should be solved:

(x1 , x2 , x3 ) ∈ :

D∇ 2 φd (x1 , x2 , x3 ) − μa φd (x1 , x2 , x3 ) + μ s φc (x1 , x2 , x3 ) = 0,

(25)

where

D=

1 3[μa + (1 − g)μs ]

(26)

is the diffusion coefficient. Eq. (25) is supplemented by boundary condition

( x1 , x2 , x3 ) ∈ :

−D n · ∇φd (x1 , x2 , x3 ) =

φd ( x 1 , x 2 , x 3 ) 2

,

(27)

where n is the outward unit normal vector. As can been seen, the diffusion approximation of light propagation in soft tissue is used in this paper. The generalized dual-phase lag Eq. (20) is supplemented by boundary condition

(x1 , x2 , x3 ) ∈ :

qt (x1 , x2 , x3 , t ) = − e n · ∇ Tt (x1 , x2 , x3 , t ) = 0.

(28)

The following initial conditions should be assumed

t=0:

Tt (x1 , x2 , x3 , t ) = Tp,

t=0:

Tb (x1 , x2 , x3 , t ) = Tp ,

 ∂ Tt (x1 , x2 , x3 , t )   =0 ∂t t=0

(29)

and (30)

where Tp is known temperature of the blood and tissue. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

5

Fig. 2. Seven-point stencil.

4. Method of solution The problem formulated is solved using the explicit scheme of finite difference method [47,50]. The uniform spatial grid of dimensions n × n × n (h is the grid spacing) is introduced and t denotes the time step. The differential quotients approximating the derivatives appearing in Eqs. (20) and (10) are presented in Appendix A. For the internal node (i, j, k) (Fig. 2) and transition t f −1 → t f , f ≥ 2, under the assumption that the metabolic heat sources Qmt , Qmb are constant, one obtains

Ce ( t + τq )

( t )

2

and

Tbf i, j,k =

Ti,fj,k =

Ce ( t + 2τq ) − G( t )2

C e τq f −2 −1 Ti,fj,k − T ( t )2 ( t )2 i, j,k e ( t + τT ) 2 f −1 e τT 2 f −2 + ∇ Ti, j,k − ∇ Ti, j,k + GTbf i, j,k + ε Qmb + (1 − ε)Qmt + Qexf −1i, j,k

t

t



(31)

G t ερb cb f −1 T T f −1 + , G t + ερb cb i, j,k G t b i, j,k

(32)

where ∇ 2 Ti,s j,k (s = f − 1 or s = f − 2) is defined in Appendix A. After some mathematical operations the tissue temperature at each internal node is calculated from the formula (i = 1, 2, …, n, j = 1, 2, …, n, k = 1, 2, …, n)

Ti,fj,k =

Ce h2 ( t + 2τq ) − Gh2 ( t )2 − 6 e t ( t + τT ) f −1 Ti, j,k Ce h2 ( t + τq )

 e t ( t + τT )  f −1 f −1 −1 −1 −1 −1 T + Ti+1, + Ti,fj−1,k + Ti,fj+1,k + Ti,fj,k−1 + Ti,fj,k+1 j,k Ce h2 ( t + τq ) i−1, j,k  e t τT  f −2 f −2 −2 −2 −2 −2 − Ti−1, j,k + Ti+1, + Ti,fj−1,k + Ti,fj+1,k + Ti,fj,k−1 + Ti,fj,k+1 j,k 2 Ce h ( t + τq )

Ce h2 τq − 6 e t τT f −2 ( t )2 f −1 − T + GTb i, j,k + ε Qmb + (1 − ε)Qmt + Qexf −1 . i, j,k i, j,k Ce ( t + τq ) Ce h2 ( t + τq ) +

(33)

Because the explicit scheme of finite difference method is considered, the stability criterion should be determined. Using the von Neumann stability analysis one obtains

2Ce h2 ( t + 2τq ) − Gh2 ( t )2 − 12 e t ( t + 2τT ) > 0

(34)

Ce h2 ( t + 2τq ) − 12 e t τT > 0.

(35)

and

To take into account the boundary conditions (28), the so-called fictitious nodes located outside the domain considered can be introduced. These nodes are located at the distance of h/2 from the boundary of domain (in the directions perpendicular to the axes x1 , x2 and x3 ). It is assumed that the temperatures at these nodes are equal to (cf. Figs. 1 and 2): T(−h/2, x2 , x3 , t f−1 ) = T(h/2, x2 , x3 , t f−1 ), T(l + h/2, x2 , x3 , t f−1 ) = T(l − h/2, x2 , x3 , t f−1 ), T(x1 ,−l/2 − h/2, x3 , t f−1 ) = T(x1 ,−l/2 + h/2, x3 , t f−1 ), T(x1 , l/2 + h/2, x3 , t f−1 ) = T(x1 , l/2 − h/2, x3 , t f−1 ), T(x1 , x2 ,−l/2 − h /2, t f−1 ) = T(x1 , x2 ,−l/2 + h/2, t f−1 ), T(x1 , x2 , l/2 + h /2, t f−1 ) = T(x1 , x2 , l/2 − h/2, t f−1 ). From the initial conditions (29) and (30) results that the tissue and blood temperatures for f = 0 and f = 1 are equal to the initial temperature: Ti,0j,k = Ti,1j,k = Tp . At the beginning, in order to determine the function Qex at the internal nodes (i, j, k) (cf. Eq. (33)) the steady-state optical diffusion Eq. (25) is solved. Differential grid with spacing h/2 is here used and the nodes are located at the vertexes of the cubes (cf. Fig. 3). Using such a differential grid, it is easier to take into account the boundary condition (27) because part of the nodes is located exactly on the boundary . In addition, one can distinguish common grid nodes for temperature and diffuse fluence rate. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM 6

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

Fig. 3. Differential grid – temperature nodes and fluence rate nodes (cross-section of the domain).

Fig. 4. The form of the function (21) for x1 = 0.

Next, in Eq. (25) the appropriate differential quotients are used (cf. Appendix A) and then for the internal node (i, j, k) (i = 1, 2, …, 2n − 1, j = 1, 2, …, 2n − 1, k = 1, 2, …, 2n − 1) one obtains

φd i, j,k =

  4D μ s h2 φc i, j,k . φd i−1, j,k + φd i+1, j,k + φd i, j−1,k + φd i, j+1,k + φd i, j,k−1 + φd i, j,k+1 + 2 24D + μa h 24D + μa h2

(36)

For the boundary nodes located on the wall x1 = 0 one has (cf. Eq. (27))

x1 = 0 :

D

φd 1,i, j − φd 0,i, j h/2

=

φd 0,i, j 2

→

φd 0,i, j =

4D φ . 4D + h d 1,i, j

(37)

In a similar manner the approximation of the boundary condition (27) for the remaining faces of the cube is formulated. The system of Eq. (36) is solved here using the iterative method. Next, for each transition t f −1 → t f , f ≥ 2 at first the blood temperature using Eq. (32) is determined, then on the basis of Eq. (33) the tissue temperature is calculated and finally the value of the Arrhenius integral (cf. formula (5)) is estimated



f i, j,k

=P

f 



exp −

l=1



E Rg Ti,l j,k

t.

(38)

5. Results of computations The uniform differential mesh with 200 × 200 × 200 nodes is introduced (edge of the cube: l = 0.02 m), so the grid spacing is equal to h = 0.0001 m. The proper estimation of source function (21) connected with the laser pulse requires the small spacing of differential grid, because this function is dependent on the spatial variables. The testing calculations confirmed that the sufficient approximation is obtained under the assumption that h ≤ 0.0001 m (cf. Appendix B). In Fig. 4 the distribution of function Qex (0, x2 , x3 ) is shown (radius of laser beam is equal to r = 0.001 m, absorption and scattering coefficients are equal to μa = 50 1/m and μs = 16900 1/m, respectively, surface irradiance of laser is equal to I0 = 5·105 W/m2 and the anisotropy factor: g = 0.952). In the Table 1 the values of parameters assumed at the stage of numerical computations are collected [51,50]. For three different variants of porosity (Table 2) the calculations are performed. The time step is equal to t = 0.01 s and duration of laser pulse equals 70 s. It was assumed that in the Arrhenius integral (38) the activation energy is equal to E = 6.67·105 J/mol, pre-exponential factor is equal to P = 1.98·10106 1/s and the universal gas constant: Rg = 8.314472 J/(mol K) [52]. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

7

Table 1 Values of parameters [51,50]. Parameter

Tissue

Blood

Thermal conductivity [W/(m K)] Specific heat [J/(kg K)] Density [kg/m3 ] Metabolic heat source [W/m3 ] Initial temperature [°C]

λt = 0.5

λb = 0.5 cb = 3770 ρ b = 1060 Qmb = 250 Tp = 37

ct = 4000 ρ t = 1000 Qmt = 250 Tp = 37

Table 2 Variants of porosity [45]. Variant

ε

G [W/(m3 K)]

τ q [s]

τ T [s]

w [kg/(m3 s)]

1 2 3

0.0041 0.0357 0.1637

34 785.174 79 102.601 96 479.910

0.46772 1.74116 5.67173

0.46771 1.74110 5.67085

1 3 5

Table 3 Maximum temperatures.

ε = 0.0041

ε = 0.0357

ε = 0.1637

Tt [°C]

Tb [°C]

Tt [°C]

Tb [°C]

Tt [°C]

Tb [°C]

67.3633

67.3233

67.0388

66.8799

65.8366

65.1198

Fig. 5. Tissue temperature history at the node N1 .

In Fig. 5 the tissue temperature history at the node N1 (0, 0, 0) obtained for three variants of porosity is presented. As can be seen, with the increase in porosity, during the heating stage the tissue temperature is lower, while during the cooling stage the temperature is higher. In Table 3 the maximum temperature obtained for each variants of calculations is shown. The GDPLE with the first variant of porosity predicts the highest temperature, while in the case of the third variant of porosity the temperature is the lowest. In Fig. 6 the history of blood and tissue temperatures at the point N1 for second and third variants of porosity is presented. It is visible, that when the porosity is lower the differences between the temperatures Tt and Tb are smaller. It is obvious that after dividing the domain shown in Fig. 1 by planes: x1 − x2 (x3 = 0) and x1 − x3 (x2 = 0), four symmetrical parts are distinguished. In Fig. 7 the tissue temperature distribution for ε = 0.0041 after the time 70 s in ¼ of the domain is presented. Fig. 8 illustrates the temperature distribution in the section x1 − x3 (x2 = 0) for time 30 s. The substantial differences in temperature distributions depending on the porosity of the tissue are observed. In addition, for larger values of porosity the smaller tissue volume is heated by laser. As was mentioned before, the destruction model based on the Arrhenius equation takes into account not only the temperature but also the time of tissue heating. In Fig. 9 the Arrhenius integral history at the node N1 is presented. It is visible that as the first the tissue destruction occurs in the case of the lowest value of porosity. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

JID: APM 8

ARTICLE IN PRESS

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

Fig. 6. Tissue and blood temperature history at the node N1 .

Fig. 7. Temperature distribution, ¼ part of the domain, porosity ε = 0.0041, time 70 s. Table 4 Tissue destruction (Vd – volume of destructed tissue).

ε = 0.0041

ε = 0.0357

ε = 0.1637

Vd [mm3 ]

%V

Vd [mm3 ]

%V

Vd [mm3 ]

%V

17.563

0.220

16.033

0.200

13.161

0.165

The results which are presented in Fig. 10 refer to the volume of destructed tissue. In Fig. 10 the Arrhenius integral spatial distribution after 200 s (end of analysis, Tt < 39 °C) is shown, while in the Table 4 the volume of destructed tissue Vd and the percent destruction of the total volume of the domain considered are collected. As might have been expected, for the lowest porosity the volume of destructed tissue is maximal. As is apparent from the results of the calculation, the change of the porosity, the basic parameter of GDPLE equation, leads to a diversification of the heating and cooling process of biological tissue. With increasing in porosity both heating and cooling of tissue proceeds more slowly (Fig. 6). The same refers to the blood temperature (Fig. 6). Another visible effect of tissue heating process described by the GDPLE equation is greater interval between the temperature of the tissue and blood (Fig. 6) – for the 3rd variant of the calculations (ε = 0.1637) and time equal to 11.7 s the differences reach 6 °C. In the final phase of the process under consideration, when the tissue reaches temperature close to the initial one, the differences between the temperatures ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

JID: APM

ARTICLE IN PRESS

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

9

Fig. 8. Temperature distribution, intersection x1 − x3 (x2 = 0), time 30 s.

Fig. 9. Arrhenius integral history at the node N1 .

obtained for different values of porosity become smaller. Relatively small differences in temperature also apply to the maximum temperatures reached by tissue and blood (Table 3 and Fig. 6). In order to estimate the porosity ε it is assumed that the vessels of diameters db are uniformly distributed in the tissue and the entire cross-section area can be considered as a set of the recurrent hexagons [45]. The diameters of the circles which are equivalent to hexagons are equal to d and the porosity can be determined as ε = db2 /d2 . On the other hand, the heat transfer coefficient α (cf. Eq. (9)) is estimated from the Nusselt number α = Nuλb /db (Nu = 4.93) which means that α is dependent on the vessels diameter db . The volumetric transfer area between tissue and blood is equal to A = 4ε /db [45], which means that A is dependent on db and d. In addition, the blood perfusion rate w is connected, among others, with the type of tissue (liver, kidney, muscle etc.). Summing up, in each variant of the computations not only the porosity ε (Table 2) but also the vessel diameter db (0.00144 m, 0.00228 m, 0.00456 m) and the blood perfusion rate w (Table 2) are different. These three parameters have an influence on the values of the coupling factor G and the relaxation times τ q , τ T . Thus, in order to analyze the differences between the tissue and blood temperatures as well as the volume of destructed tissue not only the value of porosity ε but also the values of vessels diameter db and the blood perfusion rate w should be taken into account. In the article [53] the solution of generalized dual-phase lag equation concerning the highly absorbed tissues and the strongly scattering tissues heated by laser is presented. The 1D problem is considered and the blood temperature is assumed as the constant value and it is equal to the initial temperature of tissue. In the case of scattering tissues the boundary surfaces are ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM 10

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

Fig. 10. Arrhenius integral distribution, intersection x1 − x3 (x2 = 0), time 200 s.

Fig. 11. Comparison of 1D [53] and 3D solutions.

thermally insulated and the laser heating is considered as a body heat source





Qex (x1 , t ) = μa I0 C1 exp −k1

x1

δ





− C2 exp −k2

x1

δ



p(t ),

(39)

where δ is the effective penetration depth and C1 , C2 , k1 , k2 are determined by Monte Carlo solutions [53,54]. To partially compare the results presented in the paper [53] with the results obtained using the computer code prepared by the authors the following problem has been solved. The cube of dimensions 0.05 × 0.05 × 0.05 m with thermally insulated surfaces and source function described by formula (39) has been considered. The differential mesh consisted of 125 × 125 × 125 nodes, the time step was equal to t = 0.01 s. The calculations have been done under the assumption that the porosity is equal to ε = 0.0079, the coupling factor: G = 67435 W/(m3 K), blood perfusion rate: w = 1.9822 kg/(m3 s), relaxation time: τ q = 0.4756 s, thermalization time: τ T = 0.4755 s, incident laser irradiance: I0 = 3·105 W/m2 and duration of laser pulse: 5 s. The values of the remaining parameters have been taken from [53,54]. The comparison of temperature history on irradiated surface obtained using 1D model with the temperature history at the node located nearest to the origin of coordinate system obtained using 3D model (Fig. 1) is shown in Fig. 11. It is visible, that the results are very similar. A proper estimation of the porosity of the tissue and the blood perfusion rate is very important because these values depend on many factors. As has been already shown, small changes of these parameters can significantly affect the temperature distribution. Many researchers are working on the methods of temperature measurement in living tissues, e.g. [55–58] and it can be assumed that in a near future the sufficiently accurate experimental results allow for the validation of the model.

´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

ARTICLE IN PRESS

JID: APM

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

11

6. Conclusions The heat transfer process in the tissue subjected to laser irradiation is analyzed. The problem is described by generalized dualphase lag model for the tissue sub-domain and ordinary differential equation for the blood sub-domain. The explicit scheme of finite difference method is presented and the results of computations for different values of porosity are shown. It turns out that the temperature distribution and also the Arrhenius integral are essentially different for different porosities. For a lower porosity the temperatures are higher and the degree of tissue damage is greater in comparison with a higher porosity. In the presented analysis all parameters are assumed as constant values, but, at least some of them should be regarded as dependent on the temperature or the damage of tissue. For example, as known, one of the most common markers of tissue necrosis is the loss of perfusion. In the generalized dual-phase lag model the blood perfusion rate w is introduced in an implicit way via the coupling factor G (cf. Eq. (9)). It can be assumed that for the thermally-damaged tissue, when the necrosis criterion is fulfilled, the perfusion rate w = 0, which leads to a reduction in the value of the coupling factor G. The tissue damage has also an influence on the porosity value − with the increase of the tissue damage, the value of porosity should decrease, which will affect on the thermalization time τ T , the relaxation time τ q , the effective thermal conductivity e and the effective heat capacity Ce (cf. Eqs. (15)–(18)). In literature one can find other models describing the thermal processes occurring in tissues heated by laser that take into account the presence of blood vessels, e.g. [59,60]. On the other hand, there are works (e.g. [61]) in which the parameters occurring in the mathematical model are treated as interval numbers and such an approach can be applied in the future to the problems discussed here.

Appendix A The differential quotients approximating the derivatives appearing in Eqs. (20) and (10) are the following

    

∂T ∂t

f =

−1 Ti,fj,k − Ti,fj,k

i, j,k

∂ 2T ∂t2 ∂ 2T ∂ x21 ∂ 2T ∂ x22 ∂ 2T ∂ x23

f

=



i, j,k s

=



i, j,k s

=

t −1 −2 Ti,fj,k − 2Ti,fj,k + Ti,fj,k

s s Ti+1, − 2Ti,s j,k + Ti−1, j,k j,k

h2 Ti,s j+1,k − 2Ti,s j,k + Ti,s j−1,k h2

i, j,k

s

=

(A.1)

( t )2

Ti,s j,k+1 − 2Ti,s j,k + Ti,s j,k−1 h2

i, j,k

.

(A.1)

In the above dependencies, for simplicity, the indexes t or b are omitted. The expression∇ 2 Ti,s j,k is defined as

∇ 2 Ti,s j,k =

s s Ti−1, − 2Ti,s j,k + Ti+1, j,k j,k

h2

+

Ti,s j−1,k − 2Ti,s j,k + Ti,s j+1,k h2

+

Ti,s j,k−1 − 2Ti,s j,k + Ti,s j,k+1 h2

.

(A.2)

The approximate form of Eq. (25) is as follows

 D

φd i+1, j,k − 2φd i, j,k + φd i−1, j,k

(h/2)2 − μa φd i, j,k + μ s φc i, j,k = 0.

+

φd i, j+1,k − 2φd i, j,k + φd i, j−1,k (h/2)2

+

φd i, j,k+1 − 2φd i, j,k + φd i, j,k−1



(h/2)2 (A.3)

Appendix B For the 1st variant of porosity (ε = 0.0041) the calculations were also performed under the assumption that the grid spacing is equal to h1 = 0.00005 (400 × 400 × 400 nodes). In this case, in order to fulfill the stability criteria, the time step t1 = 0.0025 s has been used. The differences between the tissue temperatures max(|Tt (h1 , t1 ) − Tt (h, t1 )|) ≤ 0.16 °C are small – cf. Fig. B.1. The calculations were also done using differential grid with 200 × 200 × 200 nodes and time steps t = 0.01 s and

t2 = 0.005 s, respectively. The tissue temperatures for both time steps are almost the same: max(|Tt (h, t) − Tt (h,

t2 )|) ≤ 1·10 −8 °C. ´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

JID: APM 12

ARTICLE IN PRESS

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

Fig. B.1. Tissue temperature history at the node N1 for different grid steps.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

A.J. Welch, The thermal response of laser irradiated tissue, IEEE J. Quantum Electron. 20 (12) (1984) 1471–1481. J. Ashley, M. Welch, Optical-Thermal Response of Laser-Irradiated Tissue, Plenum Press, New York, 1995. L.A. Dombrovsky, D. Baillis, Thermal Radiation in Disperse Systems: An Engineering Approach, Begell House, New York, 2010. S.L. Jacques, B.W. Pogue, Tutorial on diffuse light transport, J. Biomed. Opt. 13 (4) (2008) 041302-1–041302-19. M.H. Niemz, Laser – Tissue Interaction: Fundamentals and Applications, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, SpringerVerlag, Berlin, Heidelberg, New York, 2007. V.V. Tuchin, Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis, vol. PM166, second ed., SPIE Press, Bellingham, WA, 2007. A.J. Welch, M.J.C. van Gemert (Eds.), Optical-Thermal Response of Laser Irradiated Tissue, Springer, 2011. L.A. Dombrovsky, J.H. Randrianalisoa, W. Lipinski, V. Timchenko, Simplified approaches to radiative transfer simulations in laser induced hyperthermia of superficial tumors, Comput. Therm. Sci. 5 (6) (2013) 521–530. L.A. Dombrovsky, V. Timchenko, M. Jackson, Indirect heating strategy of laser induced hyperthermia: an advanced thermal model, Int. J. Heat Mass Transf. 55 (17–18) (2012) 4688–4700. L.A. Dombrovsky, V. Timchenko, M. Jackson, G.H. Yeoh, A combined transient thermal model for laser hyperthermia of tumors with embedded gold nanoshells, Int. J. Heat Mass Transf. 54 (2011) 5459–5469. M. Jaunich, S. Raje, K. Kim, K. Mitra, Z. Guo, Bio-heat transfer analysis during short pulse laser irradiation of tissues, Int. J. Heat Mass Transf. 51 (2008) 5511–5521. K. Kim, Z. Guo, Multi-time-scale heat transfer modeling of turbid tissues exposed to short-pulsed irradiations, Comput. Methods Prog. Biomed. 86 (2007) 112–123. S. Banerjee, S.K. Sharma, Use of Monte Carlo simulations for propagation of light in biomedical tissues, Appl. Opt. 49 (22) (2010) 4152–4159. J.R. Howell, R. Siegel, M.P. Mengüç, Thermal Radiation Heat Transfer, CRC Press, New York, 2010. M.F. Modest, Radiative Heat Transfer, second ed., Academic Press, New York, 2003. L.A. Dombrovsky, The use of transport approximation and diffusion-based models in radiative transfer calculations, Comput. Therm. Sci. 4 (4) (2012) 297– 315. T.J. Farrel, M.S. Patterson, B. Wilson, A diffusion theory model of spatially resolved, steady state diffuse reflectance for the non-invasive determination of tissue optical properties in vivo, Med. Phys. 19 (4) (1992) 879–888. A. Fasano, D. Hömberg, D. Naumov, On a mathematical model for laser-induced thermotherapy, Appl. Math. Model. 34 (12) (2010) 3831–3840. Z. Guo, S.K. Wan, K. Kim, Ch. Kosaraju, Comparing diffusion approximation with radiation transfer analysis for light transfer in tissues, Opt. Rev. 10 (5) (2003) 415–421. J. Mobley, T. Vo-Dinh, Optical properties of tissue, in: T. Vo-Dinh (Ed.), Biomedical Photonics Handbook, CRC Press, Boca Raton (FL), 2003 2-1–2-75. H.H. Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm, J. Appl. Physiol. l (1948) 93–122. M. Ciesielski, B. Mochnacki, Application of the control volume method using the Voronoi polygons for numerical modeling of bio-heat transfer processes, J. Theor. Appl. Mech.-Pol. 52 (4) (2014) 927–935. ´ M. Dziewonski, B. Mochnacki, R. Szopa, Sensitivity of biological tissue freezing process on the changes of cryoprobe cooling rate, Mechanika Kaunas University of Technology, 2011, pp. 82–87. K. Erhart, E. Divo, A. Kassab, An evolutionary-based inverse approach for the identification of non-linear heat generation rates in living tissues using localized meshless method, Int. J. Numer. Method. Heat Fluid Flow 18 (2008) 401–414. F. Fanjul-Vélez, O.G. Romanov, J.L. Arce-Diego, Efficient 3D numerical approach for temperature prediction in laser irradiated biological tissues, Comput. Biol. Med. 39 (2009) 810–817. ´ M. Jasinski, Investigation of tissue thermal damage process with application of direct sensitivity method, Mol. Cell. Biomech. 10 (3) (2013) 183–199. ´ M. Jasinski, Modelling of tissue thermal injury formation process with application of direct sensitivity method, J. Theor. App. Mech.-Pol. 52 (4) (2014) 947–957. ´ ´ E. Majchrzak, B. Mochnacki, M. Dziewonski, M. Jasinski, Numerical modelling of hyperthermia and hypothermia processes, Adv. Mater. Res. 268–270 (2011) 257–262. B. Mochnacki, E. Majchrzak, M. Duda, Numerical modeling of thermal processes in the living tissue domain secured with a layer of protective clothing, J. Appl. Math.Comput. Mech. 13 (1) (2014) 97–102. C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, C.R. Acad. Sci. I-Math. 247 (1958) 431–433. P. Vernotte, Les paradoxes de la theorie continue de l’equation de la chaleur, C.R. Acad. Sci. I-Math. 246 (1958) 3154–3155. F. Xu, K.A. Seffen, T.J. Lu, Non-Fourier analysis of skin biothermomechanics, Int. J. Heat Mass Transf. 51 (2008) 2237–2259. P.J. Antaki, New interpretation of non-Fourier heat conduction in processed meat, J. Heat Transf.-Trans. ASME 127 (2005) 189–193. W. Kaminski, Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, J. Heat Trans.-Trans. ASME 112 (1990) 555–560.

´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025

JID: APM

ARTICLE IN PRESS

[m3Gsc;November 2, 2015;16:16]

M. Jasi´ nski et al. / Applied Mathematical Modelling 000 (2015) 1–13

13

[35] K. Mitra, S. Kumar, A. Vedavarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, J. Heat Transf.-Trans. ASME 117 (1995) 568–573. [36] E. Majchrzak, Numerical solution of dual phase lag model of bioheat transfer using the general boundary element method, CMES-Comput. Model. Eng. Sci. 69 (1) (2010) 43–60. [37] E. Majchrzak, Application of different variants of the BEM in numerical modeling of bioheat transfer processes, Mol. Cell. Biomech. 10 (3) (2013) 201–232. [38] E. Majchrzak, Ł. Turchan, The general boundary element method for 3D dual-phase lag model of bioheat transfer, Eng. Anal. Bound. Elem. 50 (2015) 76–82. [39] A. Narasimhan, S. Sadasivam, Non-Fourier bio heat transfer modelling of thermal damage during retinal laser irradiation, Int. J. Heat Mass Transf. 60 (2013) 591–597. [40] J. Zhou, J.K. Chen, Y. Zhang, Dual-phase lag effects on thermal damage to biological tissues caused by laser irradiation, Comput. Biol. Med. 39 (2009) 286–293. [41] J. Zhou, Y. Zhang, J.K. Chen, An axisymmetric dual-phase-lag bioheat model for laser heating of living tissues, Int. J. Therm. Sci. 48 (2009) 1477–1485. [42] E. Majchrzak, Ł. Turchan, Numerical estimation of thermal dose during hyperthermia treatment using the BEM, Mechanika Kaunas University of Technology, 2011, pp. 215–219. [43] S.A. Sapareto, W.C. Dewey, Thermal dose determination in cancer therapy, Int. J. Radiat. Oncol. 10 (6) (1984) 787–800. [44] I.A. Chang, U.D. Nguyen, Thermal modeling of lesion growth with radiofrequency ablation devices, Biomed. Eng. Online 3 (27) (2004) 1–9. [45] Y. Zhang, Generalized dual-phase lag bioheat equations based on nonequilibrium heat transfer in living biological tissues, Int. J. Heat Mass Transf. 52 (2009) 4829–4834. [46] A.R.A. Khaled, K. Vafai, The role of porous media in modelling of flow and heat transfer in biological tissues, Int. J. Heat Mass Transf. 46 (2003) 4989–5003. [47] E. Majchrzak, Ł. Turchan, Numerical analysis of tissue heating using the generalized dual phase lag model, in: T. Łodygowski, J. Rakowski, P. Litewka (Eds.), Recent Advances in Computational Mechanics, CRC Press, 2014, pp. 355–362. [48] W.J. Minkowycz, A. Haji-Sheikh, K. Vafai, On departure from local thermal equilibrium in porous media due to a rapidly changing heat source: the Sparrow number, Int. J. Heat Mass Transf. 42 (18) (1999) 3373–3385. [49] B.M. Kim, S.L. Jacques, S. Rastegar, S. Thomsen, M. Motamedi, Nonlinear finite-element analysis of the role of dynamic changes in blood perfusion and optical properties in laser coagulation of tissue, IEEE J. Sel. Top. Quantum Electron. 2 (4) (1996) 922–933. [50] E. Majchrzak, Ł. Turchan, A numerical analysis of heating tissue using the two-temperature model, in: Advanced Computational Methods and Experiments in Heat Transfer XIII, WIT Transactions on Engineering Sciences, 83, WIT Press, 2014, pp. 477–488. [51] H.W. Huang, Z.P. Chen, R.B. Roemer, A counter current vascular network model of heat transfer in tissues, J. Biomech. Eng. 118 (1996) 120–129. [52] A. Szasz, N. Szasz, O. Szasz, Oncothermia: Principles and Practices, Springer, 2011. [53] N. Afrin, J. Zhou, Y. Zhang, D.Y. Tzou, J.K. Chen, Numerical simulation of thermal damage to living biological tissues induced by laser irradiation based on a generalized dual phase lag model, Numer. Heat Transf. Part A: Appl. 61 (2012) 483–501. [54] C.M. Gardner, S.L. Jacques, A.J. Welch, Light transport in tissue: accurate expressions for one-dimensional fluence rate and escape function based upon Monte Carlo simulation, Lasers Surg. Med. 18 (1996) 129–138. [55] M. Diakite, H. Odéen, N. Todd, A. Payne, D.L. Parker, Toward real-time temperature monitoring in fat and aqueous tissue during magnetic resonance-guided high-intensity focused ultrasound using a three-dimensional proton resonance frequency T1 method, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 72 (2014) 178–187. [56] G. Liu, Q. Qin, K.W. Chan, Y. Li, J.W. Bulte, M.T. McMahon, P.C. van Zijl, A.A. Gilad, Non-invasive temperature mapping using temperature-responsive water saturation shift referencing (T-WASSR) MRI, NMR Biomed. 27 (3) (2014) 320–331. [57] S. Oh, Y.C. Ryu, G. Carluccio, C.T. Sica, C.M. Collins, Measurement of SAR-induced temperature increase in a phantom and in vivo with comparison to numerical simulation, Magn. Reson. Med. 71 (5) (2014) 1923–1931. [58] M.J. Thrippleton, J. Parikh, B.A. Harris, S.J. Hammer, S.I. Semple, P.J. Andrews, J.M. Wardlaw, I. Marshall, Reliability of MRSI brain temperature mapping at 1.5 and 3 T, NMR Biomed. 27 (2) (2013) 183–190. [59] W. Dai, H. Wang, P.M. Jordan, R.E. Mickens, A. Bejan, A mathematical model for skin burn injury induced by radiation heating, Int. J. Heat Mass Transf. 51 (2008) 5497–5510. [60] R. Singh, K. Das, S.C. Mishra, Laser-induced hyperthermia of nanoshell mediated vascularized tissue – A numerical study, J. Therm. Biol. 44 (2014) 55–62. [61] B. Mochnacki, A. Piasecka Belkhayat, Numerical modeling of skin tissue heating using the interval finite difference method, Mol. Cell. Biomech. 10 (3) (2013) 133–144.

´ Please cite this article as: M. Jasinski et al., Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation, Applied Mathematical Modelling (2015), http://dx.doi.org/10.1016/j.apm.2015.10.025