Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat flux conditions on skin tissue

International Communications in Heat and Mass Transfer 39 (2012) 121–130

Contents lists available at SciVerse ScienceDirect

International Communications in Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ichmt

Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat flux conditions on skin tissue☆ H. Ahmadikia a,⁎, R. Fazlali b, 1, A. Moradi b, 2 a b

Department of Mechanical Engineering, University of Isfahan, 81746-73441, Isfahan, Iran Department of Mechanical Engineering, Bu-Ali Sina University, Hamadan, Iran

a r t i c l e

i n f o

Available online 1 October 2011 Keywords: Skin tissue Bioheat transfer Laplace transform Hyperbolic Parabolic

a b s t r a c t In this article, the parabolic (Pennes bioheat equation) and hyperbolic (thermal wave) bioheat transfer models for constant, periodic and pulse train heat flux boundary conditions are solved analytically by applying the Laplace transform method for skin as a semi-infinite and finite domain. The bioheat transfer analysis with transient heat flux on skin tissue has only been studied by Pennes equation for a semi-infinite domain. For modeling heat transfer in short duration of an initial transient, or when the propagation speed of the thermal wave is finite, there are major differences between the results of parabolic and hyperbolic heat transfer equations. The non-Fourier bioheat transfer equation describes the thermal behavior in the biological tissues better than Fourier equation. The outcome of transient heat flux condition shows that by penetrating into the depths beneath the skin subjected to heat, the amplitude of temperature response decreases significantly. The blood perfusion rate can be predicted using the phase shift between the surface temperature and transient surface heat flux. The thermal damage of the skin is studied by applying both the parabolic and hyperbolic bioheat transfer equations. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Skin is the more extensive living organ of the human body. The skin generally is made of three layers: Epidermis, dermis, and hypodermis (subcutaneous tissue). It plays a variety of important roles including sensory, thermoregulation and host defense, etc. Advances in laser, microwave and similar technologies have led to recent developments in thermal treatment of disease and injured skin tissue such as skin cancer and skin burn. Various researchers have reported studying of the bioheat transfer in the biological tissues. The most famous of them are Pennes [1], Abramson [2], Lemons et al. [3], Weinbaum et al. [4] and Xu et al. [5]. They presented different bioheat transfer models for describing thermal behavior in the biological tissues. The parabolic heat transfer equation (Pennes bioheat equation) in the biological tissue was first introduced by Pennes in 1948 [1]. After the experimental observations [6], the wave-like behavior in the heat conduction has been argued. The modification of Fourier's law presented by Cattaneo [7] and Vernott [8] is a linear extension of the Fourier equation. Liu et al. [9] originally introduced a general form of the thermal wave model of

☆ Communicated by J. Taine and A. Soufiani ⁎ Corresponding author. Tel.: + 98 311 7934517; fax: +98 311 7932746. E-mail addresses: [email protected] (H. Ahmadikia), [email protected] (R. Fazlali), [email protected] (A. Moradi). 1 Tel.: + 98 811 8257410; fax: +98 811 8257400. 2 Tel.: + 98 861 4131213; fax: +98 861 4131217. 0735-1933/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2011.09.016

bioheat transfer (TWMBT) in living tissues. Liu and Lu [10] and Lu et al. [11,12] reported that some thermal wave's effects in bioheat transfer cannot be described by the Pennes' equation. The thermal wave model is generated by considering the relaxation time for heat flux in the Pennes equation. Other types of bioheat transfer model are Dual Phase Lag (DPL) models that are generated by involving relaxation times for the heat flux and temperature in the Pennes equation [13]. The Pennes equation (Fourier's law), assumes that any thermal distribution on the body is felt throughout the body instantaneously. This means that the propagation speed of thermal distribution is infinite. This assumption is reasonable in the majority of engineering applications. It fails in the particular thermal condition or heat transfer conduction media, where heat conduction shows non-Fourier feature or thermal wave phenomena, or hyperbolic heat conduction. Generally, where the materials have a great relaxation time (such as biological tissues), the thermal wave or non-Fourier phenomena are seen in heat transfer conduction. The non-Fourier heat transfer through the soft tissues such as skin has been studied by various researchers [14–16]. Mitra et al. [17] conducted different experiments on processed meat with different boundary conditions. They also observed wave-like phenomena in heat transfer conduction. Xu et al. [13] studied non-Fourier analysis of skin biothermomechanics. They applied the Pennes equation, thermal wave and three types of DPL models for investigating the conduction heat transfer through the skin tissues subjected to constant temperature boundary condition. They applied the finite difference scheme for solving energy equation in the multi

122

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

layer tissues and to validate numerical solution scheme applied the analytical results that were obtained by Liu et al. [18]. They observed substantial discrepancies among the Pennes, thermal wave and dual phase lag models, while different DPL models give similar predictions. Liu et al. [18] presented a new thermal wave aspect on burn evaluation of the skin. They used the thermal wave and Pennes bioheat transfer models for energy equation and solved the governing equations using separation of variables for constant surface temperature heating. Several researchers have studied the temperature response and blood perfusion during surface heating [19–29]. These researches have applied the Pennes bioheat transfer equation for the energy equations [22,23,26]. Shih et al. [26] solved the Pennes equation with sinusoidal surface heat flux on the skin as semi-infinite domain. They showed that their analytical expression is suitable for describing the transient temperature response of the tissue for the whole time domain from the beginning of the periodic oscillation to the final steady periodic oscillation. Liu [30] presented the analytical solution for the Pennes equation subjected to the sinusoidal surface heat flux on the skin tissue. He realized that by measuring the phase shifting between the surface temperature and sinusoidal heat flux, the blood perfusion rate could be obtained. Xu et al. [5] applied the Pennes equation for modeling heat transfer in skin tissue and employed the Green function method for solving the energy equation. In this study, the Fourier and non-Fourier equations with constant and transient heat flux conditions on the skin surface are studied analytically by using the Laplace transform method for semi-infinite and finite domains. The thermal damage is also studied here. In the related literatures, all reported analytical solutions for bioheat transfer in biological tissues used the Pennes equation as the energy equation. The novelty of this study is generating an exact analytical solution of the thermal wave bioheat transfer for a skin tissue with the constant and transient heat flux. One of the outcomes that have been overlooked more than others regarding bioheat transfer in a biological tissue is the burn time. Liu et al. [18] calculated it by numerical solution for both Pennes and thermal wave models. This is calculated analytically in this research. Studying thermal behaviors and thermal stresses of biological tissues with dual-phase-lag bioheat transfer model is the future work of the authors.

2.1. Fourier conduction heat transfer For bioheat transfer, the Pennes equation for modeling skin tissue heat transfer is known and expressed as [1]: ∂T ∂2 T ¼ k 2 þ ρb ϖb cb ðTa −T Þ þ qmet þ qext ∂t ∂x

2  ∂T ∂ T
þ ρt ct þ τq ρb ϖb cb þ ρb ϖb cb ðT−Ta Þ 2 ∂t ∂t   2 ∂ T ∂q ∂q ¼ k 2 þ qmet þ qext þ τq met þ τq ext ∂t ∂t ∂x

τq ρt ct

ð1Þ

where, ρt, ct and k are the density, specific heat and the thermal conductivity of skin tissue, respectively; ρb and cb are the density and specific heat of the blood, respectively and ϖb is the blood perfusion rate; T and Ta are skin tissue and blood temperatures, respectively; qmet and qext are the metabolic heat generation in skin tissue and the heat generated by other heating sources, respectively. In this study, qext is considered as zero. 2.2. Thermal wave model of the bioheat transfer (TWMBT)

3. The analytical solution of thermal wave bioheat transfer equation 3.1. The constant heat flux condition on a finite domain skin surface for thermal wave model In this study, the energy equation is solved analytically in a semiinfinite and finite domain of skin with three types of the constant, cosine and unit step function heat flux boundary conditions are considered. The geometrical problem and boundary conditions for finite and semi-infinite domains of skin tissue are shown in Fig. 1. Eq. (3) is the thermal wave bioheat transfer equation, which is applied to describe the wave-like heat transfer process through the skin tissue. The initial and boundary conditions can be written by: T ðx; 0Þ ¼ Ta ;

−k

∂qðx; t Þ ¼ −k∇T ðx; t Þ: ∂t

 ∂T  ¼ ∂x x¼0

∂T  ¼0  ∂t t¼0 q0 ;

ð4Þ

∂T  ¼ 0:  ∂x x¼L

rffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi Wb cb T−Ta W c W c kWb cb ; η ¼ b b t; Λ ¼ b b τq ; x; θ ¼ k q0 ρt ct ρt ct qmet ψ ¼ qffiffiffiffiffiffiffiffi q0 Wkb cb

ð6Þ

where, Wb = ρbϖb. Eq. (3) with respect to the constant qmet and qext = 0, can be rewritten in terms of dimensionless variables as: 2

Λ

2

∂ θ ∂θ ∂ θ þ ð1 þ ΛÞ þθ¼ 2þψ ∂η ∂η2 ∂ξ

ð7Þ

where, the dimensionless initial and boundary conditions are expressed by: θðξ; 0Þ ¼ 0;

∂θ  ¼0  ∂η η¼0

ð8Þ

 ∂θ  ¼ −1; ∂ξ ξ¼0

 ∂θ  pffiffiffiffiffiffiffi ¼ 0: W c ∂ξ ξ¼ bk b L

ð9Þ

By taking the Laplace transform of Eq. (7), we have the following:

2

ð2Þ

ð5Þ

For simplicity, the dimensionless variables are defined by:

Considering the concept of finite heat propagation velocity, Vernott [7] and Cattaneo [8] reported a modified unsteady heat conduction equation as:

qðx; t Þ þ τq

ð3Þ

where, τq = α/C 2 is the thermal relaxation time, α is the thermal diffusivity, and C is the thermal wave speed in the medium [17, 31].

ξ¼

2. Heat transfer models

ρt ct

Based on Eq. (2) for heat flux including the characteristic time τq as well as the Pennes' equation, a general form of the thermal wave model of bioheat transfer (TWMBT) in living tissues is expressed by [9,13]:

∂ θ ψ −βθ ¼ − ; s ∂ξ2

2

β ¼ Λs þ ð1 þ ΛÞs þ 1:

ð10Þ

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

123

a) finite domain with constant and pulse heat b) semi-infinite domain with constant and flux boundary conditions

pulse heat flux boundary conditions

c) finite domain with transient cosine heat flux boundary condition

Fig. 1. The geometrical problem finite and semi finite domains of skin tissue and three types of boundary conditions.

By taking the Laplace transform of the boundary conditions, we have:  ∂θ  1 ¼− ;  s ∂ξ ξ¼0

 ∂θ  pffiffiffiffiffiffiffi ¼ 0: W c ∂ξ ξ¼ bk b L

ð11Þ

Now, the partial differential equation is changed into an ordinary differential equation. The exact solution of Eq. (10) with respect to the boundary conditions of Eq. (11) can be obtained by:   qffiffiffiffiffiffiffiffi  pffiffiffiffi cosh β ξ− Wkb cb L ψ   þ : θ ðξ; sÞ ¼ pffiffiffiffi pffiffiffiffiqffiffiffiffiffiffiffiffi sβ Wb cb L s β sinh β k

ð12Þ

In this study, we use the inverse theorem for the inverse Laplace transform. For complementary information about inverse theorem, see Ref. [32]. At first, the inverse Laplace transform of the first term of Eq. (12) is calculated. The poles of this function are obtained by: ¼0 pffiffiffiffi 1 β ¼ 0→s ¼ −1; s ¼ − Λ8 rffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi > > β Wb cb L ¼ 0→s ¼ −1; s ¼ − 1 > > > k ffi Λ > rffiffiffiffiffiffiffiffiffiffiffi > >
W c pffiffiffiffi Wb cb > 2 > 2 b b 2 > rffiffiffiffiffiffiffiffiffiffiffiffi ! β L ¼ λ L ¼ −λm i⇒ Λs þ ð 1 þ ΛÞs þ 1 > m < pffiffiffiffi Wb cb k k vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! sinh β L ¼ 0→ u 2 u > k kλm > > −ð1 þ ΛÞ  tð1 þ ΛÞ2 −4Λ 1 þ > > 2 > Wb cb L > > > > > ⇒s ¼ S1m ; S2m ¼ > 2Λ : ; λm ¼ mπ; m ¼ 1; 2; …

s

ð13Þ

By calculating the residues at poles, the function θ1(ξ,η) is obtained as: sffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffi !! rffiffiffiffiffiffiffiffiffiffiffiffi ! Wb cb λ i k Wb cb S1 η L ξ− 2e m cosh m L ∞ L Wb cb k k rffiffiffiffiffiffiffiffiffiffiffiffi θ1 ðξ; ηÞ ¼ rffiffiffiffiffiffiffiffiffiffiffiffi ! þ ∑ Wb cb m¼1 Wb cb S1m L ð2ΛS1m þ Λ þ 1Þ coshðλm iÞ L sinh k k sffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffi !! λm i k Wb cb η S2m η ξ− 2e cosh L − ∞ L Wb cb k e−η −Λe Λ rffiffiffiffiffiffiffiffiffiffiffiffi þ ∑ rffiffiffiffiffiffiffiffiffiffiffiffi þ Wb cb Wb cb m¼1 L ðΛ−1Þ ð2ΛS2m þ Λ þ 1Þ coshðλm iÞ: S2m L k k cosh ξ−

ð14Þ

The inverse Laplace transform of the second term of Eq. (12) is obtained by:
 η ψ e−η −Λe−Λ θ2 ðξ; ηÞ ¼ ψ þ : ð15Þ Λ−1 Finally, by combining Eqs. (14) and (15), the closely formed function θ(ξ,η) is obtained by: θðξ; ηÞ ¼ θ1 ðξ; ηÞ þ θ2 ðξ; ηÞ:

ð16Þ

3.2. The constant heat flux condition on semi-infinite skin surface for thermal wave model The skin is assumed a semi-infinite medium. Thus, the boundary condition in Eq. (11) is substituted by Eq. (17) as follows: ∂θ  ¼0 ∂ξ ξ→∞

ð17Þ

124

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

where, the exact solution of Eq. (10) considering Eq. (17) for the constant heat flux condition at the surface tissue is obtained by:

Eq. (3) with respect to the constant qmet and qext = 0 can be rewritten in terms of dimensionless variables that give:

pffiffiffi 1 − βξ ψ θ ðξ; sÞ ¼ pffiffiffiffi e þ : sβ s β

Λ

ð18Þ

The Laplace transform's tables are used for calculating the inverse Laplace transform of function θ ðξ; sÞ. The inverse Laplace transform of the first term of Eq. (18) is calculated by: θðξ; ηÞ ¼ L

−1 

22 0 pffiffiffi 13 − βξ η6 −1 e θ ðξ; sÞ ¼ ∫ 44L @ pffiffiffiffi A5 0 β

η→v

3

    −1 1 −1 ψ 7 × L 5dv þ L s η→η−v sβ

ð19Þ where, 0 0 pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 0 pffiffiffi 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 η − βξ − ΛðsþΛ1Þðsþ1Þξ −Λ − Λξ sðsþ1−Λ1 Þ e e e L−1 @ −1 B −1 Be C C pffiffiffiffi A ¼ L @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ¼ pffiffiffiffi ×L @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



ffi A: β Λ Λ s þ 1 ðs þ 1Þ s s þ 1− 1 Λ

Λ

′

 ∂2 θ
∂2 θ ′ ∂θ þ 1 þ c1 Λ þ c1 θ ¼ 2 þ φ: 2 ∂η ∂η ∂ξ

ð26Þ

The dimensionless form of surface boundary condition is expressed by:  ∂θ  ¼ − cosðηÞ: ∂ξ ξ¼0

ð27Þ

Other dimensionless boundary and initial conditions are the same as that in Section 3.1. By taking Laplace transform of Eq. (26) and considering the initial condition, we have the following:
 ∂2 θ φ ′ ′ ′ 2 ′ −β θ ¼ − ; β ¼ Λ s þ 1 þ c1 Λ s þ c1 : 2 s ∂ξ

ð28Þ

ð20Þ The boundary condition of Eq. (27) in Laplace domain is expressed To obtain the inversion of Eq. (20), the following formula from the Laplace table is applied [39]. 0 pffiffiffiffiffiffiffiffiffiffi 1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi −k sðsþaÞ 1 1 −1 @e − at pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ¼ e 2 Io a t 2 −k2 U ðt−kÞ; k≥0 L 2 sðs þ aÞ

ð21Þ

where, U is the unit step function. The inverse Laplace transform of Eq. (20) is obtained as: 0 pffiffiffi 1 η   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
− βξ − pffiffiffiffi e Λ −1 1−1 η 1 1 L−1 @e pffiffiffiffi A ¼ pffiffiffiffi e 2ð ΛÞ I0 1− η2 −ξ2 Λ U η−ξ Λ : 2 Λ β Λ

by:  ∂θ  s ¼− : ∂ξ ξ¼0 1 þ s2

ð29Þ

The exact solution of Eq. (28) with respect to the boundary and initial conditions is obtained by:
pffiffiffiffiffi pffiffiffi  s cosh β′ ξ− ω αL φ
pffiffiffiffiffipffiffiffi  þ ′ : θ ðξ; sÞ ¼ pffiffiffiffiffi

2 ′ sβ 1þs β′ ω β sinh L α

ð30Þ

ð22Þ By substituting Eqs. (15) and (22) for Eq. (19), the closely formed function θ(ξ,η) is obtained as: θðξ; ηÞ ¼ ∫

η 0


η ! v qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
  ψ e−η −Λe−Λ pffiffiffiffi e−Λ −12ð1−Λ1Þv 1 1 pffiffiffiffi e 1− v2 −ξ2 Λ U v−ξ Λ dv þ ψ þ I0 : 2 Λ Λ−1 Λ

At first, the inverse Laplace transform of the first term of Eq. (30) is calculated. The poles for this function are obtained by similar method that was described in Eq. (13) and they become:

ð23Þ s ¼ S1n ; S2n ¼

In this study, for solving the thermal wave and Pennes bioheat transfer model with transient heat flux condition on the skin surface, the cosine flux is considered on the skin surface. The geometrical problem and boundary conditions are shown in Fig. 1(c). Here, the boundary condition on the skin surface only should be adapted according to cosine heat flux condition; thus, it is expressed by: ∂T ∂x

j

x¼0

¼ q0 cosðωt Þ

ð24Þ

where, q0 is the heat flux on the skin surface, and ω is the heating frequency. Other boundary conditions are same as that in Section 3.1. Here, the dimensionless variables are defined by: ξ¼

rffiffiffiffi rffiffiffiffi ω kðT−Ta Þ ω W c qmet ′ ffiffiffiffiffiffiffi : x; θ ¼ ; c ¼ b b ; η ¼ ωt; Λ ¼ ωτq ; φ ¼ q α q0 α 1 ρt ct ω q ρcω 0

2Λ′

; λn ¼ nπ; n ¼ 1; 2; …

ð31Þ

3.3. The transient heat flux condition on a layer of skin surface for thermal wave model

−k

1 1 ; s ¼ −c1 ; s ¼ − ′ ; s ¼ −c1 ; Λ Λ′ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 2
 u u

αλn − 1 þ c1 Λ′  t 1 þ c1 Λ′ 2 −4Λ′ c1 þ 2 ωL

s ¼ i ; s ¼ −

k

ð25Þ

With calculating of the residues at poles, the function θ1(ξ,η) is obtained as: η rffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffi 

− ω α ′ eiη cosh 1 þ iΛ′ ði þ c1 Þ ξ− L e Λ′ Λ α ω  θ1 ðξ; ηÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffirffiffiffiffi  þ






ω L 1 þ Λ′ 2 1−Λ′ c1 1 þ iΛ′ ði þ c1 Þ L 2 1 þ iΛ′ ði þ c1 Þ sinh α qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffi 

ω e−iη cosh 1−iΛ′ ðc1 −iÞ ξ− L −c η c1 e 1 α þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffirffiffiffiffi  þ rffiffiffiffi





ω ω ′ ′ 2 1−iΛ ðc1 −iÞ L L 1 þ c1 −1 þ Λ′ c1 2 1−iΛ ðc1 −iÞ sinh 1 α 0 rffiffiffiffi rffiffiffiffi α   λ i α ω ξ− L 2S1n eS1n η cosh n C ∞ B L ω α C B rffiffiffiffi þ∑B C

′ A ω n¼1@ ′ 2 L 1 þ S1n 2Λ S1n þ c1 Λ þ 1 coshðλn iÞ α 1 0 rffiffiffiffi rffiffiffiffi   λ i α ω ξ− L 2S2n eS2n η cosh n C ∞ B L ω α C B þ ∑ B rffiffiffiffi C:

′ A ω n¼1@ ′ 2 L 1 þ S2n 2Λ S2n þ c1 Λ þ 1 coshðλn iÞ α

ð32Þ

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

The inverse Laplace transform of the second term of Eq. (30) is obtained by: ! η 1 e−Λ′ Λ′ e−c1 η

: − ′ þ ′ θ2 ðξ; ηÞ ¼ ϕ c1 Λ c1 −1 c1 Λ c1 −1

ð33Þ

125

By substituting Eqs. (15) and (22) for Eq. (40), the closely formed function θ(ξ,η) is obtained as: 0 1 1 − η ψ@e−η −Λe Λ A θðξ; ηÞ ¼ ψ þ 2

Λ−1

0

1 3   v 1 1 − qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   pffiffiffiffi7 Be Λ − 2 1− Λ v C
1 1 2 2 7 ffi e v −ξ Λ C I0 ðU ðη−vÞ−U ðη−v−Λi ÞÞB 1− þ∫ 6 @ pffiffiffi AU v−ξ Λ 5dv: 04 2 Λ Λ η6

Finally, by adding Eqs. (32) and (33), the close form of function θ(ξ,η) is obtained by: θðξ; ηÞ ¼ θ1 ðξ; ηÞ þ θ2 ðξ; ηÞ:

ð41Þ

ð34Þ 4. The analytical solution of the Pennes bioheat transfer equation (PBTE)

3.4. The solution of the pulse train heat flux condition on a semi-infinite skin surface for thermal wave model Usually, the common burns occur when the skin is exposed to the high intensity heat flux in a short time. Thus, to investigate the thermal behavior of skin is supposed the high intensity heat flux in a pulse time τi. Therefore, the boundary condition on skin surface is defined by: ∂T −k ∂x

j

x¼0

¼ q0 ðU ðt Þ−U ðt−τi ÞÞ:

ð35Þ

Wb cb τ: ρt ct i

b b

ð37Þ

b b

 ∂θ  1−e−sΛ i : ¼  s ∂ξ ξ¼0

ð38Þ

The dimensionless temperature in the Laplace domain for pulse train heat flux condition on the skin surface and semi-infinite domain is obtained by: −sΛ i

1−e pffiffiffiffi s β

pffiffiffi

−

e

βξ

þ

ψ : sβ

ð39Þ

The inverse Laplace transform of Eq. (39) is calculated by similar method for calculating the inverse Laplace transform of Eq. (18), described in Eqs. (19)–(22), and it is presented as follows:

−1 

2

∂ θ ψ −βP θ ¼ − ; s ∂ξ2

βP ¼ s þ 1:

ð43Þ

22

0 pffiffiffi 13 − βξ η6 −1 e θ ðξ; sÞ ¼ ∫ 44L @ pffiffiffiffi A5 0 β

" −1

× L η→v

−sΛ i

1−e s


pffiffiffiffiffiffi
qffiffiffiffiffiffiffiffi  cosh βP ξ− Wkb cb L ψ θ ðξ; sÞ ¼ pffiffiffiffiffiffi :
pffiffiffiffiffiffiqffiffiffiffiffiffiffiffi  þ Wb cb sβ P s βP sinh βP L k

ð44Þ

The poles of the first term of Eq. (44) are obtained by:

The boundary condition of Eq. (37) in the Laplace domain is expressed by:

θðξ; ηÞ ¼ L

ð42Þ

The exact solution of Eq. (43) with respect to the boundary condition of Eq. (11) is obtained by:

      ∂θ  ρt ct ρt ct η −U ðη−Λi Þ : − ξ¼0 ¼ U W c W c ∂ξ

θ ðξ; sÞ ¼

2

∂θ ∂ θ þ θ ¼ 2 þ ψ: ∂η ∂ξ

ð36Þ

Using Eqs. (36) and (6), the dimensionless boundary condition of Eq. (35) is rewritten as:

−

The Pennes bioheat transfer equation as given in Eq. (1) is applied to describe parabolic heat transfer process through the skin tissue. Substituting the dimensionless parameters for Eq. (1), hence we get:

Eqs. (4) and (5) indicate the boundary conditions. Taking Laplace transform of Eq. (42) and considering the initial conditions, we have:

The skin is assumed a semi-infinite domain. For simplicity the dimensionless parameter, Λi is defined as: Λi ¼

4.1. The constant heat flux condition on a finite domain skin surface for Pennes equation

!# η→η−v

3

  −1 ψ 7 : 5dv þ L sβ

ð40Þ

s ¼ 0; s ¼ −1; s ¼ −1; s ¼ SP ¼ −1−

kλm 2 : Wb cb L2

ð45Þ

The poles of the second term are s = 0, s = −1. Using inverse theorem method, the closely formed function θ(ξ,η) is obtained by:
qffiffiffiffiffiffiffiffi  −η cosh ξ− Wkb cb L e θðξ; ηÞ ¼
qffiffiffiffiffiffiffiffi  − qffiffiffiffiffiffiffiffi Wb cb Wb cb sinh L L k k
qffiffiffiffiffiffiffiffi
qffiffiffiffiffiffiffiffi  SP η λm i Wb cb k cosh L L ∞ 2e W c ξ− k −η qffiffiffiffiffiffiffiffi b b þ∑ : þ ψ 1−e W c b b m¼1 coshðλm iÞ SP L k

ð46Þ

4.2. The constant heat flux condition on a semi-infinite skin surface for Pennes equation For skin as a semi-infinite domain, the exact solution of Eq. (43) with respect to Eq. (17) is obtained by: pffiffiffiffi 1 ψ − βP ξ θ ðξ; sÞ ¼ pffiffiffiffiffiffi e þ : sβP s βP

ð47Þ

126

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

4.4. The solution of the pulse train heat flux condition on a semi-infinite skin surface for Pennes equation

The inverse Laplace transform of Eq. (47) is calculated by: pffiffiffiffi " # 1 e− βP ξ ψ pffiffiffiffiffiffi þ s sβP βP 2" 3 pffiffiffiffi !#

  − βP ξ η −1 e −1 1 4 5dv pffiffiffiffiffiffi ¼ ∫0 L × L s η→η−v βP η→v   ψ −1 þL sβP

θðξ; ηÞ ¼ L

−1

Here, the boundary and initial conditions are same as that in Section 3.4. The exact solution of Eq. (43) considering the mentioned boundary conditions is obtained by:

ð48Þ

−1

pffiffiffiffi !

−

e

pffiffiffiffi 1−e−sΛ i − βP ξ ψ pffiffiffiffiffiffi e þ : sβP s βP

ð56Þ

By similar method presented in Eq. (40) and applying Eq. (49), the inverse Laplace transform of Eq. (56) is calculated as follows:

where,

L

θ ðξ; sÞ ¼

βP ξ

pffiffiffiffiffiffi βP

0 pffiffiffiffiffiffiffiffiffi 1 2 pffi ! −η−ξ − ðsþ1Þξ −ξ s 4η e −1 @e −η −1 e A p ffiffiffiffiffiffiffiffiffiffiffi p ffiffi ¼L ¼ pffiffiffiffiffiffi : ¼e L πη s sþ1

ð49Þ

After substituting Eq. (49) for Eq. (48), Eq. (48) is calculated and the closely formed function θ(ξ,η) is obtained by: 2 −v−ξ 4v

ηe −η : θðξ; ηÞ ¼ ∫0 pffiffiffiffiffiffi dv þ ψ 1−e πv

2 3 2 −v−ξ 4v e −η : θðξ; ηÞ ¼ ∫0 4 pffiffiffiffiffiffi ðU ðη−vÞ−U ðη−v−Λi ÞÞ5dv þ ψ 1−e πv η

ð50Þ

4.3. The transient heat flux condition on a finite domain skin surface for Pennes equation

4.5. Thermal damage The burn evaluation is one of the most important characteristics in the bioengineering science in skin tissue. It is proven that the thermal damage begins when the basal layer temperature rises to 44 °C [33]. The basal layer is located between the epidermis and dermis. Moritz and Henriques [34, 35] first proposed the quantitative of thermal damage evaluation because the tissue damage could be represented as an integral of a chemical process rate as: t

Ω ¼ ∫0 A expð−Ea =RT Þdt Here, the boundary and initial conditions are same as that in Section 3.1. Considering the dimensionless variables Eq. (25), the Pennes equation is rewritten as: ∂θ ∂2 θ þ c1 θ ¼ 2 þ φ: ∂η ∂ξ

ð51Þ

ð57Þ

ð59Þ

where, A is a material parameter equivalent to a frequency factor, Ea is the activation energy, and R = 8.313 J/mol·K is the universal gas constant. The constants A and Ea are obtained experimentally. The values of these constant parameters for the skin tissue in different parts of the human body are presented by Xu et al. [36]. By fitting these experimental data, a linear relation was observed between Ea and Ln(A), and the relation is expressed as [36]:

Taking Laplace transform from Eq. (51) we have: Ea ¼ 21149:324 þ 2688:367 lnðAÞ: 0 ∂2 θ φ −βP θ′ ¼ − ; 2 s ∂ξ

0

βP ¼ s þ c1 :

ð52Þ

Substituting boundary condition for Eq. (52), the dimensionless temperature of skin in the Laplace domain is obtained by: θ ðξ; sÞ ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi

s cosh s þ c1 ξ− ω φ αL

pffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffi ffi p ffipffiffiffi þ 0 : 1 þ s2 s þ c1 sinh s þ c1 ω sβP αL

ð53Þ

Applying the inverse theorem method, the closely formed function θ(ξ,η) is obtained by:   rffiffiffiffi    rffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ω ω iη −i þ c1 ξ− e cosh i þ c1 ξ− L L α α rffiffiffiffi  þ rffiffiffiffi  θðξ; ηÞ ¼   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ω −i þ c1 i þ c1 L L 2 −i þ c1 sinh 2 i þ c1 sinh α α 0   rffiffiffiffi  rffiffiffiffi 1 ω ω SP η 2SP e cosh λn i ξ− L L −c1 η −c η ∞ B c1 e 1 α α C B C φ 1−e rffiffiffiffi þ ∑ B rffiffiffiffi − Cþ

ω

A c1 ω n¼1@ 1 þ c21 L 1 þ SP 2 coshðλn iÞ L α α −iη

e

cosh

=

ð54Þ

where, ! 2 αλn ; λn ¼ nπ; n ¼ 1; 2; … SP ¼ − c1 þ ωL2

ð55Þ

ð60Þ

For the first and second-degree burns, T in Eq. (59) is the basal layer temperature, while for the third-degree burns, it is the temperature at the interface between dermis and subcutaneous layers (fat). The first and second-degree burns occur when both the conditions T N 44 °C and Ω N 0.53 are satisfied in the basal layer. When Ω = 1, the second-degree burn occurs [18]. 5. Results and discussions In this study, the temperature distribution for both thermal wave (hyperbolic) and Pennes bioheat equation (parabolic) models are studied for three different surface boundary conditions. These boundary conditions are: constant heat flux, transient heat flux and pulse train heat flux conditions. The adiabatic condition is used for two different finite and semi-infinite domains. In these cases, the parabolic and hyperbolic analytical results are compared. The blood and skin tissue properties are considered same as the properties that were reported by Xu et al. [13,40]. The arterial blood and the blood perfusion rate are considered as Ta = 37 °C [37] and Wb = ρbϖb = 0.5 kg/m 3/s [26] respectively. The relaxation time is considered τq = 16 s here [17]. The temperature distribution for both the parabolic and hyperbolic bioheat transfer models for skin as finite and semi-infinite domains at the basal and the interface between dermis and fat layer (DF interface) is shown in Figs. 2 and 3, respectively. In these figures, the constant heat flux condition, q0 = 5000 W/m 2 on the skin surface is

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

Fig. 2. The comparison between temperature response for the skin as a finite and semiinfinite domain at the basal layer for both parabolic and hyperbolic models.

127

Fig. 4. Temperature response for both parabolic and hyperbolic models for different heat flux intensities and fire time duration at the basal layer.

considered. As shown in Figs. 2 and 3, there are good agreements for finite and semi-infinite domains in both the parabolic and hyperbolic models up to t = 200 s; but after this time these figures that relate to the skin as a finite and semi-infinite domain become separated. This discrepancy is related to the bottom surface boundary condition. When the skin subjects to heat flux for short time duration, supposing skin as a semi-infinite domain is a good approximation that this supposition would be an easier problem. In addition, the parabolic result is compared in Fig. 2 with that obtained by Shih et al. [26] for the same problem, where they applied a Pennes equation to study heat transfer for skin tissue as a semi-infinite domain but they ignored the term of metabolic heat generation. Good agreement is observed in Fig. 2. The discrepancy in parabolic and hyperbolic temperature response to the pulse train with different intensity and pulse time is investigated at the basal layer (see Fig. 4). In Fig. 2, the solid and dash lines represent the temperature distribution for q0 = 83,200 W/m 2 and τi = 1 s for parabolic and hyperbolic models, respectively. The dash-dot-dot and dash-dot lines also represent the temperature distribution for q0 = 5000 W/m 2 and τi = 20 s for parabolic and hyperbolic models, respectively. It appears that when the high intensity heat flux strikes to skin for a short duration of time, large discrepancies become evident amongst the predictions of the Pennes model

and thermal wave models. Conversely, when the heat flux is not intensive and the pulse time is increased, the discrepancies between the Pennes and the thermal wave models become less. The comparison between parabolic and hyperbolic models at ED and DF interfaces is investigated when the pulse train heat flux q0 = 83,200 W/m 2 in the pulse time τi = 3 s strikes the skin surface (see Fig. 5). This condition is similar to the one where the skin is exposed to the propane gas flame for 3 s [38]. The relaxation time is observed in the thermal wave model clearly, because it takes 28 s for the temperature to rise in order to be felt at DF interface, while it takes only 3 m/s for the same to happen when skin is exposed to propane gas flame (see Fig. 5). To investigate the transient heat flux condition on the skin surface, the cosine heat flux with amplitude q0 = 5000 W/m 2 is considered [26]. The temperature response on the skin surface (x = 0) subjected to a cosine heat flux with different frequencies for thermal wave model is shown in Fig. 6. This figure illustrates that with an increase in heat frequency, the temperature amplitude is decreased and the response cyclic time becomes smaller. The temperature variation of skin surface with time is shown in Fig. 7 for heat flux frequency of 0.05. In this figure, the hyperbolic model results are plotted for two different relaxation times and the parabolic model result is shown

Fig. 3. The comparison between temperature response for the skin as a finite and semiinfinite domain at the DF interface for both parabolic and hyperbolic models.

Fig. 5. The variation of temperature with time for both parabolic and hyperbolic models at the ED and DF interface when the skin subjects to the propane gas flame for 3 s.

128

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

Fig. 6. Temperature response of thermal wave model at the skin surface with time for different heat flux frequencies when q0 = 5000 W/m2.

Fig. 8. Temperature response of the parabolic and hyperbolic models at the ED and DF interface when q0 = 5000 W/m2 and ω = 0.05.

in Fig. 7. The parabolic results are like the related data presented by Shih et al. [26] with the same thermal properties. It is clear from Fig. 7 that when the relaxation time becomes smaller, the discrepancy between thermal wave model and the Pennes equation decreases. This phenomenon validates the accuracy of the analytical solution for the thermal wave model. The temperature response at the basal layer and DF interface for both the parabolic and hyperbolic models is shown in Fig. 8, where the heat flux amplitude and frequency are q0 = 5000 W/m 2 and ω = 0.05, respectively and relaxation time is 16 s. It is evident from Fig. 8 that with pervading to the skin depth the amplitude of temperature decreases strongly. The dimensionless temperature response for both the parabolic and hyperbolic models and dimensionless heat flux variation on the skin surface when q0 = 5000 W/m 2 and ω = 0.05 are shown in Fig. 9 for the initial period η = 0 toη = 8. This period has an equivalent time interval from 0 to 160 s. As mentioned before, a blood perfusion rate can be calculated by measuring the phase shift between the heat flux and temperature response at the skin surface. In the first cyclic period because there are some oscillations, the measurement of first cyclic period cannot be a correct estimate for the blood perfusion rate, since it has an unstable parameter; therefore the second cyclic period is necessary. In this study, the second cyclic period occurs in the dimensionless time

η = 7 that is equal t = 140 s, so reaching to the second cyclic period needs 2 min while with decreasing of the heat flux frequency this time becomes more. This time for heat flux frequency ω = 0.01 would be about 10 min. It is clear from Fig. 9 that there is a phase shift between the parabolic and hyperbolic models as well. Fig. 10 shows the dimensionless temperature response for hyperbolic model and heat flux variation for different values of blood perfusion rate for the dimensionless time from η = 0 to η = 20 when q0 = 5000 W/m 2 and ω = 0.05. For the smaller values of blood perfusion rate c1 = 0.01 and c1 = 0.1, the temperature distributions coincide with each other. This means that a high accuracy device is required in order to detect the phase shift between these cases. As shown in Fig. 10, for c1 = 1 the phase shift is detected clearly because the greater amounts of blood perfusion rate can carry away more heat in a short time, furthermore, the phase difference between temperature response and heat flux decreases by an increase in the blood perfusion rate. These dimensionless blood perfusion rates c1 = 0.01 and c1 = 0.1 that are applied in this study have an equivalent blood perfusion rate from 0.568 to 5.682 kg/m 3/s, when the cosine heat flux frequency is ω = 0.05. In different heat flux frequencies, these intervals for the blood perfusion rate change for medical applications.

Fig. 7. The variation of skin surface temperature with time for both parabolic and hyperbolic models when q0 = 5000 W/m2 and ω = 0.05.

Fig. 9. The dimensionless temperature of the parabolic and hyperbolic models and heat flux variation with dimensionless time at the skin surface when ω = 0.05.

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

Fig. 10. The variation of dimensionless temperature of thermal wave model with time for different values of dimensionless blood perfusion rate when ω = 0.05.

The skin burn evolution is one of the most important problems that are investigated and studied in the biomedical engineering field. The process of skin burn takes place quickly and the resulted injury often occurs at the early stage of heating. In this research, times of the first and the second-degree burns are obtained when skin tissue is exposed to the propane gas flame for 3 s. The thermal damage variation with time at the ED interface is shown in Fig. 11. It is evident from Fig. 11 that the burn time for Pennes model is about 0.3 s, while the thermal wave model predicts this time to be about 3 s. This discrepancy appears when the flash fire duration is short. There is a significant difference between the anticipated burn times of the Pennes and the thermal models. For small heat flux when the burn time duration is long, the observed difference between two models becomes negligible. Therefore, for an accurate burn time prediction under fast heating conditions with high power, results predicted by the thermal wave model can be more realistic once the value τq is determined. The burn times for these conditions were obtained numerically by Liu et al. [18] as being 0.26 and 2.89 s for the Pennes and thermal wave models, respectively. The minor discrepancy between present results and Liu's results [14] is due to the different values assigned to the skin and blood properties.

Fig. 11. The variation of thermal damage with time for both parabolic and hyperbolic models at the ED interface when the skin subjects to the propane gas flame for 3 s.

129

Fig. 12. Temperature response of thermal wave model at the ED and DF interface for different values of blood perfusion rate when the constant heat flux is q0 = 2000 W/m2.

When the high intensity heat flux strikes, the skin surface in a short time period due to the intensity of the heating process the blood perfusion rate has no effect on the temperature response. While the heat flux is low and duration time is long, increasing the rate of blood perfusion leads to a slow skin temperature increase. The effect of increasing blood perfusion rate on the temperature response that was obtained by the thermal wave model at the basal layer and DF interface for continuous heat flux q0 = 2000 W/m 2 is shown in Fig. 12 for different values of blood perfusion rates. With an increase in blood perfusion rate, the large amounts of heat can be carried away through the rate of blood perfusion. 6. Conclusions In this study, the one-dimensional bioheat transfer analysis in the skin tissue was studied. The effect of the constant, cosine and pulse train heat flux condition at the skin surface on the transient temperature response in a one-dimensional tissue was investigated. The exact solutions of the Pennes and thermal wave models with three previously mentioned heat flux condition on the skin surface were performed in this study via the Laplace transform method. The analytical expressions obtained in this work are useful for describing the tissue temperature response in a one-dimensional domain, especially for the thermal wave model, because the problems in the previous researches in this literature were solved numerically. The analytical results demonstrate that for fast heating with high power, large discrepancies are noticed amongst the anticipation of the Pennes model and thermal wave model. While for a low heat flux when the flash fire duration time is short, although the difference in the predicted temperatures from two models will finally coincide at the end, this difference at the early stage is distinctive. The analytical results for cosine heat flux show that using the phase shift between the temperature and heat flux oscillation wave can measure the blood perfusion rate. Obtaining the phase shift between the temperature and heat flux is not used in the first cyclic period because of its intrinsic unstable character. The phase shift between cosine heat flux and temperature also becomes shorter with an increase in blood perfusion rate. When the high intensity heat flux contacts the skin surface in the short time, the blood perfusion rate has no effect on the temperature response while for the low heat flux the contact to the skin surface for a long time is conversely. The analytical results for a low continuous heat flux describe that when the blood perfusion rate increases, the increasing rate of temperature in the skin tissue decreases. The thermal damage results also present the different

130

H. Ahmadikia et al. / International Communications in Heat and Mass Transfer 39 (2012) 121–130

values of burn evolution for the Pennes and thermal wave models. The presented analytical results for burn time have an excellent agreement with numerical results that were reported by Liu et al. [18] in this literature. References [1] H.H. Pennes, Analysis of tissue and arterial blood temperature in the resting forearm, Journal of Applied Physiology 1 (1948) 93–122. [2] D.I. Abramson, Circulation in the Extremities, Academic press, New York, 1967. [3] D.E. Lemons, S. Chien, L.I. Crawshaw, S. Weinbaum, L.M. Jiji, Significance of vessel size and type in vascular heat transfer, American Journal of Physiology 253 (1987) 128–135. [4] S. Weinbaum, L.M. Jiji, D.E. Lemons, Theory and experiment for the effect of vascular microstructure on surface tissue heat transfer—part I: anatomical foundation and model conceptualization, Journal of Biomechanical Engineering 106 (4) (1984) 321–330. [5] F. Xu, T.J. Lu, K.A. Seffen, Biothermomechanics of skin tissues, Journal of the Mechanics and Physics of Solids 56 (2008) 1852–1884. [6] V. Peshkov, Second sound in helium II, Journal of Physics VIII 381 (1944). [7] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, Compte Rendus 247 (1958) 431–433. [8] P. Vernotte, Les paradoxes de la theorie continue de l'equation de la chaleur, Compte Rendus 246 (1958) 3154–3155. [9] J. Liu, Z. Ren, C. Wang, Interpretation of living tissue's temperature oscillations by thermal wave theory, Chinese Science Bulletin 40 (1995) 1493–1495. [10] J. Liu, W.Q. Lu, Dual reciprocity boundary element method for solving thermal wave model of bioheat transfer, Space Medicine & Medical Engineering 10 (6) (1997) 391–395. [11] W.Q. Lu, J. Liu, Y. Zeng, Simulation of the thermal wave propagation in biological tissues by the dual reciprocity boundary element method, Engineering Analysis with Boundary Elements 22 (3) (1998) 167–174. [12] W.Q. Lu, J. Liu, Y. Zeng, Extension of the dual reciprocity boundary element method to simulate thermal wave propagation in biological tissue, Journal of Engineering Thermophysics 19 (6) (1998) 728–731. [13] F. Xu, K.A. Seffen, T.J. Lu, Non-Fourier analysis of skin biothermomechanics, International Journal of Heat and Mass Transfer 51 (2008) 2237–2259. [14] A. Grassmann, F. Peters, Experimental investigation of heat conduction in wet sand, Heat and Mass Transfer/Waerme- und Stoffuebertragung 35 (4) (1999) 289–294. [15] H. Herwig, K. Beckert, Experimental evidence about the controversy concerning Fourier or non-Fourier heat conduction in materials with a non-homogeneous inner structure, Heat and Mass Transfer/Waerme- und Stoffuebertragung 36 (5) (2000) 387–392. [16] H. Herwig, K. Beckert, Fourier versus non-Fourier heat conduction in materials with a non-homogeneous inner structure, Journal of Heat Transfer, Transactions of the ASME 122 (2) (2000) 363–365. [17] K. Mitra, S. Kumar, A. Vedavarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, Journal of Heat Transfer, Transactions of the ASME 117 (3) (1995) 568–573. [18] J. Liu, X. Chen, L.X. Xu, New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating, IEEE Transaction on Biomedical Engineering 46 (4) (1999) 420–428.

[19] J. Liu, Uncertainty analysis for temperature prediction of biological bodies subject to randomly spatial heating, Journal of Biomechanics 34 (2001) 1637–1642. [20] Z.S. Deng, J. Liu, Mathematical modeling of temperature mapping over skin surface and its implementation in thermal disease diagnostics, Computers in Biology and Medicine 34 (2004) 495–521. [21] E.Y.K. Ng, L.T. Chua, Comparison of one-and two-dimensional programs for predicting the state of skin burns, Burns 28 (2002) 27–34. [22] C. Chen, L.X. Xu, Tissue temperature oscillations in an isolated pig kidney during surface heating, Annals of Biomedical Engineering 30 (2002) 1162–1171. [23] E.H. Liu, G.M. Saidel, H. Harasaki, Model analysis of tissue responses to transient and chronic heating, Annals of Biomedical Engineering 31 (2003) 1007–1048. [24] R. Chopra, J. Wachsmuth, M. Burtnyk, M.A. Haider, M.J. Bronskill, Analysis of factors important for transurethral ultrasound prostate heating using MR temperature feedback, Physics in Medicine and Biology 51 (2006) 827–844. [25] S.C. Jiang, N. Ma, H.J. Li, X.X. Zhang, Effects of thermal properties and geometrical dimensions on skin burn injuries, Burns 28 (2002) 713–717. [26] T.C. Shih, P. Yuan, W.L. Lin, H.S. Kou, Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface, Medical Engineering & Physics 29 (2007) 946–953. [27] J. Liu, Preliminary survey on the mechanisms of the wave-like behaviors of heat transfer in living tissues, Forschung im Ingenieurwesen 66 (2000) 1–10. [28] K.C. Liu, H.T. Chen, Investigation for the dual phase lag behavior of bio-heat transfer, International Journal of Thermal Sciences 49 (2010) 1138–1146. [29] F. Xu, T.J. Lu, K.A. Seffen, Thermally-induced change in the relaxation behavior of skin tissue, Journal of Biomechanical Engineering 131 (2009) 1–10. [30] J. Liu, Estimation of blood perfusion using phase shift in temperature response to sinusoidal heating at the skin surface, IEEE Transactions on Biomedical Engineering 46 (1999) 1037–1042. [31] D.Y. Tzou, Macro- to Micro-scale Heat Transfer: The Lagging Behavior, Taylor and Francis, Washington, DC, 1997. [32] V.C. Arpaci, Conduction Heat Transfer, Addisson Wesley Publication, 1966. [33] A.M. Stoll, L.C. Greene, Relationship between pain and tissue damage due to thermal radiation, Journal of Applied Physiology 14 (1959) 373–382. [34] A.R. Moritz, F.C. Henriques, Study of thermal injuries II. The relative importance of time and source temperature in the causation of cutaneous burns, American Journal of Pathology 23 (1947) 695–720. [35] F.C. Henriques, A.R. Moritz, Studies of thermal injury, 1. The conduction of heat to and through skin and the temperatures attained therein, a theoretical and an experimental investigation, American Journal of Pathology 23 (1947) 531–549. [36] F. Xu, K.A. Seffen, T.J. Lu, Temperature dependent mechanical behaviors of skin tissue, IAENG International Journal of Computer Science 35 (2008) 1–13. [37] J. Zhou, J.K. Chen, Y. Zhang, Dual phase lag effects on thermal damage to biological tissues caused by laser irradiation, Journal of Computers in Biology and Medicine 39 (2009) 286–293. [38] W.P. Bechnke, Predicting flash fire protection of clothing from laboratory tests using second-degree burn to rate performance, Fire and Materials 8 (1984) 57–63. [39] M. Abramowitz, A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972, pp. 1026–1027. [40] F. Xu, T. Lu, Introduction to skin biothermomechanics and thermal pain, science press Beijing and Springer-Verlag Berlin Heidelberg, 2011.