Computers in Biology and Medicine 62 (2015) 65–75

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Heat treatment modelling using strongly continuous semigroups Alaeddin Malek n, Ghasem Abbasi Department of Applied Mathematics, Faculty of Mathematical Sciences, Tarbiat Modares University, P.O. Box 14115-134, Tehran, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 9 September 2014 Accepted 31 March 2015

In this paper, mathematical simulation of bioheat transfer phenomenon within the living tissue is studied using the thermal wave model. Three different sources that have therapeutic applications in laser surgery, cornea laser heating and cancer hyperthermia are used. Spatial and transient heating source, on the skin surface and inside biological body, are considered by using step heating, sinusoidal and constant heating. Mathematical simulations describe a non-Fourier process. Exact solution for the corresponding non-Fourier bioheat transfer model that has time lag in its heat flux is proposed using strongly continuous semigroup theory in conjunction with variational methods. The abstract differential equation, infinitesimal generator and corresponding strongly continuous semigroup are proposed. It is proved that related semigroup is a contraction semigroup and is exponentially stable. Mathematical simulations are done for skin burning and thermal therapy in 10 different models and the related solutions are depicted. Unlike numerical solutions, which suffer from uncertain physical results, proposed analytical solutions do not have unwanted numerical oscillations. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Thermal wave Analytical solution Semigroup theory Thermal therapy Temperature control

1. Introduction In recent years, temperature predictions for living tissues have great attraction due to its significance in basic and clinical sciences. The thermal wave model of bioheat transfer in living tissues reduces to one differential equation (for example see: [1,2]). Although this phenomenon is three-dimensional (3D), for many researchers the depth effects of the heat propagation is a matter of importance. Thus, without loss of generality, one can assume the heat energy propagation along the direction perpendicular to the skin surface. Liu et al. [3] introduced the thermal wave model to investigate physical mechanisms and the behaviours in living tissues. Traditionally, finite difference, boundary element and finite element methods are used [4–8] to model and solve different bioheat transfer problems numerically. However, if both closed form analytical and numerical solutions exist, the analytical one is preferred, since it does not depend on the dimensions of the problem. Moreover, in contrast to numerical solutions the analytical solutions are not depended on the previous grids in the domain. But for most of the numerical solutions of thermal wave model there are oscillations that are hard to dampen [9–11] and these solutions are physically doubtful results for the temperature prediction at the pulse surface heat flux. It is

n

Corresponding author. Tel.: þ 98 2182883445. E-mail addresses: [email protected] (A. Malek), [email protected] (G. Abbasi). http://dx.doi.org/10.1016/j.compbiomed.2015.03.030 0010-4825/& 2015 Elsevier Ltd. All rights reserved.

for these reasons that we aimed in this paper to present several closed form analytical solutions to the thermal wave model of bioheat transfer under different types of boundary conditions. Up to now, Haji-Sheikh et al. described a method of solution of the thermal wave equation in finite bodies using Green's functions [12]. Ahmadikia et al. [13–15] derived the solution of the hyperbolic heat equation using the Laplace transform. This paper deals with the closed analytical form solution of the thermal wave equation based on strongly continuous semigroup (C0-semigroup) theory. Proposed general solutions in this paper work for different kinds of source terms. Therefore, they are flexible in calculating the temperature distribution for various practical problems in skin burning and hyperthermia. Furthermore, the C0-semigroup solution is capable to dealing with the transient or space-dependent boundary conditions. In the recent years, Malek et al. [16,17] have performed several closed forms of analytical solutions for Pennes' and dual phase lag equation using C0-semigroup theory. By semigroup theory, an infinitesimal generator is used to establish the abstract differential equation related to the original thermal wave equation. The main results consist of computing eigenvalues and proof for the well-posedness of related operator. It is proved that an infinitesimal generator is Riesz spectral operator and the corresponding system is exponentially stable. The paper is organized as follows: in Section 2 mathematics formulation for thermal wave equation under generalized boundary conditions are described. In Section 3 semigroup formulation and closed analytical form solution are achieved. The mathematical simulation results based on semigroup theories derived for four

66

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

different boundary conditions in Section 4. Concluding remarks are given in Section 5.

2. Mathematical formulation In order to study thermal behaviour in living tissues, several models describing bioheat transfer have been developed [18]. Most of the thermal medical practices are based on the wellknown Pennes' equation. It is based on an infinitely fast propagation of temperature disturbance (Fourier's law). In contrast with Fourier's constitutive law, a modified flux model based on the finite speed of propagation of heat in the living body was suggested [19,20]. The assumption of finite speed wave propagation does not obey in the classical Fourier's law. The thermal wave theory based on non-Fourier's law can be derived from ! ! qð r ; t þ τq Þ ¼  k∇Tð r ; tÞ ð1Þ where t is the time, k is the thermal conductivity, q is the heat flux, ! ∇T is the temperature gradient, r ¼ ðx; y; zÞ is three-dimensional position vector in which z stands for the tissue depth. The relaxation time is τq ¼ α=C 2 where α is the thermal diffusivity and C stands for the heat propagation velocity. In contrast with common metals the relaxation time for heating processes is much longer than 10  8 s [10]. For the processed meat Mitra et al. [21] took τq ¼ 15 s, while in Refs. [22,23] it is shown that τq is in the range of 20–30 s for biological bodies. Pennes' bioheat transfer equation is [24] ∇q þ wb cb ðT b  TÞ þ Q m þ Q r ¼ ρc

∂T : ∂t

ð2Þ

Here, ρ, c and T denote density, specific heat and temperature of tissue, respectively. cb and wb are the specific heat and perfusion rate of blood respectively. Qm is the metabolic heat generation and Qr is the heat source for spatial heating. Tb is the arterial temperature and was regarded as a constant. Applying the firstorder Taylor's series expansion to Eq. (1) leads to ! ∂qð r ; tÞ ! ! ¼  k∇Tð r ; tÞ: ð3Þ qð r ; tÞ þ τq ∂t By taking the gradient with respect to r on both sides of (3) and substitute in (2) one can write ∇  ðk∇TÞ þ wb cb ðT b  TÞ þ Q m þ Q r
2 
 ∂T ∂Q ∂Q ∂ T ∂T : þ τq  wb cb þ m þ r ¼ ρc τq 2 þ ∂t ∂t ∂t ∂t ∂t

ð4Þ

This hyperbolic heat transfer equation is more complicated than Pennes' equation (2) since for τq ¼ 0, Eq. (4) reduces to Eq. (2). For constant thermal parameters k; Q m the z-direction form of Eq. (4) is
 ∂2 T ∂T ∂Q k 2 þ wb cb ðT b  TÞ þ Q m þ Q r þ τq  wb cb þ r ∂t ∂t ∂z
2  ∂ T ∂T ¼ ρc τ q 2 þ : ð5Þ ∂t ∂t In order to solve this equation we take  ∂Tðz; tÞ ¼ 0; Tðz; 0Þ ¼ T 0 ðzÞ and ∂t t ¼ 0

ð6Þ

in which T 0 ðzÞ for every z A ½0; l is the solution of the following ordinary differential equation [25]: 2

d T 0 ðzÞ þ wb cb ðT b  T 0 ðzÞÞ þ Q m ¼ 0 k 2 dz dT 0 ðzÞ ¼ h0 ðT f  T 0 ðzÞÞ; z ¼ 0; k dz T 0 ðzÞ ¼ T c ; z ¼ l;

ð7Þ

where Tðz; 0Þ ¼ T 0 ðzÞ is the initial heat value at starting, Tf is the ambient temperature, h0 is the coefficient corresponding to the ambient and skin surface temperature and Tc is the body core temperature at z¼l. The closed form analytical solution for Eq. (7) can be written as
  pffiffiffi pffiffiffi pffiffiffi Q h0 A cosh ð AzÞ þ Tc Tb  m sinh ð AzÞ Q wb c b k T 0 ðzÞ ¼ T b þ m þ pffiffiffi pffiffiffi pffiffiffi h0 w b cb A coshð AlÞ þ sinhð AlÞ k
 p ffiffiffi h0 Q T f T b  m sinhð Aðl  zÞÞ k wb cb ; ð8Þ þ pffiffiffi pffiffiffi pffiffiffi h0 sinhð AlÞ A coshð AlÞ þ k where A ¼ wb cb =k. In practice, different types of boundary conditions (BCs) may be considered with Eq. (5) and initial conditions (6): (a) a transient surface heat flux, (b) a surface and body core heating, (c) a cooling medium on tissue surface, (d) a surface heating with body core heat flux, which can be generalized as  ∂Tðz; tÞ k ¼ f 0 ðtÞ; Tðl; tÞ ¼ T c ; ð9aÞ ∂z z ¼ 0

Fig. 1. Schematic of a tissue model: (a) skin burns by the tissue surface heating. (b) Hyperthermia by the point heating power, heat is deposited through inserting a heating probe in the center of tumor site (z ¼ z0 ).

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

Tð0; tÞ ¼ f 1 ðtÞ;  ∂Tðz; tÞ k ∂z 

z¼0

Tð0; tÞ ¼ f 1 ðtÞ;

Tðl; tÞ ¼ T c ;

ð9bÞ

¼ hf ðf 2 ðtÞ  Tð0; tÞÞ;

Tðl; tÞ ¼ T c ;

ð9cÞ

 ∂Tðz; tÞ ¼ 0; ∂z z ¼ l

L2 ð0;lÞ

3. C0-semigroup theory for thermal wave equation Define Tðz; tÞ ¼ uðz; tÞ þ vðz; tÞ where vðz; tÞ satisfies one of the BCs (9a)–(9d). Then uðz; tÞ can be written as
 ∂2 ∂ ∂2 ρcτq 2 þ ðρc þ τq wb cb Þ þ wb cb  k 2 uðz; tÞ ∂t ∂z ∂t
 ∂ Q r ðz; tÞ ¼ 1 þ τq ∂t
 ∂2 ∂ ∂2 þ  ρcτq 2  ðρc þ τq wb cb Þ þk 2  wb cb vðz; tÞ ∂t ∂z ∂t þQ m þwb cb T b ;

ð10Þ

for the following initial conditions:   ∂uðz; tÞ ∂vðz; tÞ uðz; 0Þ ¼ T 0 ðzÞ  vðz; 0Þ; ¼  ∂t  ∂t  t¼0

t¼0

;

ð11Þ

with one of the boundary conditions:  ∂uðz; tÞ ¼ 0; uðl; tÞ ¼ 0; ∂z z ¼ 0 uðl; tÞ ¼ 0;

 ∂uðz; tÞ k ∂z 

¼ hf uð0; tÞ;

z¼0

uð0; tÞ ¼ 0;

ð12aÞ ð12bÞ

uðl; tÞ ¼ 0;

ð12cÞ

 ∂uðz; tÞ ¼ 0: ∂z z ¼ l

ð12dÞ

In order to solve this PDE, we need to study some lemma and theorems. Lemma 1. Consider the operator A0 on L2 ð0; lÞ as A0 h ¼ 

α d2 h w b c b þ h; τq dz2 ρcτq

with domain  dh are absolutely continuous; DðA0 Þ ¼ h A L2 ð0; lÞj h and dz 2

d h dz

2

Proof. For h1 A DðA0 Þ and h2 A L2 ð0; lÞ one may use integration by parts to write * + α d2 h1 wb cb 〈A0 h1 ; h2 〉L2 ð0;lÞ ¼  þ h ; h τq dz2 ρcτq 1 2 !  l Z l dh1 dh2 α d2 h2 wb cb h2  h1 ¼ þ h1  þ h dz: 2 τq dz2 ρcτq dz dz 0 0

ð9dÞ

where f 0 ðtÞ and f 1 ðtÞ are the time-dependent surface heat flux and the temperature at the skin surface respectively. Expressions for f 0 ðtÞ and f 1 ðtÞ will be defined according to the different heating style (see Fig. 1). In BC (9c), f 2 ðtÞ is time-dependent temperature for cooling medium and hf is the heat convection coefficient between the medium and the skin surface. For BCs (9a)–(9c) the body core temperature was regarded as a constant Tc, considering that the biological body tends to keep its core temperature to be stable. In BC (9d) it is assumed that heat flux vanishes for z ¼l. In the next section it is proved that the initial boundary value problem (5) and (6) with one of the BCs (9a)–(9d) leads to strongly continuous one-parameter semigroup of bounded linear operator on the specific Hilbert space.

uð0; tÞ ¼ 0;

67

ð13Þ

)

A L2 ð0; lÞ and h satisfies in one of the BCs ð12aÞ–ð12dÞ ;

This can be written in the form 〈h1 ; A0 h2 〉 if and only if h2 A DðA0 Þ thus A0 is self-adjoint. It can be shown that for all non-zero h A DðA0 Þ, A0 is positive, i.e.,

α dh dh w c ; 〈A0 h; h〉L2 ð0;lÞ ¼ þ b b 〈h; h〉L2 ð0;lÞ 4 0: ð15Þ τq dz dz L2 ð0;lÞ ρcτq From (15), it follows that A0 is coercive, i.e., 〈A0 h; h〉L2 ð0;lÞ Z

and α ¼ k=ρc. Then operator A0 is self-adjoint, positive and boundedly invertible.

for all h A DðA0 Þ;

ð16Þ

this shows that A0 is injective. Surjectivity is implied by ranðA0 Þ ¼ ðker An0 Þ ? ¼ ðker A0 Þ ? , thus A0 is invertible and bounded from the ‘Closed Graph Theorem’. For ξ A ran A0 , from the Cauchy– Schwarz inequality we have ‖A0 1 ξ‖2 r

ρc τ q ρc τ q 〈A0 A0 1 ξ; A0 1 ξ〉L2 ð0;lÞ r J ξ J J A0 1 ξ J : wb cb wb c b ρc τ

Thus J A0 1 ξ J r wb cqb J ξ J and A0 is boundedly invertible.

□

1 2

From positivity1 property of A0 the square root A0 exists. Thus, the space H ¼ DðA20 Þ  L2 ð0; lÞ with the inner product 1

1

〈u; w〉H ¼ 〈A20 u1 ; A20 w1 〉L2 ð0;lÞ þ 〈u2 ; w2 〉L2 ð0;lÞ ;

ð17Þ

T Recall where u ¼ ðu1 u2 ÞT and w1 ¼ ðw 1 w2 Þ defines a Hilbert space. 1 1 that if k A DðA0 Þ, then 〈A20 k; A20 h〉 ¼ 〈A0 k; h〉 for every h A DðA20 Þ. In the following theorem, we show how to rewrite thermal wave model as an abstract differential equation. By abstract differential equation we forget the spatial dependence, and we write the PDE as (abstract) a system of ordinary differential equations [26].

Theorem 1. The problem (10) and (11) with one of the BCs (12a)– (12d) can be formulated as an inhomogeneous abstract differential equation ! ! 8 u u > >d > ¼A þ BQ r ðtÞ þ DvðtÞ þ f > > ut > < dt ut 0 1 ! ð18Þ T 0  vð0Þ > u  > > @ ∂v A > ¼  > >  ð0Þ : ut t ¼ 0 ∂t on the Hilbert space H with the inner product 〈u; w〉H . Here u1 t ¼ du , dt f ¼ ð0 Q m þ wb cb T b ÞT and A with domain DðAÞ ¼ DðA0 Þ  DðA20 Þ, is a closed and densely defined operator on H as ! ! ! 0 I u1 u1 ¼ ; ð19Þ A  A0  βI u2 u2 τ w c

where I is identity operator and β ¼ τ1q ð q ρcb b þ 1Þ and ! 0  Qr; BQ r ¼ 1 ∂ τq ρc 1 þ τ q ∂t Dv ¼

ð14Þ

wb cb ‖h‖2L2 ð0;lÞ ρcτq

0  ∂t∂ 2  β ∂t∂  A0 2

! v:

ð20Þ

Proof. It is similar to the proof of wave equation presented in Ref. [27]. □

68

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

Table 1 The sequence of functions ϕn ðzÞ that satisfy the boundary conditions (12a)–(12d). Boundary condition

pn

ϕn ðzÞ

(12a) (12b) (12c) (12d)

ðn 1=2Þπ=l nπ=l Roots of pn cotðpn lÞ ¼  hf =k ðn 1=2Þπ=l

cos ðpn zÞ sin ðpn zÞ sin ðpn ðl  zÞÞ sin ðpn zÞ

Lemma 2. Let A be the operator on Hilbert space H introduced in Theorem 1. The adjoint of A on H with DðAn Þ ¼ DðAÞ is given by ! ! ! 0 I w1 w1 An ¼ ð21Þ A0  βI w2 w2

Lemma 3. The eigenfunctions fψ 7 n ; n ¼ 1; 2; …g forms an Riesz basis on Hilbert space H. is an orthogonal Proof. The sequence of functions fψ 7 n ðzÞg1 n¼1 sequence with inner product 〈 .,. 〉H . i.e., 8 if m an > 7n :2 τq n ρcτq 〈ψ n ; ψ  n 〉H ¼

□

pffiffiffiffiffiffiffiffi  1, define the linear operator L and XN as ! ! N X ϕn ðzÞ ϕn ðzÞ ψ ; X ¼ α ¼ n nπ i N 7n ϕn ðzÞ 7 nπl iϕn ðzÞ l ¼ N

For i ¼ L

Proof. The lemma is easily proved by using Lemma 1.

n

In the next theorem we show that A generates a contraction semigroup ðSðtÞÞt Z 0 on H; i.e., J SðtÞ J r 1 for every t Z 0. Theorem 2. The operator A in (19) is the infinitesimal generator of a contraction C0-semigroup ðSðtÞÞt Z 0 on H.

n

N X

¼

1=2

ð22Þ

then from (19) we can write u2 ¼ λu1 and  A0 u1  β u2 ¼ λu2 thus
 α ∂2 wb cb þ þ λ ð λ þ β Þ u1 ¼ 0 for u1 A DðA0 Þ: ð25Þ  τq ∂z2 τq ρc 1

Let u1 A fϕn ðzÞgn ¼ 1 be the sequence of functions on DðA0 Þ defined in Table 1. Then Eq. (25) is equivalent to
 1 w c ð26Þ λ2n þ βλn þ αp2n þ b b ¼ 0; n ¼ 1; 2; …: τq ρc In general, for every n this equation has two roots: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
rffiffiffiffiffi2
 β β 1 2 α  λ7n ¼  7  pn ; n ¼ 1; 2; …: 2 2 τq τq

ð27Þ

It is easy to see that

ψ 7 n ðzÞ ¼ ðϕn ðzÞ λ 7 n ϕn ðzÞÞT ;

N X

þ

Similarly, it can be shown that

In the following we find solution of the inhomogeneous abstract differential equation (18). We derive solution of (18) by showing that the normalized eigenfunctions of A form a Riesz basis in H. We begin by calculating the eigenvalues and eigenfunctions of A. Since ! ! u1 u1 A ¼λ ; ð24Þ u2 u2

n ¼ 1; 2; …

lies in the domain of A and satisfies Aψ 7 n ¼ λ 7 n ψ 7 n . Hence ψ 7 n and λ 7 n are corresponding eigenfunctions and eigenvalues of A. 1 The simple eigenvalues fλ 7 n gn ¼ 1 are unbounded and located on the vertical line  β =2 except for a finite number real eigenvalues. 1 One can easily see that the closure of fλ 7 n gn ¼ 1 is totally disconnected.

L2

na0

L2

jαn þ α  n j2 〈ϕn ; A0 ϕn 〉L2

n¼1

Thus by Corollary 2.2.3 in Ref. [28], A is the infinitesimal generator of a C0-semigroup ðSðtÞÞt Z 0 on H. From (22) and (23), it follows that A and An are dissipative on H. Thus ðSðtÞÞt Z 0 is a contraction C0-semigroup from Lumer–Phillips Theorem [29]. □

n

na0

n¼1

¼ Reð〈A0 w2 ; w1 〉L2 ð0;lÞ þ 〈  A0 w1  βw2 ; w2 〉L2 ð0;lÞ Þ ¼  β‖w2 ‖2L2 ð0;lÞ :

ð23Þ

αn

2

2

X

N

1=2 X

N

2

‖LðX N Þ‖H ¼ A0 α n ϕn þ αn λn ϕn



¼ N

¼ N

Reð〈Aw; w〉H Þ ¼ Reð〈A0 w2 ; A0 w1 〉L2 ð0;lÞ þ 〈 A0 w1  βw2 ; w2 〉L2 ð0;lÞ Þ

Reð〈w; An w〉H Þ ¼  β ‖w2 ‖2L2 ð0;lÞ :

na0

for arbitrary N A N and arbitrary scalars

Proof. For w A DðAÞ we show that 1=2


 l α w c λn λ  n þ p2n þ b b : 2 ρcτq τq

jαn λn þ α  n λ  n j2 ‖ϕn ‖2L2

ð28Þ

and ‖X N ‖2H

2

2

X

N

1=2 X

N nπ i

¼ A0 α n ϕn þ αn ϕn

l

¼ N

¼ N n

¼

N X

na0

n

L2

na0

L2

jαn þ α  n j2 〈ϕn ; A0 ϕn 〉L2 þ

n¼1

N X

jαn  α  n j2

n¼1

nπ 2 l

‖ϕn ‖2L2 : ð29Þ

Now, from (28) and (29), for sufficiently large values of n, one can write jαn λn þ α  n λ  n j2 nπ 2 -C a 0; jαn  α  n j2 l

ð30Þ

since the numerator and the denominator are of order n2. Thus, for every N there exist positive constants m and M, such that m J X N J r J LðX N Þ J rM J X N J . i.e., m r J L J r M. Thus, L is bounded and L  1 exists and it is bounded. Moreover, L is surjective, since the domain H is a Banach space. Therefore fψ 7 n ðzÞg1 is a Riesz n¼1 basis. Note that a Riesz basis is the image of an orthogonal basis under an invertible linear bounded operator [28]. □ From Lemma 3 one can conclude that operator A in (19) is a Riesz1 spectral operator with simple eigenvalues fλ 7 j gj ¼ 1 and corresponding eigenfunctions fψ 7 j ðzÞg1 . A representation of operator A and its j¼1 corresponding semigroup are established in the next theorem. Theorem 3 (Curtain and Zwart [28]). The Riesz-spectral operator A in (19) satisfies (a) A has the representation * + 1 X ψ ψj λj h; j Ah ¼ Jψj J Jψj J ¼ 1 j

ja0

for h A DðAÞ:

ð31Þ

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

(b) A is the infinitesimal generator of a C0-semigroup ðSðtÞÞt Z 0 if and only if supj Z 1 Reðλ 7 j Þ o 1 and ðSðtÞÞt Z 0 is given by * + 1 X ψj ψj ð32Þ eλj t ξ; SðtÞξ ¼ J ψ J J ψj J ¼ 1 j j

ja0

(c) The growth
bound of the  semigroup is given by 1 w0 ¼ inf log J SðtÞ J ¼ supReðλ 7 j Þ t40 t jZ1

ð33Þ

69

Applying S(t) on the both sides yields 0 1 0 1 T 0  vð0Þ T 0  vð0Þ A þ SðtÞ@ dv A  @ dv  ð0Þ  ð0Þ dt dt ! ! Z t uðsÞ ¼ A þ BQ r ðsÞ þ DvðsÞ þ f ds ut ðsÞ 0 Z t Sð  sÞðBQ r ðsÞ þ DvðsÞ þ f Þ ds: SðtÞ 0

Proof. See [28], Theorem 2.3.5.

□

From Theorem 3 solution ðu ut ÞT for inhomogeneous abstract differential equation (18) can be written as ! ! Z t T 0 ðzÞ  vðz; 0Þ u Sðt  sÞðBQ r ðz; sÞ þ Dvðz; sÞ þ f Þ ds: ¼ SðtÞ þ  ∂vðz;0Þ ut ∂t 0 ð34Þ

From (35) one can write 0 1 ! T 0  vð0Þ Z t uðtÞ Aþ Sðt  sÞðBQ r ðsÞ þ DvðsÞ þ f Þ ds; ¼ SðtÞ@ dv ut ðtÞ  ð0Þ 0 dt which proves the assertion.

□

Our next goal is to show ðSðtÞÞt Z 0 in (32) is exponentially stable, i.e., the solution for system (18) tends to zero exponentially fast as t-1.

T Let Q r ðz ,. Þ; vðz ,. Þ A C 1 ð0; t f Þ and ðT 0 ðzÞ  vðz; 0Þ  ∂vðz;0Þ ∂t Þ A D 1=2 ðA0 Þ  L2 ð0; lÞ, we call (34) a strongly solution for (18). If Q r ðz ,. Þ; vðz ,. Þ A H 1 ð0; t f Þ (H1 is the Sobolev space) and ðT 0 ðzÞ  1 T 2 vðz; 0Þ  ∂vðz;0Þ Þ A DðA Þ 0  L 2 ð0; lÞ, we call (34) a mild solution for ∂t (18).

Definition 1. The C0-semigroup ðSðtÞÞt Z 0 on the Hilbert space H is exponentially stable if there exist positive constants M and α such that

Theorem 4. Any arbitrary mild solution of abstract differential equation (18) can be written uniquely in the form (34).

Theorem 5. The C0-semigroup ðSðtÞÞt Z 0 in (32) is exponentially stable.

Proof. Any arbitrary ðu ut ÞT as a mild solution of (18) must satisfy the integrated form of (18) ! ! 0 T 0  vð0Þ 1 Z ! t uðtÞ uðsÞ @ A dv þ A ¼ þ BQ r ðsÞ þ DvðsÞ þ f ds: ut ðtÞ ut ðsÞ  ð0Þ 0 dt

Proof. The coefficients of the quadratic equation (26) are positive, thus all solutions of (27) have negative real parts. The eigenvalues do not converge to zero since w0 ¼ λ1 of (33) is strictly smaller than zero. This proves the theorem. □

ð35Þ Applying ASð  tÞ on the both sides and integrating yields 0 1 T 0  vð0Þ Z t @ Ads dv ASð sÞ  ð0Þ 0 dt ! Z t uðsÞ ¼ ASð  sÞ ds ut ðsÞ 0 ! ! ! Z s Z t uðτÞ ASð sÞ A þ BQ r ðτÞ þ DvðτÞ þ f dτ ds  ut ðτÞ 0 0 ! ! Z t Z t uðsÞ uðsÞ ASð  sÞ A ds þ Sð  tÞ ¼ ut ðsÞ ut ðsÞ 0 0 þ BQ r ðsÞ þ DvðsÞ þ f Þds ! ! Z t uðsÞ Sð  sÞ A þ BQ r ðsÞ þ DvðsÞ þ f ds  ut ðsÞ 0 ! ! Z t uðsÞ A ¼ Sð  tÞ þ BQ r ðsÞ þ DvðsÞ þ f ds ut ðsÞ 0 Z t  Sð  tÞðBQ r ðsÞ þ DvðsÞ þ f Þds 0

and T 0 vð0Þ

Sð  tÞ Z ¼ Sð  tÞ t

 0

t

A 0

Z

ð0Þ  dv dt

! þ

uðsÞ ut ðsÞ

!

T 0  vð0Þ ð0Þ  dv dt

! !

þBQ r ðsÞ þDvðsÞ þ f ds

Sð  sÞðBQ r ðsÞ þ DvðsÞ þ f Þ ds

J SðtÞ J r Mexpð  αtÞ

for t Z 0:

The analytic solution (34) involves imaginary coefficient

λj ; j ¼ 71; 72; … and it is not useful to depict the graphs numerically. In the next theorem separation of variables is used to give real expression for the complex expression (32). Theorem 6. A series expression for the semigroup ðSðtÞÞt Z 0 described by (32) is equivalent to 0 1 1 X γ j t þ d e  γ j t Þϕ ðc e j j j B C ! B C u0 j¼1 B

  C SðtÞ ¼ e  ðβ=2Þt B X C 1
u1 B C β β @  γ j eγ j t  dj  γ j e  γ j t ϕj A  cj 2 2 j¼1 ð36aÞ pffiffiffiffiffiffiffiffi for the case j 2  τq β j Z 2pj ατq while ! u0 SðtÞ ¼ e  ðβ=2Þt u1 0 1 1 X ðaj cos ðωj tÞ þ bj sin ðωj tÞÞϕj B C B C j¼1 B 


 C B X C B 1 C β β @ bj þ aj ωj sin ðωj tÞ ϕj A bj ωj  aj cos ðωj tÞ  2 2 j¼1 ð36bÞ pffiffiffiffiffiffiffiffi for the case j 2  τq β j o 2pj ατq where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
rffiffiffiffiffi2


rffiffiffiffiffi2  β 1 2 α α β 1 2   γj ¼  pj pj  ; ωj ¼ ; 2 τq 2 τq τq τq and aj ¼ 〈u0 ; ϕj 〉L2 ð0;lÞ ; bj ¼

1

ωj

β

ð〈u1 ; ϕj 〉L2 ð0;lÞ þ 〈u0 ; ϕj 〉L2 ð0;lÞ Þ; 2

70

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

cj ¼

1 2γ j





β 2

1 dj ¼  2γ j

  þ γ j 〈u0 ; ϕj 〉L2 ð0;lÞ þ 〈u1 ; ϕj 〉L2 ð0;lÞ ;



β 2



   γ j 〈u0 ; ϕj 〉L2 ð0;lÞ þ 〈u1 ; ϕj 〉L2 ð0;lÞ :

ð37Þ

Proof. Substituting uðz; tÞ ¼ MðzÞNðtÞ in the homogeneous form of (10) and rearrangement leads to the equation

τq ρcN″ ðtÞ þ ðτq wb cb ÞN 0 ðtÞ NðtÞ

¼

″

kM ðzÞ  wb cb MðzÞ : MðzÞ

ð38Þ

Since the left-hand side depends only on t and the right-hand side depends only on z. It follows that each side is equal to a constant 1 2 ðwb cb þ kpj Þ that is derived from taking MðzÞ A fϕj ðzÞgj ¼ 1 in Table 1 corresponding to the suitable boundary condition, we have the following equation:

τq ρcN″j ðtÞ þ ðτq wb cb ÞN 0j ðtÞ þ ðwb cb þ kp2j ÞN j ðtÞ ¼ 0:

ð39Þ

Solution for Eq. (39), for each value of n, is N j ðtÞ ¼ e  ðβ=2Þt ( aj cos ðωj tÞ þ bj sin ðωj tÞ  cj eγ j t þ dj e  γ j t

pffiffiffiffiffiffiffiffi if j 2  τq β j o 2pj ατq otherwise

ð40Þ

By applying the initial conditions (11), aj ; bj ; cj and dj can be calculated. Thus the explicit solution of the homogeneous problem (10) and (11) is uðz; tÞ ¼ e

 ðβ=2Þt



∂ ð1 þ τq ÞQ r ; ϕj ∂t τ q ρc L2 ð0;lÞ !
2 

∂ ∂ þ β þ A0 v; ϕj þ〈Q m þ wb cb T b ; ϕj 〉L2 ð0;lÞ ds ϕj  ∂t ∂t 2 L2 ð0;lÞ 

8 1 X > > > ðaj cos ðωj tÞ þ bj sin ðωj tÞÞϕj ðzÞ > > > > ðcj eγ j t þ dj e  γ j t Þϕj ðzÞ > > :j¼1

otherwise

pffiffiffiffiffiffiffiffi

ατq

ð41Þ

and explicit expression of C0-semigroup (32) is given as ! ! u0 u SðtÞ ¼ : □ u1 ut

ð43bÞ pffiffiffiffiffiffiffiffi for the case j 2  τq β j o 2pj ατq . In some practical problems, the function vðz; tÞ does not have the second derivative with respect to time. Thus the solution (43a) and (43b) does not work. This happens when one faces the pulse heat flux at the skin surface with BC (9a) (see Section 4.3.3). In this case, variational boundary value problem (VBVP) corresponding to the boundary value problem (5) with the BC (9a) is used successfully. Note that solutions based on VBVP demands less smoothness [30]. By using the transformation T 1 ðz; tÞ ¼ Tðz; tÞ  T c , Eq. (5) can be changed to the following initial boundary value problem:
 ∂2 T 1 ∂T 1 ∂Q r k 2 þ wb cb ððT b  T c Þ  T 1 Þ þ Q m þ Q r þ τq  wb cb þ ∂t ∂t ∂z
2  ∂ T 1 ∂T 1 ð44Þ ¼ ρc τ q 2 þ ∂t ∂t T 1 ðz; 0Þ ¼ T 0 ðzÞ  T c

j¼1

γj

ð42Þ

〈T 0 ðzÞ  vðz; 0Þ; ϕj 〉L2 ð0;lÞ

!


 β ∂vðz; 0Þ ; ϕj  sinhðγ j tÞ þ γ j coshðγ j tÞ  sinhðγ j tÞ ϕj ∂t 2 L2 ð0;lÞ Z 1 X 1



t

1




∂





e  ðβ=2Þðt  sÞ sinhðγ j ðt  sÞÞ 1 þ τq Q r ; ϕj γ τq ρc ∂t L2 ð0;lÞ j¼1 j 0 !
2 

∂ ∂ þ A v;  þ β ϕ þ 〈Q þ w c T ; ϕ 〉 dsϕj 0 m b b b L ð0;lÞ j j 2 ∂t ∂t 2 L2 ð0;lÞ

þ

ð43aÞ pffiffiffiffiffiffiffiffi

for the case j 2  τq β j Z 2pj ατq while uðz; tÞ ¼ e  ðβ=2Þt


β

1 X 1 j¼1

ωj

〈T 0 ðzÞ  vðz; 0Þ; ϕj 〉L2 ð0;lÞ

sin ðωj tÞ þ ωj cos ðωj tÞ



2 !

∂vðz; 0Þ ; ϕj sin ðωj tÞ ϕj  ∂t L2 ð0;lÞ Z t 1 X 1 þ e  ðβ=2Þðt  sÞ sin ðωj ðt  sÞÞ j¼1

ωj

0

 ∂T 1 ðz; tÞ ¼0 ∂t t ¼ 0

ð45Þ

ð46Þ

Theorem 7. A VBVP corresponding to the problem (44) with the BC (46) is in the following form: Find a function T 1 ðz; tÞ that belongs to Hilbert space L2 ð0; lÞ and that satisfies the equation aðT 1 ; ϕÞ ¼ 〈φ; ϕ〉L2 ð0;lÞ

From (34), the explicit solution for inhomogeneous problem (10) and (11) yields 1 X 1

and

with the boundary conditions  ∂T 1 ðz; tÞ ¼ f 0 ðtÞ; T 1 ðl; tÞ ¼ 0 k ∂z z ¼ 0

where aðT 1 ; ϕÞ ¼ k

uðz; tÞ ¼ e  ðβ=2Þt

1



ð47Þ

∂T 1 dϕ ; ∂z dz L2 ð0;lÞ



∂T 1 ∂2 T 1 þ ρc τ q 2 ; ϕ þ wb cb T 1 þ ðτq wb cb þ ρcÞ ∂t ∂z L2 ð0;lÞ

∂Q 〈φ; ϕ〉L2 ð0;lÞ ¼ wb cb ðT b  T c Þ þ ðQ m þQ r Þ þ τq r ; ϕ ∂t L2 ð0;lÞ  l ∂T 1 þk ϕ ∂z 0 for all functions ϕ such that ϕ A H 1 ð0; lÞ and satisfies in the BC (12a), a : H 1  H 1 -R is a bilinear form and φ : H 1 -R is a linear functional. Proof. By multiplying ϕ and integrating by parts from Eq. (44) the theorem is proved. □ P Let T 1 ðz; tÞ ¼ n T 1n ðtÞϕn ðzÞ where functions ϕn are selected from Table 1. From Eq. (47) for ϕ ¼ ϕn ; n ¼ 1; 2; …, we have

 2 d T 1n ðtÞ α 2 wb cb 1 τ q wb cb dT 1n ðtÞ þ 1 T ¼  p þ ðtÞ  1n n 2 dt ρc τ q τq ρc τq dt 1 þ 〈φ; ϕn 〉L2 ð0;lÞ ð48aÞ ρc τ q T 1n ð0Þ ¼ 〈ðT 0 ðzÞ  T c Þ; ϕn 〉L2 ð0;lÞ ;

dT 1n ðtÞ ¼0 dt

ð48bÞ

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

The second-order ordinary differential equation (48a) and (48b) can be solved using C0-semigroup theory according to what is said in this section. In Section 4, several mathematical simulations are considered and solved by C0-semigroup theories to show the practicality and efficiency of this technique in various medical thermotherapy and skin burning.

4. Results and discussions In this section, for the finite tissue domain, the presented solutions will incorporate relatively to situations such as volumetric heating on the surface or inside the biological body for the transient or space-dependent boundary conditions. The solution based on whatever is derived in Sections 2 and 3 is straightforward and easy to utilize for thermal wave models that may encounter in tissue property measurement, cancer hyperthermia, laser surgery and evaluation of the thermal injury. In practice, one computes the temperature responses subject to the constant, sinusoidal, pulse and point heating in the volume or on boundary of the tissue. Numerical simulations are done using the examples appeared in Refs. [10,11,25] with the following tissue properties [3,25]: ρ ¼ 1 1000 kg m  3 , c ¼ cb ¼ 4200 J kg =1C, T b ¼ T c ¼ 37 1C, l ¼ 0:03 m 1 k ¼ 0:5 Wm =1C, wb ¼ 0:5Wm  3 =s, Q m ¼ 3380 Wm  3 , h0 ¼ 10 Wm  2 =1C, T f ¼ 25 1C, hf ¼ 100 Wm  2 =1C and τq ¼ 20 s. In cancer hyperthermia, skin burning is implemented using microwave, radio frequency, and electromagnetic apparatus. In these cases one must include the heating source Qr in the formulation. Although many different types of heating source Qr may be used, in the following three most practical ones are studied.

4.1. Heat transfer in tissue subject to exponential heating source When one deals with microwave, laser or ultrasound the heating source constructed from the Beer's law [31] the heat flux qr ðz; tÞ ¼ P 0 ðtÞexpð  ηzÞ decays exponentially with the distance from the skin surface. In these cases, the heating source is expressed by [32] Q r ðz; tÞ ¼ 

∂qr ¼ ηP 0 ðtÞexpð  ηzÞ; ∂t

ð49Þ

4.1.1. Surface adiabatic condition with BC (9a) In laser–tissue interaction [32], for η ¼ 200 m  1 , Fig. 2 depicts the temperature distributions of a biological body where the surface is insulated (see BC (9a), f 0 ðtÞ ¼ 0) subject to two different spatial heatings P 0 ðtÞ ¼ 250 Wm  2 and P 0 ðtÞ ¼ 250 þ 200 cos ð0:02tÞ Wm  2 . Fig. 2(a) and (b) depicts the case of constant and sinusoidal spatial heating respectively. In Fig. 2(b) there is an intercross for temperature as time increases in sinusoidal heating. In Fig. 3 effects of the scattering coefficient η ¼ 20; 50; 100 and 200 m  1 to the temperature response at skin surface are given. The higher temperature is expected for the larger coefficient. In Fig. 3(b) it is shown that for sinusoidal spatial heating, one expects temperature oscillation too. 4.1.2. Pulse heat flux at skin surface with BC (9a) This case may happen when one faces the eye surgery for a short period of time [11]. The surface pulse heat flux ( 1000 Wm  2 t o 1200 s f 0 ðtÞ ¼ ð50Þ 0: t Z 1200 s is applied on the boundary of the tissue. In the presence of laser heating source Q r ðz; tÞ as it is given by Eq. (49), constant spatial heating P 0 ðtÞ ¼ 250 Wm  2 , η ¼ 200 m  1 is used in Fig. 4. In Fig. 4 it is shown that at t ¼1200 s when the surface pulse was stopped the highest temperature happens. This is a time that for small blood perfusion wb ¼0.5, the tissue temperature will reach steady state gradually (Fig. 4(a)). This phenomenon is quickly for wb ¼4 (Fig. 4(b)). Moreover, in this problem the highest temperature occurs at the skin surface. For t Z 1200 the temperature decreases and as we get closer to the body core, temperature value is less. For time close to beginning, as it is shown in Fig. 4 in the higher depth one expects more delay time for a heat flux to respond to the temperature gradient (non-Fourier law). Fig. 5 displays that temperature oscillation due to sinusoidal spatial heating P 0 ðtÞ ¼ 250 þ200 cos ð0:02tÞ Wm  2 for η ¼ 200 m  1 . These simulation results help in the thermal injury analysis, to better understand how to apply various pulse heating. 4.1.3. Sinusoidal surface heat flux; BC (9a) In estimating the blood perfusion, one may use the repeated irradiation from regulated laser; i.e., sinusoidal surface heat flux. Fig. 6 depicts the surface temperature transient when tissue was simultaneously subjected to sinusoidal surface heat flux f 0 ðtÞ ¼ 1000 þ 500 cos ð0:02tÞ Wm  2 and sinusoidal spatial laser heating P 0 ðtÞ ¼ 250 þ 200 cos ð0:02tÞ Wm  2 , η ¼ 200 m  1 for different depths z¼0, 0.012 and 0.024 m. In Fig. 6, oscillation at skin

41

41

40

40

39

39 TemperatureT °C

TemperatureT °C

where P 0 ðtÞ and η are the spatial heating on skin surface and the scattering coefficient respectively.

38 37

t 0 t 500s

36

71

38 37

t 0 t 500s

36

t 1000s 35 34 0.000

t 1000s 35

t 2000s t 2400s 0.005

0.010

0.015

0.020

0.025

0.030

34 0.000

t 2000s t 2400s 0.005

0.010

zm

Fig. 2. Temperature distribution at (b) P 0 ðtÞ ¼ 250 þ 200 cos ð0:02tÞ Wm  2 .

different

0.015

0.020

0.025

0.030

zm

times

for

η ¼ 200 m  1

and

surface

adiabatic

condition

f 0 ðtÞ ¼ 0

subject

to

(a)

P 0 ðtÞ ¼ 250 Wm  2 ;

72

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

41

200

40

100

40

37

200

36

TemperatureT °C

TemperatureT °C

38

100 20 0

2000

4000

6000

8000

20

38

37

36

50

35 34

50

39

39

35

10000

0

500

1000

t s

1500

2000

2500

t s

Fig. 3. Effect of scattering coefficient η ¼ 20; 50; 100 and 200 m  1 on temperature response at skin surface when f 0 ðtÞ ¼ 0 subject to (a) P 0 ðtÞ ¼ 250 Wm  2 ; P 0 ðtÞ ¼ 250 þ 200 cos ð0:02tÞ Wm  2 .

z 0

60

48

z 0.012m

z 0.024m TemperatureT °C

TemperatureT °C

50

45

z 0.012m

46

z 0.024m

55

z 0

40

44 42 40 38

35

36 0

500

1000

1500

2000

2500

0

500

1000

t s

1500

2000

2500

t s

Fig. 4. Transient temperature at three position when pulse heat flux f 0 ðtÞ is applied (subject to constant spatial heating P 0 ðtÞ) for different blood perfusion: (a) wb ¼ 0:5 Wm  3 =s; (b) wb ¼ 4 Wm  3 =s.

z 0

60

z 0

48

z 0.012m

z 0.012m

TemperatureT °C

TemperatureT °C

46

z 0.024m

55

50

45

40

z 0.024m

44 42 40 38

35

36 0

500

1000

1500

2000

2500

0

500

t s

1000

1500

2000

2500

t s

Fig. 5. Transient temperature at three position when pulse heat flux f 0 ðtÞ is applied (subject to sinusoidal spatial heating P 0 ðtÞ) for different blood perfusion: (a) wb ¼ 0:5 Wm  3 =s; (b) wb ¼ 4 Wm  3 =s.

surface is shown while the relaxation time is evident in the depth z¼ 0, 0.012 and z ¼0.024 meter.

by a heating probe such that the cooling medium is applied on the surface (see Fig. 1(b) and BC (9a)). The point heating source is formulated by

4.2. Heat transfer in tissue subject to point heating source

Q r ðz; tÞ ¼ P 1 ðtÞδðz  z0 Þ;

In clinics, the point heating is applied by power heating apparatus where heat in the deep tumor site z ¼ z0 is deposited

where δðz  z0 Þ is the Dirac delta function and P 1 ðtÞ is the point heating power.

ð51Þ

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

4.2.1. Temperature response subject to point heating source under BC (9c) In cancer hyperthermia or thermal comfort analysis one may use BC (9c). Fig. 7(a) illustrates the steady state temperature distribution using (51) with point heating power P 1 ðtÞ ¼ 2500 Wm  3 under cooling temperatures f 2 ðtÞ ¼ 0; 5; 15 and 25 1C for z0 ¼ 0:021 m and t¼500 s. In Fig. 7(b) the influence of heating power to the temperature distribution when tissue is subjected to point heating source Q r ðz; tÞ ¼ P 1 ðtÞδðz  0:021Þ and the constant temperature of cooling medium is f 2 ðtÞ ¼ 15 1C for the t¼ 500 s is shown. 4.3. Surface heating via boundary with no heating source Qr To analyse tissue burns caused by the surface heating via the boundary, the heating source is vanished since heat is incident on the surface tissue.

4.3.1. BC (9a): constant surface heat flux In thermal injury analysis, when one uses a hot plate, the constant heating at the skin surface is encountered [31]. In practice, for the temperature increase one needs the higher surface heating. This is valuable for thermal comfort evaluation. In Fig. 8(a), for constant surface heatings f 0 ðtÞ ¼ 1000 Wm  2 , f 0 ðtÞ ¼ 500 Wm  2 and f 0 ðtÞ ¼ 200 Wm  2 the temperature response at depth z¼0.01208 m are plotted. As it is shown in the beginning, there is a temperature increasing delay and after a while, temperature will increase. 70

z 0

65

z 0.012m z 0.024m

55 50 45 40 35 0

500

1000

1500

2000

2500

t s

Fig. 6. Effect of sinusoidal surface heat flux f 0 ðtÞ and sinusoidal spatial heating P 0 ðtÞ on temperature response at three depths z ¼ 0; 0:012 and 0.024 m.

4.3.2. BC (9a): sinusoidal surface heat flux In this case, the effects of the heating frequency w and blood perfusion wb to the temperature response is considered at the depth z¼0.00208 m for f 0 ðtÞ ¼ 1000 þ 500 cos ðwtÞ Wm  2 . In Fig. 8(b) it is shown that for larger blood perfusion, the shorter time is needed to come into situation that constant oscillation magnitude for temperature happens. 4.3.3. BC (9a): pulse surface heat flux Problem 1. Here, the heating source Qr ¼0 and heat flux (50) are applied for BC (9a). Fig. 9(a) and (b) shows the temperature for different depths for wb ¼0.5 and wb ¼4. In Fig. 9(b) for large blood perfusion, it is shown that the lowest temperature is at the skin surface when the tissue temperature reaches to steady state. This is a different phenomenon with the case when spatial heating exists (see Section 4.3). Problem 2. In this problem Qr ¼ 0, l ¼0.01208 m and heat flux is considered as ( 40000 Wm  2 t o 20 s f 0 ðtÞ ¼ ð52Þ 0: t Z 20 s Fig. 10(a) plots the elevated temperature distributions at different times. Fig. 10(b) shows the effects of the relaxation time τq on the elevated temperature at z ¼0.00208 m. Although Honner [9] stated that it is hard to dampen the numerical oscillations of this kind of problems caused by the pulsed surface heat flux, Liu in year 2008 presented a modified discretization scheme and claimed no numerical oscillations [10]. However, in Liu's Fig. 9, some oscillations and negative values can be observed. These kinds of numerical oscillations that are physically doubtful result do not happen by the proposed method in this paper. This is because we use the analytical close form solution based on C0-semigroups theory. Fig. 10(a) and (b) implies that the proposed solution is efficient and accurate for Problem 2. 4.3.4. BC (9d): constant surface temperature heating The 3D plot of elevated temperature Tðz; tÞ  T 0 ðzÞ profile is shown in Fig. 11. It is shown that under constant heating 12 1C, the tissue temperature can be very different at tissue depth z ¼ 0:00208 m from the tissue surface z ¼0. At the initial stages there is a lag in heating and then takes a jump which caused by a wave front resulted from change in temperature from T 0 ð0Þ to T 0 ð0Þ þ 12 1C at the skin surface. Elevated temperatures predicted by the thermal wave model tend to the steady state after a long time and corresponds to the depth of the point.

P1 0

70

P1 500Wm

50

TemperatureT °C

60 40

30

f2 0°C

TemperatureT °C

TemperatureT °C

60

73

50

3

P1 1000Wm

3

P1 2500Wm

3

P1 5000Wm

3

40

f2 5°C f2 15°C

20

30

f2 25°C 0.000

0.005

0.010

0.015 z m

0.020

0.025

0.030

0.000

0.005

0.010

0.015

0.020

0.025

0.030

z m

Fig. 7. Temperature distribution in tissue as (a) cooling medium temperature f 2 ðtÞ on tissue changed and (b) point-heating power P 1 ðtÞ on tissue changed.

74

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

f0 1000 Wm

48

2

f0 200 Wm

2

60

55 Temperature T °C

Temperature T °C

46

f0 500 Wm

2

44 42

wb 0.5 & w 0.02 wb 0.5 & w 0.01

50

wb 4 & w 0.02 wb 4 & w 0.01

45

40 40

38

0

1000

2000

3000

4000

5000

0

500

1000

t s

1500

2000

2500

t s

Fig. 8. (a) Effect of constant surface heat flux f 0 ðtÞ on the temperature response at depth z¼ 0.01208 m; (b) effect of heating frequency and blood perfusion on the temperature response at the depth z¼ 0.00208 m subject to sinusoidal surface heat flux f 0 ðtÞ.

z 0

z 0 55

46

z 0.012m

50

45

z 0.024m

44 Temperature T °C

Temperature T °C

z 0.024m

z 0.012m

40

42

40

38

36

35 0

500

1000

1500

2000

2500

0

500

1000

t s

1500

2000

2500

t s

Fig. 9. Transient temperature at three position when pulse heat flux is applied for different blood perfusion: (a)wb ¼ 0:5 Wm  3 =s; (b) wb ¼ 4 Wm  3 =s.

t 100 s t 250 s

25

°C

50

T0 z

t 180 s

30

elevated Temperature T z,t

elevated Temperature T z,t

T0 z

°C

35

20 15 10 5

40

30

tq 5 s

20

t q 10 s

10

t q 20 s t q 30 s

0

0 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0

zm

50

100

150

t s

Fig. 10. (a) Elevated temperature Tðz; tÞ  T 0 ðzÞ for τq ¼ 20 s at various times (b) effects of the relaxation time τq on the elevated temperature at z¼ 0.00208 m.

4.3.5. BC (9b): sinusoidal surface heating In medical clinics, the repeated laser irradiation may cause this kind of situation. In Fig. 12 it is shown that near the surface z ¼ 8  10  5 m the temperature profile is similar to the sinusoidal surface heating Tð0; tÞ ¼ T 0 ð0Þ þ5 sin ð0:2π tÞ1C. However, there is delay in temperature profile at the depth z¼0.00208 m. Note that the sinusoidal surface heating may be beneficial in the wave theory in bioheat transfer in living tissues when the thermal wave model

(5) and (6) are used. As it is shown in Fig. 12 a longer heat penetration in skin depth can be predicted.

5. Conclusions In this paper, two methods of solving thermal wave equation are derived based on strongly continuous semigroups theory to analyse

A. Malek, G. Abbasi / Computers in Biology and Medicine 62 (2015) 65–75

300

200

References

t (s) 100

0 10

5

0

0.02

0.03

0.01

0.00

z m

Fig. 11. Elevated temperature Tðz; tÞ  T 0 ðzÞ predicted for the case of constant surface temperature heating Tð0; tÞ ¼ T 0 ð0Þ þ 12 1C, for z A ½0; l and t A ½0; 300.

Temperature T °C

40

35

30

z 8 10 5 m z 0.00208m 0

10

20

75

30

40

t s

Fig. 12. Temperature predicted by thermal wave model during sinusoidal heating Tð0; tÞ ¼ T 0 ð0Þþ 5 sin ð0:2πtÞ 1C at two different depths z ¼ 8  10  5 m and z ¼ 0:00208 m.

the various thermotherapy or skin burning in the clinical studies. In practice, results from this research can help the physician to raise the temperature of the disease's tissue to a therapeutic value, and then thermally destroy it, (see examples in Sections 4.1–4.3). In order to make a reasonable comparison, the data for conducted calculations are taken from the practitioners in the thermal engineering departments. This mathematical simulation benefits to medical therapy by clear analysis of where, when, how much, how and in which style, one may apply the heat sources on the patient's living tissue. From the mathematical point of view, when the pulse surface heat flux is applied, the hybrid method based on C0-semigroup and the corresponding variational boundary value problem are used to derive the closed form analytical solution. Different ways (exponential, point and adiabatic sources) of actual thermal loads are modelled through various initial-boundary value problems. Spurious unwanted oscillations near the discontinuity caused by numerical solutions for the parabolic–hyperbolic equations do not coincide with the problems physical properties. The closed form analytical solutions proposed in this paper do not suffer from the unwanted oscillations near and far from the discontinuities. The mathematical simulations show the efficiency of the proposed solutions for different practical cases. Conflict of interest None declared. Acknowledgments The authors wishes to express their gratitude to the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper.

[1] L. Jing, R. Zepei, W. Cuncheng, Interpretation of living tissue's temperature oscillations by thermal wave theory, Chinese Sci. Bull. 40 (17) (1995) 1493–1494. [2] A. Malek, G. Abbasi, Optimal control for Pennes' bioheat equation, Asian J. Control 17 (6) (2015) 1–12. [3] J. Liu, X. Chen, L.X. Xu, New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating, IEEE Trans. Biomed. Eng. 46 (4) (1999) 420–428. [4] A. Malek, Z. Bojdi, P. Golbarg, Solving fully 3D microscale dual phase lag problem using mixed-collocation, finite difference discretization, ASME J. Heat Transf. 134 (9) (2012) 094501–094506. [5] A. Malek, S.H. Momeni-Masuleh, A mixed collocation finite difference method for 3D microscopic heat transport problems, J. Comput. Appl. Math. 217 (1) (2008) 137–147. [6] S.H. Momeni-Masuleh, A. Malek, Hybrid pseudo spectral-finite difference method for solving a 3D heat conduction equation in a submicroscale thin film, Numer. Methods: Part D. E. 23 (5) (2007) 1139–1148. [7] 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 (3) (2009) 286–293. [8] F. Xu, M. Lin, T.J. Lu, Modeling skin thermal pain sensation: role of non-Fourier thermal behavior in transduction process of nociceptor, Comput. Biol. Med. 40 (5) (2010) 478–486. [9] M. Honner, Heat waves simulation, Comput. Math. Appl. 38 (9) (1999) 233–243. [10] K.-C. Liu, Thermal propagation analysis for living tissue with surface heating, Int. J. Therm. Sci. 47 (5) (2008) 507–513. [11] J. Liu, Preliminary survey on the mechanisms of the wave-like behaviors of heat transfer in living tissues, Forsch. Ingenieurwes. 66 (1) (2000) 1–10. [12] A. Haji-Sheikh, J. Beck, Green's function solution for thermal wave equation in finite bodies, Int. J. Heat Mass Transf. 37 (17) (1994) 2615–2626. [13] H. Ahmadikia, A. Moradi, R. Fazlali, A. Parsa, Analytical solution of non-Fourier and fourier bioheat transfer analysis during laser irradiation of skin tissue, J. Mech. Sci. Technol. 26 (6) (2012) 1937–1947. [14] H. Ahmadikia, R. Fazlali, A. Moradi, Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat flux conditions on skin tissue, Int. Commun. Heat Mass 39 (1) (2012) 121–130. [15] R. Fazlali, H. Ahmadikia, Analytical solution of thermal wave models on skin tissue under arbitrary periodic boundary conditions, Int. J. Thermophys. 34 (1) (2013) 139–159. [16] A. Malek, G. Abbasi, Optimal control solution for Pennes' equation using strongly continuous semigroup, Kybernetika 50 (4) (2014) 530–543. [17] H. Heidari, H. Zwart, A. Malek, Controllability and Stability of 3D Heat Conduction Equation in a Submicroscale Thin Film, Department of Applied Mathematics, University of Twente, Enschede, The Netherlands, 2010, pp. 1–21. [18] H. Arkin, L. Xu, K. Holmes, Recent developments in modeling heat transfer in blood perfused tissues, IEEE Trans. Biomed. Eng. 41 (2) (1994) 97–107. [19] M.N. Özisik, D.Y. Tzou, On the wave theory in heat conduction, ASME J. Heat Transf. 116 (3) (1994) 526–535. [20] D.Y. Tzou, Macro- to Microscale Heat Transfer: The Lagging Behavior, Taylor & Francis, Washington, DC, 1996. [21] K. Mitra, S. Kumar, A. Vedavarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, ASME J. Heat Transf. 117 (3) (1995) 568–573. [22] W. Kaminski, Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, ASME J. Heat Transf. 112 (3) (1990) 555–560. [23] A. Luikov, Analytical Heat Diffusion Theory, Academic Press, New York, 1968. [24] H.H. Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm, J. Appl. Physiol. 1 (2) (1948) 93–122. [25] Z.S. Deng, J. Liu, Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies, J. Biomech. Eng. 124 (2002) 638–649. [26] I. Titeux, Y. Yakubov, Application of Abstract Differential Equations to Some Mechanical Problems, vol. 558, Springer, Netherlands, 2003. [27] E. Hille, R.S. Phillips, Functional Analysis and Semi-groups, vol. 31, AMS, New York, 1957. [28] R.F. Curtain, H. Zwart, An Introduction to Infinite-Dimensional Linear Systems Theory, Text in Applied Mathematics, vol. 21, Springer-Verlag, New York, 1995. [29] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Applied Mathematical Sciences, vol. 44, Springer-Verlag, New York, 1983. [30] B.D. Reddy, Introductory Functional Analysis with Applications to Boundary Value Problems and Finite Elements, Text in Applied Mathematics, vol. 27, Springer-Verlag, New York, 1998. [31] K.R. Diller, Modeling of bioheat transfer processes at high and low temperatures, Adv. Heat Transf. 22 (1992) 157–357. [32] J.H. Li, H. Liang, Laser Medicine-Applications of Laser in Biology and Medicine (in Chinese), Science Press, Beijing, 1989.

©2015 Elsevier