International Journal of Heat and Mass Transfer 62 (2013) 153–162

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Study of the one dimensional and transient bioheat transfer equation: Multi-layer solution development and applications D.B. Rodrigues a,b,⇑, P.J.S. Pereira a,c,⇑, P. Limão-Vieira a, P.R. Stauffer b, P.F. Maccarini b a

CEFITEC, Departamento de Física, Faculdade de Ciências e Tecnologia, Universidade Nova de Lisboa, 2829 516 Caparica, Portugal Department of Radiation Oncology, Hyperthermia Division, PO BOX 3085 Duke University Medical Center, Durham, NC 27710, USA c Department of Mathematics, Instituto Superior de Engenharia de Lisboa, Rua Conselheiro Emídio Navarro 1, 1959 007 Lisboa, Portugal b

a r t i c l e

i n f o

Article history: Received 30 April 2012 Received in revised form 3 November 2012 Accepted 13 November 2012 Available online 22 March 2013 Keywords: Bioheat transfer equation Method of separation of variables Eigenfunctions and eigenvalues Bessel functions Finite Element Method Hyperthermia applications

a b s t r a c t In this work we derive an analytical solution given by Bessel series to the transient and one-dimensional (1D) bioheat transfer equation in a multi-layer region with spatially dependent heat sources. Each region represents an independent biological tissue characterized by temperature-invariant physiological parameters and a linearly temperature dependent metabolic heat generation. Moreover, 1D Cartesian, cylindrical or spherical coordinates are used to define the geometry and temperature boundary conditions of first, second and third kinds are assumed at the inner and outer surfaces. We present two examples of clinical applications for the developed solution. In the first one, we investigate two different heat source terms to simulate the heating in a tumor and its surrounding tissue, induced during a magnetic fluid hyperthermia technique used for cancer treatment. To obtain an accurate analytical solution, we determine the error associated with the truncated Bessel series that defines the transient solution. In the second application, we explore the potential of this model to study the effect of different environmental conditions in a multi-layered human head model (brain, bone and scalp). The convective heat transfer effect of a large blood vessel located inside the brain is also investigated. The results are further compared with a numerical solution obtained by the Finite Element Method and computed with COMSOL Multiphysics v4.1Ó. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The heat transfer in living tissues, known as bioheat transfer, is a complex phenomenon that depends on the thermodynamics of the biological system, its thermal constitutive parameters and the thermal response to external stimulus, e.g., electromagnetic or ultrasonic waves used in cancer treatments [1–4]. The study of bioheat transfer is especially relevant to the field of thermal medicine, since experimental temperature data is not extensively available. Temperature measurement techniques are mostly invasive as well as expensive and provide a limited number of measurement points. Non-invasive temperature measurement techniques, such as magnetic resonance thermal imaging, allow volumetric temperature measurements. However, they are limited due to its high cost and low thermal resolution [3,5]. Several therapeutic applications based on the knowledge of bioheat transfer involve either raising or lowering temperature from ⇑ Corresponding authors. Address: CEFITEC, Departamento de Física, Universidade Nova de Lisboa 2829 516 Caparica, Portugal. Tel.: +351 212 948 576x10521; fax: +351 212 948 549 (D.B. Rodrigues). E-mail addresses: [email protected], [email protected] (D.B. Rodrigues). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.11.082

normal body temperature, namely, hyperthermia [1–3] and hypothermia [4], respectively. Hyperthermia may be defined as raising the temperature of a certain region of the body above normal for a defined period of time, usually between 30 and 90 min [1]. The most common techniques to induce hyperthermia are based on heat deposition from electromagnetic [3,5] or ultrasound sources [2], where the biological tissues convert the absorbed energy into heat causing a temperature increase. Another approach to heat generation involves injecting magnetic nanoparticles immersed in fluid into the target tissue to absorb energy at a higher rate than the surrounding tissue from an externally applied electromagnetic field [6–9]. This technique is known as magnetic fluid hyperthermia (MFH). The efficacy of hyperthermia for cancer therapy is dependent on the delivery of well-controlled moderate heating (approximately 42 °C) to the entire tumor volume without overheating the surrounding critical normal tissues [3,10]. This technique takes advantage of the rapid neoplastic cell growth, which makes it more sensitive to an increase of temperature [11]. To optimize new hyperthermia based procedures, it is essential to develop a simplified but accurate model to estimate the temperature distribution and highlight the overall effect of the various parameters.

154

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

Nomenclature i r⁄ t⁄ T i T 0i T a Pi

tissue layer index spatial coordinate (m) time (s) temperature (K) initial temperature (K) arterial blood temperature (K) internal heat generation (W m3)

In 1948, Pennes was the first to propose and validate experimentally an analytical bioheat transfer model with a heat loss term due to blood perfusion [12]. Besides perfusion, Pennes’ model also accounted for thermal storage, heat conduction and heat generation caused by internal and/or external sources. Other accurate bioheat transfer models have been suggested [13,14]. However, Pennes’ model is the most widely used because of its simplicity and acceptable accuracy if no large thermally significant blood vessels are close to the analyzed heated region [12,13]. Solutions of the Pennes’ bioheat equation were obtained in regions with Cartesian, cylindrical and spherical geometries [6–8,15–18]. Durkee et al. derived an analytical solution of the classical bioheat equation using eigenfunctions for spherical and Cartesian coordinates in the reference [15] and cylindrical coordinates in reference [16]. In both cases, a constant heat source term was used. In reference [17] Durkee et al. used Green functions to solve the classical bioheat equation for time dependent sources. Continuity boundary conditions to the temperature and heat flow were imposed at the interfaces. Moreover, Neumann and Robin boundary conditions at the inner and outer surfaces were assumed, respectively. Bagaria and Johnson [7] used the method of separation of variables to obtain a transient and 1D solution to estimate the temperature in two concentric spherical regions. They assumed that the tumor was located in the inner region containing magnetic nanoparticles only with a polynomial distribution. On the other hand, a source term described by an exponential function was validated experimentally by Salloum et al. [9]. The purpose of this work is to derive an analytical solution to the transient and one-dimensional Pennes’ bioheat equation in a multi-layer region with generic spatially dependent heat sources. Each region represents an independent biological tissue (e.g., skin, fat or muscle) characterized by temperature-invariant physiological parameters and linearly temperature dependent metabolic heat generation. Moreover, 1D Cartesian, cylindrical or spherical coordinates are used to define the geometry and continuity boundary conditions are imposed to the temperature and heat flow between adjacent layers. The inner and outer surfaces satisfy equations with adaptable parameters that allow one to define Dirichlet, Neumann and/or Robin boundary conditions. This bioheat transfer model, which makes use of a formalism previously described by Rodrigues et al. [19], is applied to obtain the theoretical temperature profiles in the tumor bed and surrounding healthy tissue using two spatially dependent heat source terms to simulate a MFH treatment. We further explore the potential of this model to study the effect of different environmental conditions in a multi-layered human head model (brain, bone and scalp). The convective effect of a large blood vessel located inside the brain is also investigated assuming a laminar and fully thermally developed blood flow. Furthermore, we use two approaches to validate the analytical solution. In the first one, we determine the error associated with the truncated Bessel series that defines the transient solution

mass density (kg m3) specific heat capacity (J kg1 K1) thermal conductivity (W m1 K1) blood perfusion (s1) mass density of blood (kg m3) specific heat capacity of blood (J kg1 K1) metabolic heat generation (W m3)

qi ci ki

xbi qb cb Q mi

whereas in the second one we compare the analytical and numerical solutions. These numerical solutions are obtained using the Finite Element Method (FEM) and computed with COMSOL Multiphysics v4.1Ó. 2. Mathematical formulation 2.1. Pennes bioheat transfer equation The bioheat transfer equation in a multi-layer region is given by

qi ci

  @T i   ðr ; t Þ ¼ ki r2 T i ðr  ; t Þ þ xbi qb cb T a  T i ðr  ; t Þ @t  þ Q mi ðr  ; t  Þ þ Pi ðr  Þ

ð1Þ

where

r2 T i ðr  ; t Þ ¼ Q mi ðr  ; t Þ

¼

   1 @ G @T i   r ðr ; t Þ rG @r  @r

Q m0i

þ

Q msi T i ðr  ; t  Þ

ð2Þ ð3Þ

with 1 6 i 6 n; r i1 6 r  6 ri ðn 2 NÞ and G = 0, 1 and 2 for problems with 1D Cartesian, cylindrical and spherical symmetric geometries, respectively. Note that the mathematical method prescribed here  works for a constant Q mi ¼ Q m0i or linearly temperature dependent metabolic heat generation given by (3) with a slope denoted by Q msi [20]. The tissue described by (1)  is controlled by heat    temperature  storage qi ci @T i =@t , thermal conduction ki r2 T i , dissipation   of through blood flow ðxbi qb cb T a  T i Þ and heat generation  heat  P i , which represents the contribution from volumetric heat generation, converted from some other form of energy such as electromagnetic, ultrasonic or other modes of heating. Metabolic heat   generation Q mi is another type of heat input resulting from biochemical conversion of energy within tissue [11,12]. As example, a region with n layers and a spherically symmetric geometry is presented in Fig. 1. These layers correspond to biological tissues (e.g., skin, fat and muscle) and r 0 can be equal to zero or r0 > 0 m to take into account air or liquid body regions as well as catheters like those used in transurethral prostatic microwave thermotherapy [21]. 2.2. Boundary and initial conditions Boundary conditions of first, second and third kinds to the temperature at the inner and outer surfaces are assumed (see (4) and (5)). Temperature and heat flow must satisfy continuity boundary conditions at the tissue interfaces (see (6)–(9)). An initial spatially dependent temperature is also considered (see (10)).  Inner surface of 1st layer (i = 1)

Ain

  @T 1     r ; t þ Bin T 1 r0 ; t  ¼ C in @r  0

ð4Þ

155

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

with



 2 

ai ¼ ki = qi ci xbm rn 

bi ¼ q

ri ðrÞ

ð16Þ

  b c b bi  Q msi =ð i c i bm Þ       ¼ bi b c b T a  T r0 þ Q m0i

x x q

q x

þ

Q msi T r0

þ

Pi ðrÞ

ð17Þ     = qi ci xbm T r0 ð18Þ

To solve (15) by the method of separation of variables, the solution of this nonhomogeneous equation is decomposed in two components as follows:

T i ðr; tÞ ¼ T i ðr; tÞ þ T^ i ðrÞ

ð19Þ

where T i and T^ i are the transient (homogeneous) and steady state (nonhomogeneous) components, respectively. Note that this analytical method works for linear and homogeneous equations with constant coefficients (see (20)). Thus, we have

@T i ðr; tÞ ¼ ai r2 T i ðr; tÞ  bi T i ðr; tÞ @t ai r2 T^ i ðrÞ  b T^ i ðrÞ ¼ ri ðrÞ

ð20Þ ð21Þ

i

2.4. Boundary and initial conditions to the transient and steady state problems Fig. 1. Cross section of a region with n layers and a spherical symmetric geometry (Li  ith tissue layer).

 Outer surface of nth layer (i = n)

Aout

  @T n     r ; t þ Bout T n r n ; t  ¼ C out @r  n

ð5Þ

     r i1 ; t  ¼ T i1 ri1 ; t   @T   @T    ri1 ; t ki i ri1 ; t ¼ ki1 i1 @r @r

 Inner surface of 1st layer (i = 1)

@T 1 ðr 0 ; tÞ þ Bin T 1 ðr0 ; tÞ ¼ 0 @r dT^ 1 Ain ðr 0 Þ þ Bin T^ 1 ðr0 Þ ¼ C in dr

Ain

 Inner interface of ith layer (i – 1)

T i

The boundary and initial conditions to the transient and steady state equations (see (20) and (21)) can be obtained using (4)–(10) and (19) as well as a dimensionless analysis giving:

ð6Þ ð7Þ

ð8Þ ð9Þ

 Initial condition

T i ðr  ; t ¼ 0Þ ¼ T 0i ðr  Þ

ð10Þ

Note that parameters in (4) and (5) can be chosen in order to obtain Dirichlet, Neumann or Robin boundary conditions. 2.3. Solution methodology To obtain the dimensionless bioheat transfer equation the following relations are introduced:

  T i ðr; tÞ ¼ T i ðr; tÞ  T r0 =T r0 r¼r



ð11Þ

=rn

ð12Þ

t ¼ t xbm

ð13Þ

where xbm is the average perfusion rate within all tissue layers and T r0 is the initial reference temperature given by

T r0 ¼

n 1X 1 n i¼1 r i  r i1

Z

r i

r i1

T 0i ðr Þdr



ð14Þ

@T i ðr; tÞ ¼ ai r2 T i ðr; tÞ  bi T i ðr; tÞ þ ri ðrÞ @t

Aout Aout

@T n @r

!

d Tb n dr

þ Bout T n ðr n ; tÞ ¼ 0

ð24Þ

þ Bout Tb n ðrn Þ ¼ C out

ð25Þ

r¼r n

!

r¼r n

 Inner interface of ith layer (i – 1)

T i ðr i1 ; tÞ ¼ T i1 ðr i1 ; tÞ T^ i ðr i1 Þ ¼ T^ i1 ðr i1 Þ ! ! @T i @T i1 ¼ ki1 ki @r @r r¼ri1 r¼r ! ! i1 dT^ i dT^ i1 ¼ ki1 ki dr dr r¼ri1

ð26Þ ð27Þ ð28Þ

ð29Þ

r¼r i1

 Outer interface of ith layer (i – n)

T i ðr i ; tÞ ¼ T iþ1 ðr i ; tÞ T^ i ðr i Þ ¼ T^ iþ1 ðr i Þ ! ! @T i @T iþ1 ¼ kiþ1 ki @r @r r¼ri r¼r ! ! i dT^ i dT^ iþ1 ki ¼ kiþ1 dr dr r¼ri

Thus, one can obtain the following equation:

ð23Þ

 Outer surface of nth layer (i = n)

 Outer interface of ith layer (i – n)

  T i r i ; t  ¼ T iþ1 ðr i ; t  Þ  @T      @T   ri ; t ki i ri ; t  ¼ kiþ1 iþ1 @r @r 

ð22Þ

ð30Þ ð31Þ ð32Þ

ð33Þ

r¼ri

 Initial condition to the transient component

ð15Þ

T i ðr; t ¼ 0Þ ¼ T 0i ðrÞ  T^ i ðrÞ

ð34Þ

156

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

where

  Ain ¼ 1; Bin ¼ Bin r n =Ain ; C in ¼ C in rn =Ain T r0  Bin   Aout ¼ 1; Bout ¼ Bout r n =Aout ; C out ¼ C out rn =Aout T r0    T 0i ðrÞ ¼ T 0i ðrÞ  T r0 =T r0

T i ðr; tÞ ¼ Ri ðrÞWi ðtÞ ð35Þ  Bout ð36Þ ð37Þ

ð44Þ

The functions Ri and Wi satisfy the ordinary differential equations,

ai r2 Ri ðrÞ þ k2i Ri ðrÞ ¼ 0

ð45Þ

 dWi ðtÞ  þ bi þ ai k2i Wi ðtÞ ¼ 0 dt

ð46Þ

2.5. Analytical solution to the steady state problem

The solutions of (45) and (46) can be written as follows:

2.5.1. Generic solution to the steady state equation The solution of the linear and nonhomogeneous ordinary differential Eq. (21) is given by:

Rim ðrÞ ¼ aim r eðGÞ J eðGÞ ðkim rÞ þ bim reðGÞ Y eðGÞ ðkim rÞ    Wim ðtÞ ¼ dim exp  bi þ ai k2im t

T^ i ðrÞ ¼ Hi ðrÞ þ Si ðrÞ

where aim, bim and dim are unknown coefficients and k2im is a constant of separation. Note that Je(G) and Ye(G) are the Bessel functions of the first and second kinds, respectively. Using the superposition principle the solution of the homogeneous transient equation is defined by

ð38Þ

where Hi is the solution of the homogeneous equation,

ai r2 Hi ðrÞ  bi Hi ðrÞ ¼ 0

ð39Þ

and Si is the particular integral of (21). Solving Eq. (39), one can obtain

Hi ðrÞ ¼ ui r eðGÞ IeðGÞ ðci rÞ þ v i r eðGÞ K eðGÞ ðci rÞ

T i ðr; tÞ ¼

and ai > 0 as well as e (G) = 0.5(G  1). Note that Im(z) denotes the imaginary part of the complex number z. In Eq. (40), ui and vi are unknown coefficients and Ie(G) as well as Ke(G) are modified Bessel functions of the first and second kinds, respectively. Moreover, the order of these modified Bessel functions is denoted by e (G) and given by 1/2, 0 and 1/2 for Cartesian, cylindrical and spherical coordinates, respectively. One can show in a straightforward manner the following relationship:

dK eðGÞ dIeðGÞ 1 IeðGÞ ðci rÞ ðci rÞ  K eðGÞ ðci rÞ ðci rÞ ¼  r dr dr

ð41Þ

Using (21), (38), (40) and (41) as well as the method of variation of parameters the particular integral Si is given by

Si ðrÞ ¼

reðGÞ K eðGÞ ðci rÞ

Z

ai 

r

eðGÞ

IeðGÞ ðci rÞ

ai

r eðGÞþ1 IeðGÞ ðci rÞri ðrÞdr Z

r eðGÞþ1 K eðGÞ ðci rÞri ðrÞdr

u1

3

2

vn

0

0

0  ^ i1 /

0  ^ i2 /

ð49Þ

Z n X ki i¼1

a

2 i

ri

rG Rim ðrÞRil ðrÞdr ¼

ri1



0

if m – l

Nm

if m ¼ l

ð50Þ

with

Nm ¼

Z n X ki

a2 i¼1 i

ri

ri1

rG R2im ðrÞdr

ð51Þ

The proof of (50) is given in Appendix B. From (47)–(49) and the continuity boundary conditions (see (26), (28), (30) and (32)), one can show

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi b1  bi þ a1 k21m =ai

ð52Þ

dim ¼ dm

ð53Þ

ð42Þ

31 2 3 ^c3 0  0 0 6 7    7 7 67 6^ 7 ^ i4    0 0 7 7 6U 7 / 7 6 i7 7 7 ^ ^ i1 x ^ i2 x ^ i3 x ^ i4    0 0 7 6 x X 6 i7 7 6 7       5 45 ^c6  0 0 0 0    ^c4 ^c5

^c1 ^c2    6  7 6   7 6 6 7 6 6 6 vi 7 6 0 0    7 6 6 6u 7 ¼ 6 0 0  6 iþ1 7 6 7 6 6 4  5 4  

þ1 X Rim ðrÞWim ðtÞ

2.6.2. Eigenvalue problem One can prove that the transverse eigenfunction Rim satisfies the orthogonality condition

kim ¼

2.5.2. Calculation of the unknown coefficients ui and vi b i satisfies a system of 2n equations described by The function T the boundary conditions (23), (25), (27), (29), (31) and (33). This way, the following matricial system is deduced:

2

ð48Þ

m¼1

ð40Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where ci ¼ Im bi =ai when bi > 0 or ci ¼ bi =ai when bi 6 0

ð47Þ

0  ^ i3 /

ð43Þ where the coefficients of the matrices on the right hand side of (43) are presented in Appendix A. To solve this matricial system one can use, e.g., Maple™ software. 2.6. Analytical solution to the homogeneous transient problem 2.6.1. Method of separation of variables To solve (20) by the method of separation of variables, one can assume the following relationship:

To determine the eigenvalues k1m a system of 2n equations is defined using the boundary conditions (22), (24), (26), (28), (30), (32) and (47)–(49). To obtain solutions other than the trivial solution, one must impose



c1m





0

0





0

c2m



0









0



/i1m

0

   xi1m

/i2m  i2m x

 0

 

 0

0

 0

0

0



0









/i3m

/i4m



0

xi3m xi4m   

0

 0

 0

     c4m

0

  



0

¼0 0



  

c

5m

ð54Þ The coefficients of this determinant can be found in Appendix A. To calculate the eigenvalues the complex Newton method is used. The roots number of (54) is confirmed by the Cauchy’s argument principle. To reduce the round-off errors all numerical calculations are made with double precision using Maple™ software. 2.6.3. Determination of the coefficients aim and bim Using the eigenvalues calculated previously and the continuity boundary conditions to the temperature and heat flow (see (30) and (32)), one can establish the following recurrence formula to calculate the coefficients aim and bim:

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162



aiþ1;m biþ1;m

"

¼

 i4m  i3m / /   i4m xi3m x

#1 "

 i1m /  i1m x

157

#  i2m  aim / ; with 1 6 i 6 n  1  i2m bim x ð55Þ

From (22) and (47)–(49) the coefficients a1m and b1m are related by b1m ¼ c1 m a1m =c2 m . Note that a1m is arbitrary but any two sets of eigenfunctions obtained with different values of a1m are proportional to each other which are valid solutions to the eigenvalue problem. Furthermore, one does not need to consider a1m in the transient solution since the coefficient dm can be rewritten as ~m ¼ dm a1m (see (47)–(49)and (53)). d 2.6.4. Determination of the coefficient dm To determine the layer invariant coefficient dm the initial condition (34) and (50) are used giving

dm ¼

n 1 X ki Nm i¼1 a2i

Z

ri

r G Rim ðrÞðT 0i ðrÞ  Tb i ðrÞÞ dr

ð56Þ Fig. 2. Eigenvalues of the transient solution obtained here.

r i1

From (52) the eigenvalues can be real and/or complex numbers for a generic problem. Thus, one should consider the real part of Ti as a final and physically meaningful solution to the biological problem. To obtain a time evanescent solution to the homogeneous and transient problem for any set of physiological parameters (ai > 0 and bi > 0) and assuming real eigenvalues kim, we consider k2im as a constant of separation (see (44)–(49) and (52)). Otherwise, if one chooses k2im as a constant of separation, a different branch structure for the eigenfunctions is obtained. In this case, slightly different forms of these equations are deduced giving

Rim ðrÞ ¼ aim r eðGÞ IeðGÞ ðkim rÞ þ bim reðGÞ K eðGÞ ðkim rÞ    Wim ðtÞ ¼ dim exp  bi  ai k2im t

ð57Þ ð58Þ

The transient and homogeneous equation does not include a source term to supply energy for the biological system. Therefore, time evanescent solutions must be obtained. Note that (58) is not a time evanescent solution for any set of physiological parameters ai and bi. 3. Results and model validation 3.1. MFH technique with a polynomial heat source term To optimize the transient and 1D temperature distribution in the tumor, a spatially dependent heat source term was considered by Bagaria and Johnson [7]. They assumed concentric spherical regions to the tumor and healthy tissue. Magnetic nanoparticles were contained in the inner region (tumor) to simulate a MFH cancer treatment with a source term described by the following polynomial equation:

P1 ðr  Þ ¼ Q 1 þ Q 2 r þ Q 3 r2

After some algebra the particular integral to the steady state nonhomogeneous problem (see (18), (42) and (59)) can be evaluated giving

 expðc1 rÞ  2Q 2 þ ðsinhðc1 rÞ þ coshðc1 rÞÞ rQ 1 c21 ra1 c41     þQ 2 2 þ r 2 c21 þ Q 3 6r þ r 3 c21

S1 ðrÞ ¼

ð60Þ

  where Q j ¼ Q j =q1 c1 xb1 T r0 ; 1 6 j 6 3. From Fig. 2, one can conclude that the eigenvalues of the transient homogeneous problem satisfy approximately a linear distribution for 1 6 m 6 100 with m 2 N. This means that the time decay of the solution to the transient and homogeneous problem is more accentuated as m increases. An advantage of this linear relationship is to avoid further numerical calculations of the eigenvalues. Thus, the computational time decreases significantly. In Fig. 3, one can compare the temperature profiles obtained here (denoted by V-solution) with those of Bagaria and Johnson at the time t⁄ = 100 s and t⁄ = 300 s. The steady state temperature distributions are also compared. As expected, a good agreement between the temperature profiles is found. According to Dewey [22], the ideal therapeutic temperature is higher than 42 °C. In the case studied here, one can see that the temperature of a small part of the spherical tumor (radius 1 cm) is not in the ideal therapeutic range. However, the hyperthermia

Steady State

ð59Þ

To compare our results with those presented by Bagaria and Johnson, the following set of parameters is chosen as in [7]: q1 ¼ 1100 kg m3, q2 ¼ 1000 kg m3 ; c1 ¼ c2 ¼ cb ¼ 4200 J kg1 ; K1 k1 ¼ 0:55 W m1 K1 ; k2 ¼ 0:50 W m1 K1 ; qb xb1 ¼ qb xb2 ¼ 1 kg m3 s1 ; T a ¼ T r0 ¼ 310 K;Q m1 ¼ Q m2 ¼ 0 W m3 ; Q 1 ¼ 1  105 W m3 ; Q 2 ¼ 8  107 W m4 ; Q 3 ¼ 7  109 W m5 ; r 0 ¼ 0 m;r1 ¼ 0:01 m and r 2 ¼ 0:04 m. Note that the subscripts 1 and 2 are associated with the tumor and the surrounding tissue, respectively. The boundary conditions to the temperature at the inner and outer surfaces are defined by the following parameters: Ain ¼ Aout ¼ 0 ^ W m2 ðC ^ 2 RÞ W m1 K1 ; Bin ¼ Bout ¼ 1 W m2 K1 ; C in ¼ C and C out ¼ 310 W m2 .

Fig. 3. Comparative results between Bagaria et al. and those obtained here (Vsolution).

158

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

ei ðr; 0; MÞ ¼ T 0i ðrÞ  Tb i ðrÞ 

M X

dm Rim ðrÞ

ð63Þ

m¼1

Fig. 4. Temperature profiles inside and outside the tumor obtained here at different radius values (cm).

In Fig. 5, one can observe the absolute value of the truncation error jei(r⁄,0;M)j as a function of the radius (cm) in both layers. From this figure, one can conclude that jei (r⁄,0;M)j has a global maximum at r ¼ r 0 , which is caused due to the imposition that temperature must be bounded at this point. Since this error function is a result of combined Bessel functions with particular integral given by (60), a decay rate of jei(r⁄,0;M)j with distance is observed. However, since we impose continuity boundary conditions at r  ¼ r 1 , the error increases and as a result we find two relative maximums at both sides of this interface. Moreover, as the order M increases more Bessel functions with different arguments are combined and the number of zeros increases. As expected, as the order of the Bessel series M increases the error decreases. Thus, the value of M can be chosen depending on the accuracy required. 3.2. MFH technique with an exponential heat source term

technique induces a tumor sensitization effect caused by heat. Thus, one can take advantage of the tumor sensitization to enhance the efficacy of other techniques such as radiotherapy and chemotherapy [1,3]. Fig. 4 shows the transient temperature for different locations inside and outside the tumor. In hyperthermia treatments, both therapeutic temperature and treatment time are crucial factors that contribute to cell death. The literature clearly shows that significant cell killing as well as sensitization to radiation and/or chemotherapy occurs when cells are heated to temperatures higher than 42 °C for one hour or more [1,22]. 3.1.1. Error analysis The transient solution is given by an infinite Bessel series (see (49)) and must be truncated at some order say m = M to obtain an approximate solution. Consequently, a truncation error e(r,t;M) must be evaluated, which is defined by the expression:

T i ðr; tÞ ¼

M X

Rim ðrÞWim ðtÞ þ ei ðr; t; MÞ

ð61Þ

m¼1

From Fig. 2, one can conclude that the eigenvalues increase with increasing m and since the time dependent function Wim must be time evanescent (see (48) and (53)) the maximum truncation error occurs at t = 0. From (34) we have T i ðr; t ¼ 0Þ ¼ T 0i ðrÞ ¼ T^ i ðrÞ. Thus, one can write

jei ðr; t; MÞj 6 jei ðr; 0; MÞj where ei(r,0;M) can be rewritten as follows:

ð62Þ

Salloum et al. measured the temperature elevation in the muscle tissue of rat hind limbs induced by intramuscular injections of magnetic nanoparticles during in vivo MFH experiments [9]. Instead of a polynomial distribution, they proposed a Gaussian distribution for the internal source term:

  P1 ðr  Þ ¼ A1 exp r2 =r 2 £

ð64Þ

where r£ is a parameter that determines how far nanoparticles spread from the injection site and A1 represents the maximum strength of the volumetric heat generation. These parameters were calculated using an inverse heat transfer analysis and the bioheat equation. Similar to Salloum et al. [9], we choose the following set of parameters: A1 ¼ 169180 W m3 ; r£ ¼ 0:0078 m;qnp ¼ 5240 kg m3 ; qb ¼ q2 ¼ 1000 kg m3 ; cnp ¼ 670 J kg1 K1 ; c2 ¼ 4180 J kg1 K1 ; k1 ¼ k2 ¼ 0:5 W m1 K1 ;Q m1 ¼ Q m2 ¼ 0 W m3 ; xb1 ¼ xb2 ¼ 0:00129 s1 and T a ¼ T r0 ¼ 308:3 K. Note that the subscripts 1 and 2 refer to the muscle tissue with and without nanoparticles, respectively. The nanoparticles index is np and a constant thermal conductivity is taken into account. Moreover, a homogeneous conical geometry was used to model rat hind limbs in [9]. However, in this work we assume a cylindrical geometry with two distinct regions. The radius of the inner and outer regions are r1 ¼ 0:0078 m and r2 ¼ 0:0212 m, respectively. In this case r £ ¼ r 1 and r 0 ¼ 0 m. Here, the mass density and the specific heat capacity of the inner region are given by [8]

q1 ¼ uqnp þ ð1  uÞq2 c1 ¼ ucnp þ ð1  uÞc2

ð65Þ ð66Þ

where u = 0.0033 is the volume fraction of nanoparticles in the tissue, which is calculated using the data given in [9]. The boundary conditions to the temperature at the inner and outer surfaces are defined by the following parameters: b W m2 ðC ^ 2 RÞ A ¼ Ain ¼ 0 W m1 K1 ; Bin ¼ 1 W m2 K1 ; C in ¼ C out     k2 ; Bout ¼ hair ; C out ¼ hair T air ; T air ¼ 302:5 K and hair = 13 W m2 K1. One can show in a straightforward manner the following relationships ðm 2 N0 Þ:

W m ðrÞ ¼ Z

Fig. 5. Error analysis for the transient temperature solution at t⁄ = 0 s.

Z

r2mþ1 P1 ðrÞdr ¼ 

m P1 ðrÞ X m! 2mþ22k 2k r r 2 k¼0 k! £

ð67Þ

 2l 2mþ2 m X k1 A1 r £ m!  2 2  r 2mþ2 P1 ðrÞ X W m ðrÞ ðm!Þ r £ dr ¼ Ei 1;r =r £ þ r kðl!Þ r£ 4 4 k¼1 l¼0 ð68Þ

159

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

a theoretical analysis [7]. As criteria, we impose here that the strength of the volumetric heat generation at r0 and r 1 as well as   the average power per unit of area r 0 6 r  6 r 1 of both heat sources are the same. Using these criteria we have: Q 1 ¼ 169180 W m3 ; Q 2 ¼ 5:53  106 W m4 and Q 3 ¼ 1:05 109 W m5 . This way, one can show that the temperature curves almost coincide for all r0 6 r  6 r 2 . From (4) and the values of the parameters Ain ; Bin and C in of both source terms, the inner surface boundary conditions are defined such that the temperature is bounded at r  ¼ r0 . Therefore, the v1 coefficients in (40) are obtained imposing that the limit of the steady state and nonhomogeneous solutions (see (38)) as r⁄ approaches rþ 0 exists and is finite.

Fig. 6. Salloum experimental data [9] versus the analytical solution obtained here (V-solution).

where the function denoted integral,   2  by Ei  is the exponential  ; A1 ¼ A1 =q1 c1 xb1 T r0 and r £ ¼r £ =r 1 . P1 ðrÞ ¼ A1 exp  r2 r=r1 r £ From (18), (42), (64), (67) and (68) the particular integral of the steady state nonhomogeneous problem is deduced giving

S1 ðrÞ ¼ K 0aðc11 rÞ

þ1 X m¼0

 I0 ðac11 rÞ

c2m 1

4m ðm!Þ2

þ1 X m¼0

c2m 1

W m ðrÞ

4m ðm!Þ2

h

   R W m ðrÞ Cd ðm þ 1Þ  ln c21 r þ

W m ðrÞ dr r

i

ð69Þ where Cd is the digamma function. In Fig. 6, the steady state temperature distribution derived here is compared with that of Salloum et al. [9] obtained by experimental measurements. From the symmetry of the problem, the steady state temperature is an even function and it is presented as a function of the distance between medial and lateral skin of the limb. The graph of this even function is translated of r2 units to the right in the abscissas axis in order to correctly describe this skin-to-skin distance. A good agreement between the analytical and experimental results is achieved in the full region of analysis. The results of Figs. 3 and 6 suggest that the temperature distribution associated with the exponential and polynomial heat source terms can be similar for a certain set of parameters (see (59) and (64)). In order to obtain similar temperature profiles, we recalculate Bagaria’s parameters (see (59)) since these were obtained from

hb l = 0 W m − 2 K − 1

3.3. Temperature distribution in a multi-layered human head with the effect of large blood vessels In this final example, the model developed in this work is used to study the steady state temperature distribution in the human head. The analytical solution is compared with a FEM numerical solution implemented in COMSOL Multiphysics v4.1Ó using a parallel sparse direct linear solver (PARDISO) with a mesh constituted by 4090 elements (maximal element size of 91.3 lm). Four layers are   brain white matter r1 ¼ 6:7 cm , brain gray considered:     matter r 2 ¼ 8:5 cm , bone r 3 ¼ 8:9 cm and scalp r 4 ¼ 9:3 cm . The outer surface is surrounded by air at T air ¼ 302:5 K and different convection coefficients hair are compared: 0, 5 and 15 W m 2 K1. From these parameters we can deduce the outer boundary condition from Bout ¼ hair ; C out ¼ hair T air and Aout ¼ kscalp , where kscalp = 0.49 W m1 K1 is the thermal conductivity of the scalp [23]. To further explore the potential of this model, we also study the convective effect of a large blood vessel located inside the brain with radius rbl = 0.17 cm. In this case, the inner radius r0 has to be the same as rbl in order to account for the presence of the vessel. A constant convective heat transfer coefficient is assumed at the vessel surface and the blood is considered to be at 37 °C, which corresponds to body core temperature. Assuming a laminar and fully thermally developed blood flow, the convective heat transfer coefficient is given by the following expression [24]:

hbl ¼ 1:83kbl =r bl

ð70Þ

where kbl = 0.5 W m1 K1 is the thermal conductivity of blood. All the other physical and physiological properties of the biological tissues can be found in [23].

h air = 0 W m − 2 K

−1

h air = 5 W m − 2 K

−1

hbl = 538.23 W m − 2 K − 1

h air = 15 W m − 2 K − 1

Fig. 7. Temperature in a multi-layered model of the head with different air and blood convection coefficients, where the layers correspond to: brain white matter (0–6.7 cm), brain gray matter (6.7–8.5 cm), bone (8.5–8.9 cm) and scalp (8.9–9.3 cm). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

160

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

The three cases shown in Fig. 7 model three distinguishable typical daily situations: an isolated environment as when we use an insulating hat (hair = 0 W m2 K1), an interior room environment with low but steady air circulation, e.g., an operating room (hair = 5 W m2 K1) and a windy exterior environment (hair = 15 W m2 K1). A first conclusion that can be drawn is the ability of the human head to preserve the brain temperature at 37 °C regardless of air convection coefficients used. The inclusion of a large vessel in this model provides the opportunity to analyze the effect of convection from blood vessels that are larger than arterioles and venules. This is not accounted in the classic Pennes’ model. The steady state temperature of white matter is 37.3 °C without the presence of large blood vessels (see Fig. 7). Moreover, the presence of the blood vessel of rbl = 0.17 cm located in the center of the brain, causes that the temperature near the wall drops down to 37 °C. This model aims to study the effect of a single vessel and not the complex blood vessel network in the brain [25]. However, the uniform white matter temperature of 37.3 °C (without the blood vessel) and the drop of 0.3 °C near the blood vessel, suggest that the highly vascularized brain should be at 37 °C. This conclusion agrees with the results presented in Konstas [23]. From Fig. 7, one can observe a good agreement between analytical and FEM numerical solutions. In general, numerical solutions have enough accuracy and can be developed for more complex geometries. The geometry flexibility and no development of complex mathematical formalisms are perhaps the major advantages of this kind of solutions. However, there are cases where convergence and stability problems arise leading to the divergence and blow-up of these solutions. It is well known that the violation of a CFL (Courant-Friedrichs-Lewy) type condition leads to unstable modes [26]. In some situations, one can avoid these problems by increasing the number of iterations, i.e., smaller mesh element sizes and time steps are required. Thus, the computational time increases accordingly. This may imply a non-acceptable running time of numerical algorithms. These unstable modes may also be intrinsic to the partial differential equations or caused by more complex geometries especially with steepest gradients. Note that the spatial and time dependence as well as the type of some terms or coefficients in the equations can modify the nature, i.e., the convergence and stability as well as the linearity of the heat transfer system. Therefore, these divergent and unstable modes cannot be avoided in some more complex problems for all space and time domain. Moreover, the numerical solutions are approximate in its nature, meaning quantitative, and do not provide qualitative information about the problem. From a mathematical point of view, both quantitative and qualitative analysis of solutions are preferable for more understanding of complex problems. In fact, the analytical solutions allow the study of both quantitative and qualitative response of the biological system to a given heat source term or coefficient in the partial differential equations. This is highlighted, e.g., in the Eqs. (60) and (69) where the temperature response to a polynomial and exponential heat sources is given by a combination of explicit and well defined mathematical functions. On the other hand, analytical solutions also allow real and space–time dependences of temperature and its gradients. These insights are given by an analytical approach and ultimately lead to a better understanding of the biological problem.

terms are considered. Boundary conditions of first, second and third kinds to the temperature at the inner and outer surfaces are obtained from an appropriate choice of parameters. The biological model is further enhanced by introducing an important thermoregulation mechanism, namely, the temperature dependent metabolic heat rate. The complexity of the problem presented here is highly reduced by the layer-independent coefficient, dm, and by the dependence of all eigenvalues on the one associated with the first layer. These eigenvalues also show a linear dependence on the order of the infinite Bessel series denoted by m. This linear dependence reduces the computational time since no numerical computation is required. Moreover, the mathematical formalism developed here shows that the full solution is defined by two branches, where one is time evanescent and the other is not. We stress that time evanescent solutions must be obtained since the transient and homogeneous equation does not have an external heat source term to supply energy to the biological system. Thus, the number of branches is reduced from two to one. This means that after a long time the full solution does not depend on time and is given approximately by the solution of the steady state (nonhomogeneous) equation which only depends on the spatial coordinate. The derived solution is applied to study two distinct spatially dependent heat sources that characterize the heating induced by magnetic nanoparticles subjected to an electromagnetic field. The closed-form of the analytical solution for such polynomial and exponential heat sources is deduced. A good agreement is achieved between the results provided by these analytical solutions and those obtained theoretically and experimentally by other authors. It is proved that the polynomial and exponential profiles can be interchanged by choosing the correct parameters associated with these source terms. We can conclude that these polynomial and exponential functions are equivalent and valid to describe the heat source in magnetic fluid hyperthermia. In hyperthermia treatments, both therapeutic temperature and treatment time are crucial factors in establishing cell death. The spatial transient model developed here is suitable for analyzing both these physical quantities. It is difficult to achieve an ideal temperature profile within the entire tumor and its healthy tissue surroundings, but tumor sensitization can still be induced. Thus, one can take advantage of the tumor sensitization to apply other techniques, e.g., radiotherapy or chemotherapy, in order to obtain efficacious cancer therapy. The inability to generate and maintain ideal temperature distributions in heterogeneous perfused tissue is the reason why hyperthermia is considered as an adjuvant rather than standalone cancer therapy. The versatile inner and outer boundary conditions introduced here are used to study different environmental conditions as well as the effect of a large blood vessel inside a human brain, which complements the classical bioheat model. This study shows the ability of the human head to preserve the brain temperature around 37 °C in a wide range of environmental conditions. The blood vessel could only be accounted due to the flexible inner boundary condition, which allows to have r0 = 0 or r0 > 0. As future work one can model a problem with an inner surface to simulate a hyperthermia or ablation catheter immersed in a diseased tissue to study the efficacy of different heating techniques.

4. Conclusions An analytical solution to the Pennes’ bioheat equation is derived to calculate the heat transfer in a multi-layer perfused tissue. Multi-layer regions with 1D Cartesian, cylindrical and spherical symmetric geometries as well as spatially dependent heat source

Appendix A. A-Definition of the coefficients in the Eqs. (43) and (54) The coefficients of the matrices on the right hand side of (43) are defined by

161

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

^c1 ¼ Ain ðr eðGÞ IeðGÞ ðc1 rÞÞ0r¼r þ Bin r0eðGÞ IeðGÞ ðc1 r0 Þ 0 eðGÞ

0 1 rÞÞr¼r 0

eðGÞ Bin r 0 K eðGÞ ð 1 r 0 Þ

K eðGÞ ðc

þ

^c3 ¼ C in 

Ain ðS1 ðrÞÞ0r¼r0

 Bin S1 ðr 0 Þ

^c4 ¼ Aout ðr

eðGÞ

^c5 ¼ Aout ðr

eðGÞ

^c2 ¼ Ain ðr

^c6 ¼ C out  ^ i1 ¼ /

IeðGÞ ðc

0 n rÞÞr¼r n

þ

c

Bout r neðGÞ IeðGÞ ð n ; r n Þ

c

c

c

eðGÞ ri IeðGÞ ð i ri Þ

c

ðA:4Þ ðA:5Þ ðA:6Þ ðA:7Þ

^ i2 ¼ r eðGÞ K eðGÞ ðc r i Þ / i i ^ i3 ¼ /

eðGÞ ri IeðGÞ ð iþ1 r i Þ

^ i4 ¼ /

eðGÞ ri K eðGÞ ð iþ1 ri Þ

ðA:8Þ

c

ðA:9Þ

c

ðA:11Þ

0 i rÞÞr¼r i

IeðGÞ ðc

ðA:12Þ

^ i2 ¼ ki ðr eðGÞ K eðGÞ ðci rÞÞ0r¼r x i

ðA:13Þ

^ i3 ¼ kiþ1 ðr eðGÞ IeðGÞ ðciþ1 rÞÞ0r¼r x i

ðA:14Þ

^ i4 ¼ kiþ1 ðr eðGÞ K eðGÞ ðciþ1 rÞÞ0r¼r x i

ðA:15Þ

^ i ¼ kiþ1 ðSiþ1 ðrÞÞ0  ki ðSi ðrÞÞ0 X r¼r i r¼r i

ðA:16Þ

where the prime symbol denotes differentiation with respect to r. The coefficients of the determinant in (54) are given by

c1m ¼ Ain ðr

eðGÞ

J eðGÞ ðk1m rÞÞ0r¼r0

þ

eðGÞ Bin r0 J eðGÞ ðk1m r 0 Þ eðGÞ

c2m ¼ Ain ðr eðGÞ Y eðGÞ ðk1m rÞÞ0r¼r0 þ Bin r 0 c4m ¼ Aout ðr

eðGÞ

J eðGÞ ðknm rÞÞ0r¼r0

þ

Y eðGÞ ðk1m r 0 Þ

ðA:18Þ ðA:19Þ

eðGÞ

eðGÞ

/i2m ¼ /i3m ¼

Y eðGÞ ðknm r n Þ

J eðGÞ ðkim ri Þ

eðGÞ ri Y eðGÞ ðkim r i Þ eðGÞ ri J eðGÞ ðkiþ1;m r i Þ eðGÞ

/i4m ¼ r i

ðA:17Þ

eðGÞ Bout r 0 J eðGÞ ðknm rn Þ

c5m ¼ Aout ðreðGÞ Y eðGÞ ðknm rÞÞ0r¼r0 þ Bout r 0 /i1m ¼ r i

Y eðGÞ ðkiþ1;m r i Þ

xi1m ¼ ki ðreðGÞ JeðGÞ ðkim rÞÞ0r¼ri xi2m ¼ ki ðreðGÞ Y eðGÞ ðkim rÞÞ0r¼ri xi3m ¼ kiþ1 ðreðGÞ JeðGÞ ðkiþ1;m rÞÞ0r¼ri xi4m ¼ kiþ1 ðreðGÞ Y eðGÞ ðkiþ1;m rÞÞ0r¼ri

ðA:20Þ ðA:21Þ ðA:22Þ ðA:23Þ ðA:24Þ ðA:25Þ ðA:26Þ ðA:27Þ ðA:28Þ

Appendix B. Proof of the orthogonality condition (50) Let Rim and Ril ðm; l 2 NÞ be transverse eigenfunctions satisfying Eq. (45). Thus, we have

1 rG 1 rG

  d dRim k2 rG þ im Rim ¼ 0 dr dr ai   d k2il G dRil r þ Ril ¼ 0 dr dr ai

ðB:1Þ

!     Ril d Rim d k2im k2il G dRim G dRil R R ¼0 r r   G þ rG dr dr r dr dr ai ai il im

Ril

ri1

þ

Z

ri

r i1

    Z ri d dRim d dRil dr  dr rG rG Rim dr dr dr dr r i1 ! k2 k2 r G im  il Ril Rim dr ¼ 0

ai

ai

ðB:5Þ

Multiplying the above equation by ki and summing over all i(1 6 i 6 n), one can deduce

r n  X dRim dRil i ki rG Ril  ki rG Rim dr dr ri1 i¼1 Z n r X ki  2  i G þ kim  k2il r Rim Ril dr

ai

i¼1

ri1

¼0

ðB:6Þ

From (B.6) and applying interface boundary conditions (26), (28), (30) and (32) as well as (47)–(49), we have

h

kn rG Rnl dRdrnm  kn r G Rnm dRdrnl

i r¼r n

h i  k1 r G R1l dRdr1m  k1 rG R1m dRdr1l

r¼r0

þ

n X  2  R ri G ki 2 þ a kim  kil ri1 r Rim Ril dr ¼ 0 i¼1

i

ðB:7Þ From the outer surface boundary condition (24) and (47)–(49), we obtain

 dRnm Aout ¼0 þ Bout Rnm dr r¼r n  dRnl Aout ¼0 þ Bout Rnl dr r¼r n

ðB:8Þ ðB:9Þ

Multiplying (B.8) by kn r Gn Rnl ðr n Þ and (B.9) by kn r Gn Rnm ðrn Þ and subtracting, one can show

 dRnm dRnl Aout kn rGn Rnl ¼0  Rnm dr dr r¼rn

ðB:10Þ

Now, different cases can be analyzed. If Aout – 0 (B.10) can be rewritten as follows:

 dRnm dRnl kn rGn Rnl ¼0  Rnm dr dr r¼rn

ðB:11Þ

Note that (B.11) is valid for Aout – 0 and Bout – 0 as well as Aout – 0 and Bout = 0. Moreover, if Aout = 0 and Bout – 0, one can obtain Rnm(rn) = 0 and Rnl (rn) = 0 using (B.8) and (B.9), respectively. Thus, (B.11) is also valid for this last case. A similar methodology can be applied using the inner surface boundary condition (22) and (47)–(49), giving

 dR1m dR1l k1 rG0 R1l ¼0  R1m dr dr r¼r0

n X ki  i¼1

ðB:12Þ

ðB:4Þ

ai

k2im  k2il



Z

ri

rG Rim Ril dr ¼ 0

ðB:13Þ

r i1

From (52) and a similar expression to kil, one can show

ðB:3Þ

Multiplying (B.3) by rG and integrating from ri1 to ri, one can show ri

!  r Z ri dRim dRil i k2 k2 r G Ril þ r G im  il Ril Rim dr ¼ 0  r G Rim dr dr ri1 ai ai r i1

Finally, one can use Eqs. (B.7), (B.11) and (B.12) to obtain

ðB:2Þ

Multiplying (B.1) by Ril and (B.2) by Rim and subtracting, we obtain

Z

Applying the method of integration by parts twice in the first integral of (B.4), the following expression is obtained:

ðA:10Þ

^ i ¼ Siþ1 ðr i Þ  Si ðr i Þ U ^ i1 ¼ ki ðr x

ðA:2Þ ðA:3Þ

K eðGÞ ð n rÞÞ0r¼rn þ Bout r neðGÞ K eðGÞ ð n ; r n Þ Aout ðSn ðrÞÞ0r¼rn  Bout Sn ðrn Þ

eðGÞ

ðA:1Þ

k2im  k2il ¼

 a1  2 k  k2 ai 1m 1l

ðB:14Þ

Note that k21m  k21l – 0 for m – l. Using this last result as well as (B.13) and (B.14), one can deduce the orthogonality condition (50) giving

Z n X ki i¼1

a

2 i

ri

r i1

r G Rim ðrÞRil ðrÞdr ¼ 0 ðm – lÞ

ðB:15Þ

162

D.B. Rodrigues et al. / International Journal of Heat and Mass Transfer 62 (2013) 153–162

References [1] R.W.Y. Habash, D. Krewski, R. Bansal, H.T. Alhafid, Principles, applications, risks and benefits of therapeutic hyperthermia, Front. Biosci. 3 (2011) 1169–1181 (Elite edition). [2] T.W.H. Sheu, M.A. Solovchuk, A.W.J. Chen, M. Thiriet, On an acoustics-thermalfluid coupling model for the prediction of temperature elevation in liver tumor, Int. J. Heat Mass Transfer 54 (2011) 4117–4126. [3] O.I. Craciunescua, P.R. Stauffer, B.J. Soher, C.R. Wyatt, O. Arabe, P. Maccarini, S.K. Das, K.S. Cheng, T.Z. Wong, E.L. Jones, M.W. Dewhirst, Z. Vujaskovic, J.R. MacFall, Accuracy of real time noninvasive temperature measurements using magnetic resonance thermal imaging in patients treated for high grade extremity soft tissue sarcomas, Med. Phys. 36 (2009) 4848–4858. [4] Y.J. Wang, L. Zhu, A.J. Rosengart, Targeted brain hypothermia induced by an interstitial cooling device in the rat neck: experimental study and model validation, Int. J. Heat Mass Transfer 51 (2008) 5662–5670. [5] J. Gellermann, W. Wlodarczyk, A. Feussner, H. Fahling, J. Nadobny, B. Hildebrandt, R. Felix, P. Wust, Methods and potentials of magnetic resonance imaging for monitoring radiofrequency hyperthermia in a hybrid system, Int. J. Hyperther. 21 (2005) 497–513. [6] M.A. Giordano, G. Gutierrez, C. Rinaldi, Fundamental solutions to the bioheat equation and their application to magnetic fluid hyperthermia, Int. J. Hyperther. 26 (2010) 475–484. [7] H.G. Bagaria, D.T. Johnson, Transient solution to the bioheat equation and optimization for magnetic fluid hyperthermia treatment, Int. J. Hyperther. 21 (2005) 57–75. [8] C.T. Lin, K.C. Liu, Estimation for the heating effect of magnetic nanoparticles in perfused tissues, Int. Commun. Heat Mass Transfer 36 (2009) 241–244. [9] M. Salloum, R.H. Ma, L. Zhu, An in-vivo experimental study of temperature elevations in animal tissue during magnetic nanoparticle hyperthermia, Int. J. Hyperther. 24 (2008) 589–601. [10] P. Schildkopf, O.J. Ott, B. Frey, M. Wadepohl, R. Sauer, R. Fietkau, U.S. Gaipl, Biological rationales and clinical applications of temperature controlled hyperthermia-implications for multimodal cancer treatments, Curr. Med. Chem. 17 (2010) 3045–3057. [11] C. Song, I. Choi, B. Nah, S. Sahu, J. Osborn, Microvasculature and Perfusion in Normal Tissues and Tumors, Seegenschmiedt Mh, Thermoradiotherapy and Thermochemotherapy, Springer-Verlag, Berlin/New York, 1995. [12] H.H. Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm, J. Appl. Physiol. 1 (1948) 93–122.

[13] M. Stanczyk, J. Telega, Modelling of heat transfer in biomechanics-a review Part.I soft tissues, Acta Bioeng. Biomechan. 4 (2002) 31–61. [14] J. Fan, L. Wang, A general bioheat model at macroscale, Int. J. Heat Mass Transfer 54 (2011) 722–726. [15] J.W. Durkee, P.P. Antich, C.E. Lee, Exact-solutions to the multiregion timedependent bioheat equation 1: solution development, Phys. Med. Biol. 35 (1990) 847–867. [16] J.W. Durkee, P.P. Antich, Characterization of bioheat transport using an exact solution to the cylindrical geometry, multiregion, time-dependent bioheat equation, Phys. Med. Biol. 36 (1991) 1377–1406. [17] J.W. Durkee, P.P. Antich, Exact-solutions to the multiregion time-dependent bioheat equation with transient heat-sources and boundary-conditions, Phys. Med. Biol. 36 (1991) 345–368. [18] J. Okajima, S. Maruyama, H. Takeda, A. Komiya, Dimensionless solutions and general characteristics of bioheat transfer during thermal therapy, J. Therm. Biol. 34 (2009) 377–384. [19] D.B. Rodrigues, P.J.S. Pereira, P.M. Limão-Vieira, P.F. Maccarini, Analytical solution to the transient 1D bioheat equation in a multilayer region with spatial dependent heat sources, in: Proceedings of the Eighth IASTED International Conference on Biomedical Engineering, Acta Press, Innsbruck, 2011, pp. 96–103. [20] K.N. Rai, S.K. Rai, Heat transfer inside the tissues with a supplying vessel for the case when metabolic heat generation and blood perfusion are temperature dependent, Heat Mass Transfer 35 (1999) 345–350. [21] L. Zhu, L.X. Xu, N. Chencinski, Quantification of the 3-D electromagnetic power absorption rate in tissue during transurethral prostatic microwave thermotherapy using heat transfer model, IEEE Trans. Biomed. Eng. 45 (1998) 1163–1172. [22] W.C. Dewey, Arrhenius relationships from the molecule and cell to the clinic, Int. J. Hyperther. 25 (2009) 3–20. [23] Konstas, M. Neimark, A. Laine, J. Pile-Spellman, A theoretical model of selective cooling using intracarotid cold saline infusion in the human brain, J. Appl. Physiol. 102 (2007) 1329–1340. [24] J.J.W. Lagendijk, The influence of bloodflow in large vessels on the temperature distribution in hyperthermia, Phys. Med. Biol. 27 (1982) 17–23. [25] F.H. Netter, Atlas of Human Anatomy, Elsevier, Philadelphia, 2006. [26] N.D. Lopes, P.J.S. Pereira, L. Trabucho, A numerical analysis of a class of generalized Boussinesq-type equations using continuous/discontinuous FEM, Int. J. Numer. Meth. Fluids 69 (2012) 1186–1218.