Finite integral transform-based analytical solutions of dual phase lag bio-heat transfer equation

Accepted Manuscript

Finite integral transform-based analytical solutions of dual phase lag bio-heat transfer equation Sumit Kumar , Atul Srivastava PII: DOI: Reference:

S0307-904X(17)30371-2 10.1016/j.apm.2017.05.041 APM 11793

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

27 July 2016 12 May 2017 23 May 2017

Please cite this article as: Sumit Kumar , Atul Srivastava , Finite integral transform-based analytical solutions of dual phase lag bio-heat transfer equation, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.05.041

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights  Analytical solution of non-Fourier bio heat transfer equation. Finite integral transform.



Time independent/dependent boundary conditions.



Coupling of transient RTE and non-Fourier heat conduction models.



Thermal analysis of laser-irradiated biological tissue phantoms.

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT Finite integral transform-based analytical solutions of dual phase lag bio-heat transfer equation Sumit Kumar, Atul Srivastava* Department of Mechanical Engineering, Indian Institute of Technology Bombay, Mumbai 400076, India

Abstract Finite integral transform (FIT)-based analytical solution of dual phase lag (DPL) bio-heat

CR IP T

transfer equation has been developed. One of the potential applications of the developed analytical approach is in the field of photo-thermal therapy wherein the interest is in determining the thermal response of the laser-irradiated biological samples. In order to demonstrate the applicability of the generalized analytical solutions, three problem definitions have been formulated; (1) time independent boundary conditions (constant surface

AN US

temperature heating), (2) time dependent boundary conditions (medium subjected to sinusoidal surface heating) and (3) biological tissue phantoms subjected to short-pulse laser irradiation. In the context of the case study involving biological tissue phantoms, the FITbased analytical solutions of Fourier as well as non-Fourier heat conduction equations have been coupled with the numerical solution of transient form of radiative transfer equation

M

(RTE) to determine the resultant temperature distribution. Performance of the FIT-based approach has been assessed by comparing the results of the present study with those reported

ED

in the literature. A comparison of DPL-based analytical solutions with those obtained using the conventional Fourier and hyperbolic heat conduction models has been presented. Relative influence of relaxation times associated with the temperature gradients (τT) and heat flux (τq)

PT

on the resultant thermal profiles has also been discussed. To the best of the knowledge of the authors, the present study is the first successful attempt towards developing complete FIT-

CE

based analytical solution(s) of the non-Fourier heat conduction equation(s), which have subsequently been coupled with the numerical solution of transient form of RTE. The work

AC

finds its importance in a range of areas such as material processing, photo-thermal therapy etc.

Keywords: Dual phase lag model, Finite integral transform technique, Analytical solution, Bio-heat transfer.

*

Corresponding author. Tel.: +91 22 25767531; fax: +91 22 25726875. E-mail addresses: [email protected], [email protected] (A. Srivastava)

2

ACCEPTED MANUSCRIPT

Nomenclature constant

Greek symbols

cv

specific heat of the tissue

γ

characteristic root

C, F, G function

θ

temperature difference, T-Tb

H

inverse of length, h/k

̃̅

transformed temperature

h

heat transfer coefficient

ν, η

eigenvalues

k

thermal conductivity

ρ

density of the medium

L

length

ω

blood perfusion/characteristic root

N

normalization integral

QM

metabolic heat generation

Subscripts

Qr

external heat source

b

blood

⃗

heat flux

H

transient

amplitude of sinusoidal profile temperature

t

Time

W

width

X, Y

eigenfunction

x, y

cartesian spatial coordinates

ED

1. Introduction

M

metabolic

NH

particular solution

n

previous time instant

n+1

current time instant

p, m

integer

r

external heat source

SS

steady state

M

T

AN US

constant heat flux

CR IP T

A, B

PT

For the past two centuries, Fourier law of heat conduction has found its application in a range of engineering and scientific problems for understanding the phenomena of heat transfer.

CE

However, it has been demonstrated by a group of researchers that the standard Fourier law of heat conduction has its limitations in applications involving small-length scales (e.g. thin film

AC

deposition), extremely high heat fluxes and temperature gradients, laser-material interaction wherein short duration laser pulses are employed, non-homogenous structures such as biological samples etc. [1-6]. One of the primary reasons for these limitations is the inherent assumption of infinite speed of thermal front propagation through the medium, and hence the model assumes that the heat flux and temperature gradients setup across a material volume occur at the same instant of time [7]. However, in applications involving samples with nonhomogeneous internal structures, e.g. biological samples, it has been experimentally demonstrated that the speed of thermal wave propagation through the medium is finite [8, 9]

3

ACCEPTED MANUSCRIPT and hence the Fourier law of heat conduction cannot accurately predict the thermal response of such samples. In order to circumvent the above mentioned limitations of the conventional Fourier model-based heat conduction equation, a range of researchers have proposed certain modifications in this basic model. In one of the major developments, Cattaneo [10] and Vernotte [11] introduced an additional phase lag term for the heat flux. The resultant model

conduction equation and is expressed as ⃗ Although, the introduction of

⃗⃗

CR IP T

was subsequently named as Cattaneo-Vernotte heat conduction equation or hyperbolic heat

(1)

in Equation (1) takes care of the assumption of

infinite thermal wave propagation speed, the modified form of equation does not consider the

AN US

effects of micro-structural interactions in non-homogeneous medium such as biological tissue samples. In view of these limitations of Cattaneo-Vernotte model, Tzou [12] modified the hyperbolic heat conduction equation by introducing another phase lag term for temperature gradients, and the resultant final form has been termed as dual phase lag (DPL) model, which can be expressed in the following form of equation: ,

denotes the phase lag time for heat flux and

-

(2)

corresponds to the phase lag time

ED

where

⃗⃗

M

⃗

associated with the temperature gradients. The DPL model, given by Equation (2), is considered to be the generalized form of non-Fourier heat conduction model. In the limits , the DPL model transforms into hyperbolic model, while for

PT

of

, one gets

its conversion into the standard Fourier law of heat conduction.

CE

Among other application areas such as nanoscale transient heat transfer, rapid solidification of a pure metal particle, two-layer slab with imperfect interface, ultrafast laser

AC

heating of gold films [13-16], the non-Fourier heat conduction models (DPL/hyperbolic) have been employed for investigating the thermal response of laser-irradiated biological samples in the context of photo-thermal therapy. Efforts have been made by various researchers towards developing appropriate numerical and/or analytical models for solving the DPL/hyperbolic heat conduction equations for the determination of temperature distribution. A survey of literature reveals that these non-Fourier heat conduction equations have primarily been solved using numerical approaches and some of the notable studies in this context have now been discussed. Zhou et al. [17] numerically investigated the thermal damage of a laser-irradiated

4

ACCEPTED MANUSCRIPT biological tissue sample using one-dimensional DPL model. The non-Fourier bio-heat transfer equation was discretized using finite volume method (FVM). Another study was reported by Zhou et al. [18] in which DPL model based bio-heat transfer equation in twodimensional axisymmetric biological tissues subjected to surface laser heating or body laser heating was solved using FVM to determine the temperature distribution. Along similar lines, Narasimhan and Sadasivam [19] used FVM for solving the DPL model based bio-heat transfer equation for one-dimensional model of human eye during retinal laser-irradiation.

CR IP T

Udayraj et al. [2] analysed heat transfer phenomena through fabric-human skin system subjected to flash fire as well as radiant heat exposure. The authors employed finite difference method (FDM) scheme for solving the DPL model based bio-heat transfer equation. In recent works of the authors of the present study, the DPL based bio-heat transfer equation was solved using FVM [20] and Lattice Boltzmann Method (LBM) [21] for

AN US

understanding the thermal response of laser-irradiated two-dimensional tissue phantoms. A group of researchers have also made use of hybrid methods for transforming the transient problems into their steady form using the technique of Laplace transform and then solved the steady state problem using standard numerical methods such as FVM, FDM, etc.

M

The solution thus obtained in terms of Laplace variable is then converted into physical domain numerically or using the standard tables available for the inversion of Laplace

ED

transform. Chen and Lin [22] solved the one-dimensional hyperbolic heat conduction model using one such hybrid scheme that is based on Laplace transform and control volume method for determining the numerical oscillation free temperature field distribution. Liu [23]

PT

investigated the thermal behaviour of one-dimensional skin sample subjected to constant, sinusoidal or step surface heating within the frame work of hyperbolic model based bio-heat

CE

transfer equation. Lee et al. [24] applied conjugate gradient method based inverse algorithm to solve one-dimensional inverse hyperbolic heat conduction problem for estimating the

AC

unknown time-dependent surface heat flux in living skin tissue from the known temperature history at the measurement locations. The authors solved the direct problem using a hybrid scheme for temperature determination when the time-dependent heat flux, properties and initial and boundary conditions are known. Liu and Cheng [25] investigated the thermal response of tissue samples subjected to instantaneous surface heating. The authors solved the hyperbolic bio-heat transfer equation using Laplace transform-based hybrid method. Liu and Chen [26] made use of DPL model for understanding the thermal behaviour of biological tissue during magnetic hyperthermia treatment. The authors employed the hybrid numerical scheme based on the Laplace transform, change of variables, and the modified discretization 5

ACCEPTED MANUSCRIPT technique in conjunction with hyperbolic shape function for solving the one-dimensional DPL based bio-heat transfer problem. In one of the recent studies, Liu and Wang [27] employed the hybrid method to solve the one-dimensional DPL model based bio-heat transfer equation for investigating the thermal response of laser-irradiated biological samples. The wide application of numerical approaches and/or Laplace transform-based hybrid schemes for solving the generalized non-Fourier bio-heat transfer equations, as discussed above, has primarily been motivated by the level of complexity associated with the governing

CR IP T

equations, due to which formulating the complete analytical solution is highly challenging. In spite of these challenges, some efforts have still been made to develop analytical solutions for the Fourier/non-Fourier forms of bio-heat equations. In this context, Liu et al. [28] developed analytical solution for the one-dimensional Fourier-based bio-heat transfer equation for skin subjected to sinusoidal heating. Shih et al. [29] obtained the exact solution of the Pennes bio-

AN US

heat transfer equation in one-dimensional semi-infinite skin surface subjected to sinusoidal heating using Laplace transform. Liu et al. [30] solved the one-dimensional hyperbolic heat conduction model based bio-heat transfer equation using separation of variables to study the phenomena of heat transfer in skin sample subjected to instantaneous heating. Saleh and Al-

M

Nimr [31] employed the variational formulation of Laplace transforms to solve the hyperbolic heat conduction equation. Ahmadikia et al. [32] solved one-dimensional Fourier and non-

ED

Fourier (hyperbolic) models based bio-heat transfer equations for skin samples using Laplace transform. The skin samples were treated as semi-infinite and finite domains and were subjected to constant, periodic and pulse train heat flux boundary conditions. Recently,

PT

Askarizadeh and Ahmadikia [33] developed the analytical solution of one-dimensional DPL model based bio-heat transfer equation using Laplace transform to investigate the thermal

CE

behaviour of skin tissues subjected to constant, periodic and pulse train heat flux boundary conditions. In another analytical study, Askarizadeh and Ahmadikia [34] worked with two-

AC

dimensional Fourier and non-Fourier (hyperbolic and DPL) models based bio-heat transfer equations in cylindrical coordinate systems. The authors solved these equations using Laplace transform for determining the thermal response of skin tissue surface exposed to constant and pulse train heat flux. Poor et al. [35] employed the separation of variables in combination with Duhamel’s theorem for solving the one-dimensional DPL bio-heat transfer equation to determine the temperature distribution inside skin tissue subjected to different boundary conditions such as constant, cosine and pulse train heat flux conditions. Torabi and Zhang [36] obtained the analytical solution of two-dimensional DPL model based heat conduction equation for cylindrical geometry subjected to constant and time varying heat flux boundary 6

ACCEPTED MANUSCRIPT conditions using Duhamel’s theorem in conjunction with the concept of separation of variables. Lam [37] presented a generalized analytical solution in the form of an infinite series for determining the temperature response of a finite thin film that was subjected to heating by a time-varying and spatially-decaying internal heat source. The solution was obtained using the superposition technique in conjunction with solution structure theorems and Fourier series expansion. The authors showed that using the developed generalized analytical solution, one can generate temperature profiles inside the physical domain for the

CR IP T

classical Fourier as well as non-Fourier heat conduction models (e.g. Cattaneo-Vernotte and DPL). In one of the other studies, Lam and Fong [38] investigated the thermal behaviour of a solid material subjected to surface as well as internal heating by solving one-dimensional dual phase lag model-based heat conduction equation. The overall solution methodology employed the principle of superposition and the method of Fourier series expansion in

AN US

conjunction with solution structure theorems.

The studies based on analytical approaches for solving the Fourier and/or non-Fourier form of bio-heat transfer equations have primarily made use of the inversion of Laplace transform as part of the hybrid schemes. Though the inversion of Laplace transform for one-

M

dimensional cases can be achieved with the help of standard Laplace inversion tables, one has to turn to numerical approaches in the case of two-dimensional studies. The mathematical

ED

complexities associated with Laplace transform-based schemes lie in the fact that for multidimensional transient problems (e.g. a two-dimensional tissue sample), the partial derivative with respect to time can be removed from the governing partial differential equation using

PT

Laplace transformation. If the resultant partial differential equation is homogeneous in nature, one can apply methods such as separation of variables for obtaining the complete analytical

CE

solution. However, if the Laplace transformed partial differential equation is not homogenous in nature, the solution methodology becomes even more challenging. In addition, inversion

AC

from the Laplace domain into the physical domain may be mathematically challenging if the standard Laplace transform table is not available. These complexities with Laplace transformbased hybrid scheme have motivated the researchers for exploring other transform techniques, and in this context, finite integral transform (FIT) has attracted considerable attention in recent times. The inherent advantage of FIT lies in the fact that it directly converts a partial differential equation into ordinary differential equation, which can subsequently be solved without much mathematical complexity. Moreover, compared to Laplace transform, the inversion is also relatively straight forward. Hence the implementation of FIT for handling problems involving time independent/dependent boundary conditions and 7

ACCEPTED MANUSCRIPT source terms becomes relatively easier in comparison with Laplace transform-based schemes. In view of these advantages, a range of researchers have made use of FIT for obtaining analytical solutions of Fourier as well as hyperbolic heat conduction equations. Some of these notable studies include the work reported by Serfaty and Cotta [39] wherein the authors developed a hybrid numerical-analytical solution using the generalized integral transform technique for nonlinear diffusion problems with variable equation coefficient. Cotta and Milkhailov [40] employed a generalized integral transform technique to solve nonlinear heat

CR IP T

diffusion problems. Singh et al. [41] developed an exact analytical solution of Fourier based heat conduction equation for asymmetric transient heat conduction in a multilayer annulus subjected to time-dependent boundary conditions using FIT. Frankel et al. [42] presented a general one-dimensional temperature and heat flux formulation for hyperbolic heat conduction equation in composite medium using generalized FIT technique. Abdel-Hamid

AN US

[43] obtained an analytical solution of one-dimensional hyperbolic model based heat conduction equation in a finite medium subjected to periodic heat flux boundary conditions using FIT technique. Monteiro et al. [44] analysed heat transfer in a finite slab subjected to boundary conditions of prescribed heat flux and convection heat transfer by solving one-

M

dimensional hyperbolic based heat conduction equation using the generalized integral transform technique. Tang and Araki [45] derived the analytical solution of one-dimensional DPL model based heat conduction equation in finite rigid slab irradiated by short pulse lasers

ED

using Green’s function method and FIT technique. Cotta et al. [46] employed the generalized integral transform technique to solve the Pennes bio-heat transfer equation to determine

PT

temperature distribution inside the one-dimensional biological tissue with spatially varying thermo-physical properties.

CE

The above presented literature on FIT shows that though its importance has been realized by various researchers, the technique has not yet been employed for developing the

AC

exact analytical solution of two-dimensional generalized non-Fourier heat transfer equation i.e. DPL model in the context of problems wherein the conventional Fourier-based heat conduction model does not hold good. As mentioned earlier, one of these areas is bio-heat transfer. The concerned studies reported till date have only been limited to either Fourier and hyperbolic heat conduction equations or one-dimensional DPL-based heat conduction models only. With this background, the present work reports the development and implementation of a generalized analytical solution of two-dimensional DPL-based bio-heat transfer equation using FIT technique. The two-dimensional biological tissue phantom has been subjected to time independent as well as time dependent boundary conditions. The analytical solutions 8

ACCEPTED MANUSCRIPT obtained for each case have been compared with the numerical results available in the literature for the same operating parameters. The DPL-based temperature distributions within the body of the tissue phantom as obtained using FIT have been compared with those of Fourier as well as hyperbolic model based bio-heat transfer equations. The effect of phase lags terms in the form of relaxation times associated with temperature gradients and heat flux i.e. τT and τq on the resultant thermal profiles have also been presented based on the results obtained using FIT technique.

CR IP T

2. Mathematical formulation

As mentioned in the previous section, the present work is concerned with the development of FIT-based analytical solution of the two-dimensional DPL model for determining the thermal response of laser-irradiated biological tissue phantoms. The mathematical expression for the DPL model has earlier been presented as Equation (2). Under the influence of various heat

AN US

source terms (e.g. metabolic heat generation, blood perfusion, and external heat source such as laser), the local temperature distribution inside the biological tissue samples can be determined using the following bio-heat transfer equation:

(3)

M

⃗

where Tb and QM respectively denote the arterial temperature and metabolic heat generation

ED

while Qr is the external heat source.

In order to express the DPL-based bio-heat conduction equation purely in terms of

PT

temperature, Equations (2) and (3) can be combined to derive the following form of equation: (

)

[

*

]

+

CE

(4)

Under the assumptions of Tb = constant, QM = constant, and

defined as

, the

AC

governing equation (Equation (4)) in two-dimensional Cartesian coordinate system along with the initial and boundary conditions can be expressed as (

)

*

(

)

(

)+ (5)

Initial conditions: (6)(a) (6)(b) 9

ACCEPTED MANUSCRIPT Boundary conditions: (7)(a) (7)(b) (7)(c) (7)(d)

CR IP T

where A1, A2, A3, A4 and B1, B2, B3, B4 are constants. These constants have been used for converting one form of the boundary condition into another form. For example, when constants “A” become zero, the mixed boundary condition gets converted into the Dirichlet form of boundary condition (for instance, absolute values of temperature specified at the given boundary of the domain). On the other hand, as constants “B” are made equal to zero,

AN US

one can achieve the Neumann form of boundary condition from the mixed boundary condition. It is also to be noted here that the boundary conditions represented through Equation 7 above are the mixed boundary conditions based on the conventional Fourier law of heat conduction. The Fourier form of these boundary conditions have been employed primarily due to the limitation of finite integral transform technique in successfully

M

transforming the convective boundary conditions that have relaxation times parameters associated with them. In general, the non-Fourier forms of boundary condition equations

ED

contain mixed derivative term (derivative with respect to time and space as well). In view of this limitation of FIT, in the present work, the thermal analysis of the laser-irradiated

PT

biological phantoms has been performed with convective boundary conditions without considering the effects of relaxation time parameters.

CE

2.1 Finite integral transform (FIT) technique Following Hahn and Ozisik [47], the central idea behind FIT is based on the removal of space derivative from the partial differential equation using the concept of integral transformation.

AC

Subsequently, the governing partial differential equation can be reduced to ordinary differential form in time variable for the transformation of the quantity of interest e.g. temperature, as in the present case. The ordinary differential equation can be solved subjected to the transformed initial conditions and the result is inverted successively to obtain the solution for the desired temperature distribution. The integral transform pair for a given function

with respect to the x variable can be expressed as below [47]:

Inversion formula:

10

ACCEPTED MANUSCRIPT (

∑

) (

̅(

)

)

(8)(a)

Integral transform: ̅(

)

(

∫

)

(8)(b)

and the integral transform pair for the function ̅ (

) with respect to variable y can be

expressed as ̅(

̃̅ (

∑

)

CR IP T

Inversion formula:

)

Integral transform: ̃̅ (

eigenvalues (

),

integral ( ( )

̅(

∫

)

(9)(b)

AN US

The

)

(9)(a)

eigenfunctions ( (

)

) and

the

normalization

) for different combinations of boundary conditions can be found in

[47]. Using the technique of integral transform, the governing equation (Equation (5)) is first

M

transformed with respect to variable x (Equation (8)(b)) and then with respect to the variable

repeated here. ̃( ̅

)

̃( ̅

̃( ̅

)

PT

where

ED

y (Equation (9)(b)). Details of this procedure can be found in [47] and hence have not been

(

CE AC (

∫

)

∫

(

*

)

(

(10)

)

(11)(a)

)

(11)(b)

(

) (

+

)

) (

*

∫

(

)

)

(

∫ ∫

(

)*

∫

(

)

+ (

)

(

+

)

*

∫ (

)* (

)

)

+

+

∫ (12) 11

ACCEPTED MANUSCRIPT The first four terms on the right hand side of Equation (12) represent the transformed boundary conditions while the fifth term represents the transformed metabolic heat generation. The sixth term represents the transformed external heat source. If one or more of the boundaries of tissue phantom are subjected to constant temperature boundary conditions, there is a possibility of the values of A becoming zero. In order to take care of these limitations, necessary modifications have been made in Equation (12). For example, in the case of the left boundary of tissue phantom being subjected to

(

by

)

.

)

CR IP T

(

constant wall temperature boundary condition, the term

has been replaced

A look at Equation (10) reveals that it is a second order non-homogeneous differential and

are the constant parameters. The

AN US

equation with respect to only one variable t while

general solution of second order differential equation can be obtained in the following form: ̃̅ (

[

)

where

]

√( )

and

[

̃( ̅̅̅̅̅

]

)

(13)

are the characteristics roots of Equation (10), A and B are

̅̅̅̅̅( any arbitrary constants and ̃

M

) is the particular solution of Equation (10).

The initial conditions can also be transformed as )

̃( ̅

)

∫ ∫

∫

ED

̃̅ (

∫

(

(

̅̃ (

)

̃̅ (

)

) (14)(a) ) (14)(b)

PT

Before finding the values of A and B, the total time span has first been divided into N parts.

then ̃̅ (

AC

If

CE

Therefore, initial conditions (t = tn) can be given as ̃̅ (

) and

) becomes ̅̃ (

̃ ̅(

) while

) ̃ ̅(

(15) )

̃( becomes ̅

).

Equation (13) suggests that depending upon the relative magnitudes of ( ) and , one

can have the case of either overdamping [( )

] or underdamping [( )

]. Hence the

values of A and B are obtained using the initial conditions (Equation (15)) for the overdamping or underdamping cases, and then substituted into Equation (13). Once the transformed temperature distribution ( ̃̅ (

)) at time instant

12

has been

ACCEPTED MANUSCRIPT obtained, it has been inverted from transformed temperature ( ̃̅ ( (

)) to temperature

) using the inversion formula given by Equations (8)(a) and (9)(a), an exercise

which leads to the following expressions: * ̃̅ (

)

} ∑

(

∑

(

̃( ̅̅̅̅̅̅

)

(

)

∑

]

∑ ̃( ̅̅̅̅̅

∑

(

∑ (

∑

)

+[

̃( ̅

)

(16)(a)

)

)

+[

̃( ̅

] )

̃( ̅̅̅̅̅̅

)

) ( )

√

, and

+[

) ( )

(

̃( ̅̅̅̅̅

)

√( )

where

)

}

* ̃̅ (

)

*

) (

)

(

(

{

̃( ̅̅̅̅̅

(

∑

) (

) )

ED

∑

)+ *

}

(

{

) }+

{

(

∑

)+ *(

{

)

*

)

̃( ̅̅̅̅̅

)

CR IP T

) (

AN US

(

∑

M

∑

]

(16)(b)

( )

It is to be noted here that the subscript "n" denotes the previous time instant while subscript

PT

"n+1" represents the current time instant. 2.2 Solution for Fourier law-based bio-heat transfer equation

CE

As mentioned earlier, for

, the DPL model gets converted into the standard

Fourier model of heat conduction. Therefore, multiplying Equation (10) by

AC

substituting

where

The expression for

, results into the following first order differential equation ̃( ̅

(

and

(

and

)

̃ ̅(

(

)

)

)

(17)(a) (17)(b)

) is obtained by substituting

into Equation (12).

The general solution of the first order differential equation is given below ̃̅ (

)

∫(

(

13

)

)

(18)

ACCEPTED MANUSCRIPT where A is any arbitrary constant. The value of A is obtained by using the following initial ̃( condition: (t = tn), ̅

). Once the transformed temperature ( ̃̅ (

obtained by solving Equation (18), it is inverted back to temperature (

)) is ) by using

the inversion formula (Equations (8)(a) and (9)(a)) leading to the following expression: ∑

∑

(

)

(

̃̅ (

)

)

(19)

3. Problem statement

CR IP T

Application of FIT based generalized analytical solution has been demonstrated in the context of two-dimensional rectangular domain with dimension

, as schematically shown in

Figure 1. As shown, the left (x = 0) and right (x = L) walls of the rectangular domain have been maintained under constant wall temperature boundary conditions (Tb = 37oC). Three different test cases have been formulated based on the nature of boundary conditions imposed

AN US

on the bottom (y = 0) and the top wall (y = W) of the physical domain. (0, W) y

M

x

(0, 0)

(L, 0)

(a)

10

6

PT

y (mm)

8

4

12

10

8

y (mm)

ED

12

6 311

4 313

311

CE

2

315

315

2

317

317

0

200

319 321

400

319

600

800

1000

0

x (mm)

321

0

200

400

600

800

1000

x (mm)

(b)

AC

0

(c)

Figure 1: (a) Schematic of physical domain under consideration; Indicative temperature contours for Case 1 (constant surface temperature heating) at two different time instants of t = 50 s (b) and t = 100 s (c) have also been shown as representative cases.

The following forms of boundary conditions on the top and bottom walls of the physical domain have been employed: 1)

,

(Case 1: Constant surface temperature heating)

14

ACCEPTED MANUSCRIPT

2)

,

3)

(Case 2: Sinusoidal surface heating)

,

(Case 3: Laser-irradiated biological

tissue phantom) As mentioned above, for each of the three test cases, the boundary conditions on the left and right walls of the two-dimensional domain have been fixed in the following form:

CR IP T

(20)(a) (20)(b)

The derivations of the FIT based analytical solution for these three test cases have now been presented in the following sub-sections: 3.1 Test case I (Constant surface temperature heating)

AN US

As part of the first case study, one of the dimensions of the physical domain schematically shown in Figure 1 has been stretched to a fairly large value in comparison with the other dimension, resulting into L = 1 m and W = 0.01208 m. The case study with relatively higher value of aspect ratio (L/W) has been undertaken so as to compare the results obtained using the FIT-based analytical approach with those reported by Liu [23] in the context of one-

M

dimensional skin samples wherein the authors employed the hybrid scheme for solving the governing equation. The comparative study has been performed under the same operating

source (

ED

parameters. Following [23], the metabolic heat generation has been neglected, external heat ) is zero and the thermo-physical properties of the tissue and blood employed in

the analysis are as follows:

PT

ρ = 1000 kg∙m-3, k = 0.2 W∙m-1∙K-1, cv = cb = 4200 J∙kg-1∙K-1, Wb = ρbωb = 0.5 kg∙m-3∙s-1. The resultant DPL-based governing equation has now been analytically solved subjected to

CE

the following boundary and initial conditions: (21)(a)

Top horizontal surface:

(21)(b)

AC

Bottom horizontal surface:

The value of

has been fixed at 12

, as employed by Liu in [23].

Initial conditions: (22)(a) (22)(b) The eigenvalues, eigenfunctions and norm are given below [47] ; (

)

; ( ) 15

(23)(a)

ACCEPTED MANUSCRIPT

and

;

;

(23)(b)

It is to be noted here that Equation (23)(a) is applicable for all the three test problems considered in the present work because the boundary conditions on the left and right walls of the domain are same in all the cases considered. Substituting the values of C1, C2, C3 and C4 into Equation (12), the expression for ) becomes (

[

)

]

(24)

CR IP T

(

The particular solution of the second order ordinary differential equation (Equation (10)) can be shown to be ̃( ̅̅̅̅̅

)

(

)

[

]

(25)

AN US

The right hand side of Equation (25) is independent of time, therefore the derivative of particular solution with respect to time is zero. ̃( ̅̅̅̅̅̅

Substituting

̃ ̅(

̃ ̅(

),

̃( ̅̅̅̅̅

)

)

)

and

M

ED

∑

{

∑

AC ∑

∑ (

(

(

}+

CE

)

∑

PT

∑

∑

̃( ̅̅̅̅̅̅

)

) into Equation (16), the analytical form of temperature distribution

inside the tissue phantom becomes

(

̃( ̅̅̅̅̅

(26)

)

)

(

∑

[

]

*

)

[

] (

[

]

*(

)

{

)

}

( )

)

(

)

[

(27)(a)

] (

+[

)

]

(

)

( ) (27)(b)

Substituting the value of

(

) from Equation (24) into Equation (18), and integrating

the resultant form of equation would enable the calculation of the arbitrary constant A after

16

ACCEPTED MANUSCRIPT making use of the initial condition. Finally, the analytical solution for the Fourier model based bio-heat transfer can be expressed as

∑ ∑

∑

(

)

∑

[

(

)

[

] (

)

] (

(28)

)

It is to be noted that Equations (27) and (28) do not satisfy the boundary conditions. In order

CR IP T

to circumvent this problem, the closed form has been derived and produced below (the details of the derivation of closed form solution can be seen in Appendix A.1): ∑ [

(

(

)

)

] (

[

] (

)

(29)

)

AN US

∑

∑

where

3.2 Test case II (Time dependent boundary conditions: Sinusoidal surface heating) A test case wherein one of the horizontal surfaces of the two-dimensional tissue phantom has been subjected to time-dependent boundary condition has now been considered. Similar

M

heating configuration is generally employed for estimating the blood perfusion terms, as also reported by Liu [23]. For the present test problem, parameters like geometrical dimensions of

ED

the physical domain, thermo-physical properties of the tissue and blood have been kept identical to those presented in the previous section. The bottom horizontal surface of the

PT

tissue phantom has been subjected to sinusoidal surface heating while the temperature at the

is a constant term,

(30)(a) (30)(b) is the amplitude of the sinusoidal surface heating and

is the

AC

Here

CE

top wall has been set equal to Tb = 37oC i.e.

angular frequency. The steady state temperature distribution inside the tissue phantom has been used as the initial condition (Equation (31)); (31)(a) (31)(b) where and

defines the steady state temperature distribution inside the tissue phantom has been defined in the form of

. 17

ACCEPTED MANUSCRIPT The steady state temperature distribution has been obtained when the bottom wall of the tissue phantom has been subjected to a constant heat flux ( ) condition. Assume

as (32)

The governing equation (Equation (33)) and the corresponding boundary conditions (Equation (34)) for

are given below

CR IP T

(33) And

(34)(a)

(34)(b)

AN US

(34)(c)

(34)(d)

Using Equations (33) and (34), the steady state temperature distribution inside the body of the tissue phantom as obtained using the FIT techniques can be shown to be of the following form (details of the development of the solution can be seen in Appendix A.2) ∑

[

]

*

+

M

(

)

)

ED

+

(

(

)

∑

(

)

[

]

*

(35)

The governing equation (Equation (36)), initial conditions (Equation (37)) and boundary

PT

conditions (Equation (38)) for (

CE

(

are given below

)

*

)+

(

) (36)

Initial conditions:

AC

(37)(a) (37)(b)

Boundary conditions: (38)(a) (38)(b) (38)(c) (38)(d)

18

ACCEPTED MANUSCRIPT The eigenvalues, eigenfunctions and norm are [47] ;

;

(

)

(39)

Equations (36) and (5) are similar, hence the procedure for obtaining the solution of is the same as that for

, as has already been described in Section 3.1, and

hence it is not being re-produced here to avoid repeatability. Substituting the values of C1, C2, C3 and C4 into Equation (12), the expression for ) becomes (

CR IP T

(

[

)

]

[

]

(40)

The particular solution of second order ordinary differential equation (Equation (10)) is ̃( ̅̅̅̅̅

[

)

]

*,

-

)

[

]

*

AN US

Differentiating Equation (41) with respect to t, we get ̃( ̅̅̅̅̅̅

,

̃( Substitute ̅

-

̃( ̅

),

̃( and ̅̅̅̅̅

(

,

-

(

,

)

)

̃( ̅̅̅̅̅

)

-

+(42)

̃( ̅̅̅̅̅̅

)

) into Equation (16) to obtain

+(41)

)

and

M

. Once

(Equation (35)) have been obtained, the final temperature distribution ) inside the tissue phantom can be determined by substituting

ED

(

and

into Equation (32) as follows: ∑

∑

,

(

∑ )

-

(

( )

]

∑

[

]

)

)

(

(

)

[

[

]

[

* }

) }

{

*,

}

+

(

)

(43)(a)

*

+ ]

]

{

)

( )

)

∑

]

∑ + *(

{

+*

[

(

* }+

(

*

(

)

+

{

∑ +

∑ +

]

*

[ )

[

∑

(

]

∑

CE

)

AC

∑

(

[

)

PT

+

(

[

]

19

*

(

)

∑

( +*

) +[

[

]

*

ACCEPTED MANUSCRIPT ] [

∑ ∑

,

(

)

(

* [

∑

-

+

Substituting

]

(

) )

]

+*

+

*,

(

)

-

( )

)

(43)(b)

CR IP T

∑

(

(

into Equation (40), subsequently

) into Equation (18) and

integrating, the arbitrary constant A can be calculated using the initial condition. Once for Fourier model based bio-heat transfer equation is known, the final temperature distribution has been obtained by substituting the expressions of

∑

(

)

[

]

*

AN US

Equation (32) which becomes:

+

+

∑

* *

+ +

[

[

]

*,

-

M

∑

∑

]

*

ED

∑

and

+

∑

(

)

(

)

,

(

(

-

)

)

+

[

(

]

into

*

)

(44)

3.3 Laser-irradiated biological tissue phantom

PT

The problem of thermal response of laser-irradiated biological tissue phantoms using the DPL-based bio-heat transfer equation has numerically been investigated in one of our recent

CE

works [20]. In the present section, the same problem as considered in [20], has been solved using the FIT-based analytical approach. A two-dimensional tissue phantom in the form of a 2 mm has been considered as the physical domain.

AC

square enclosure with dimensions 2 mm

The tissue phantom has been irradiated with short laser pulses with a pulse width of tp = 4.6667 ps. Note that short pulse lasers with similar parameters (for instance, Nd-Yag operating at fundamental wavelength) find frequent applications in medical community at clinical levels [48-51]. The optical and thermo-physical properties of the tissue phantom have been kept identical to those employed in [20]. The choice of these values of optical and thermo-physical properties has been based on the fact that the tissue phantom considered in the present study mimics the soft biological tissues with high water content. The top wall of the tissue phantom has been subjected to the following form of mixed boundary conditions: 20

ACCEPTED MANUSCRIPT

Here

Lower Horizontal wall:

(45)(a)

Top horizontal wall:

(45)(b)

. Equation 45(b) denotes the convective boundary condition at the top horizontal

surface of the tissue phantom that is subjected to short pulse laser-irradiation. It is to be mentioned here that the boundary condition at the top horizontal wall of the domain (as indicated by Equation (45)(b)) is of Fourier form as it does not take into account

CR IP T

the relaxation times (T and q) associated with the process. This simplification in the nature of the boundary condition has been adopted taking into account the limitation of FIT-based analytical approach for such applications, as has also been mentioned earlier in Section 2. Initial conditions:

(46)(a)

AN US

(46)(b)

The eigenvalues, eigenfunctions and norm are given below [47] ,

(

and

) (

)

(47)

Substituting the values of C1, C2, C3 and C4 into Equation (12), the expression for ) becomes

M

(

[

)

ED

(

]

(48)

The particular solution of the second order ordinary differential equation (Equation (10)) is ̃( ̅̅̅̅̅

[

]

[

]

(49)

PT

)

The right hand side of Equation (49) is independent of time, therefore the derivative of

AC

CE

particular solution with respect to time is zero i.e.

̅( Substituting ̃

̃( ̅̅̅̅̅̅

),

̃( ̅

)

)

(50)

̃( ̅̅̅̅̅̅

̃( ̅̅̅̅̅̅

)

)

̃( and ̅̅̅̅̅

)

into Equation (16), the generalized DPL-based temperature distribution inside the tissue phantom at the current time instant

∑

[

becomes

] (

(

)[

]] *( 21

)

)

{

∑

[ ̃̅ (

∑ }

)

ACCEPTED MANUSCRIPT

( ∑

̃( ̅

∑

)

[

∑

*

{

) {

}+ }

{

( }

(

)[ ̃ ̅(

*

+

]

(

(51)(a)

[ ̃̅ (

∑

+[

]] * )

∑

)

( )

)

(

)

)

CR IP T

( ∑

+

] (

∑

)

( )

)

(51)(b)

where

) is the temperature distribution at the previous time

(

) from Equation (48) into Equation (18), and after

.

Substituting the value of

AN US

instant

and ̃̅ (

integration, the arbitrary constant A has been calculated using ̃̅ (

). Thus, the

]

(

ED

[

∑

M

analytical solution for Fourier model based bio-heat transfer equation becomes

)[

]]

(

∑

)

(

)

∑

[ ̃̅ (

)

(52)

PT

It is to be mentioned here that the present test case pertains to the physical situation of laserirradiated biological tissue phantoms. Hence, one needs to estimate the external heat source

CE

term (that appears in the Fourier/non-Fourier bio-heat transfer equation) based on radiative heat flux as part of the solution methodology to determine the thermal response (temperature

AC

distribution) of the tissue phantoms. In the present analysis, the light propagation through biological tissue phantom has been mathematically modeled using the transient radiative transfer equation (RTE). The transient form of RTE has been numerically solved using discrete ordinate method (DOM). Details of the solution methodology along with the boundary conditions employed for solving the governing equation (RTE) have been presented in one of our recent publications [51], and hence have not been repeated in the present manuscript. The present work couples the numerical solution of the transient form of RTE (obtained in [20]) with the FIT-based analytical solution of Fourier and non-Fourier

22

ACCEPTED MANUSCRIPT model-based bio-heat transfer equations (developed as part of the present work) for further thermal analysis. Details of the solution procedure of the coupling of transient RTE Fourier as well as non-Fourier bio-heat transfer equations may be seen elsewhere [20, 21]. Following the approach discussed in [20, 21], the temperature rise due to a single pulse (

) has first

been obtained and thereafter this temperature rise has been added into the temperature distribution already obtained for the previous time instant: )

̃̅ (

)

where (

)

∫

(

)

(53)

CR IP T

̃̅ (

∫

(54)

Trapezoidal rule has been employed for evaluating the double integral appearing on the right hand side of Equation (54).

AN US

4. Results and discussion

The present section discusses the results obtained using the generalized FIT-based analytical solutions developed for each of the three test cases considered above. The analytical solutions of the Fourier and/or non-Fourier (hyperbolic as well as DPL-based heat conduction equations) given in the previous section have first been verified against the results reported in

M

the literature [23] for constant surface temperature heating. The developed generalized analytical solution is also valid for time dependent boundary conditions, which has been

ED

demonstrated using the problem definition of Test case II (i.e. sinusoidal surface heating). Finally, the results pertaining to temperature distribution obtained for the case of laser-

PT

irradiated biological tissue phantom have been compared with those obtained using numerical approach reported in [20]. The effect of the number of eigenvalues, comparison of

CE

temperature distributions predicted using Fourier and non-Fourier (hyperbolic and DPL) heat conduction models, effect of two phase lag terms associated with heat flux (τq) and temperature gradient (τT) on temperature distribution inside the biological tissues subjected to

AC

the combinations of different boundary conditions have been investigated and reported. 4.1 Constant surface temperature heating Results of the benchmark studies for the verification of analytical solution developed

for the test case of constant surface temperature heating have been discussed. The results based on the analytical solution have first been verified against the results reported in the literature [23] for one-dimensional biological skin sample subjected to constant surface heating. It is to be mentioned that the work reported in [23] has made use of the hybrid method for thermal analysis. In the present work, four different sets of eigenvalues i.e. 23

ACCEPTED MANUSCRIPT ,

,

and

have first been considered. Values of relaxation times of

heat flux (τq) and temperature gradients (τT) have been kept as 20 s and 0 s respectively. Figure 2 shows the effect of the number of eigenvalues on the temporal profiles of temperatures at the spatial location of x = 0.5 m, y = 0.00208 m. It is to be seen from the figure that the temporal profiles of temperature are almost invariant beyond number of eigenvalues. An error analysis to quantify the effect of the number of eigenvalues on the final solution revealed that the maximum error between 20×26 and 40×52 sets of

CR IP T

eigenvalues is ≈ 11% while that between 40×52 and 60×78, the error was found to be ≈ 6%. Similarly, the maximum error between the results of 40×52 and 80×104 sets of eigenvalues is ≈ 10%. Hence

number of eigenvalues have been chosen for further analysis.

The analytical solution based temperature distribution obtained using the selected combination of eigenvalues (

) has been compared with the results reported in

present study and that reported in [23].

AN US

the literature [23]. A reasonably good agreement is to be seen between the results of the

The temporal profiles of temperatures as predicted using the developed analytical solution shown in Figure 2 reveal the presence of oscillations in the region of sharp

M

discontinuity (jump) i.e. between t = 30 s to 40 s. For better clarity, this portion of the profiles has been zoomed out and shown in the inset of Figure 2. The observed oscillations are

ED

significantly strong for the lower values of the number of eigenvalues chosen. The figure clearly demonstrates the effect of increasing number of eigenvalues (in x as well as in y directions) and it is to be seen that the oscillations tend to get suppressed at higher number of

PT

eigenvalues. For instance, the temperature profile corresponding to

shows

the reduced strength of these oscillations and is in reasonably good match with the results

CE

reported in [23]. The plausible reasons for such oscillations in the analytical solutions that are observed in the regions of sharp discontinuities have been discussed later in this section. It is

AC

also to be mentioned here that with increasing number of eigenvalues, the computational time also significantly increases. Hence the results discussed further in the present work pertain to an eigenvalues combination of

as a trade off between the accuracy of the

analytical solution and the computational time.

24

ACCEPTED MANUSCRIPT 10 20x26 40x52 60x78 80x104 Liu (2008)

8

 (K)

6

4

4

2

0

0 0

0

50

10

20

100

30

150

t (s)

40

50

200

60

70

CR IP T

2

80

250

300

Figure 2: Effect of number of eigenvalues on temporal profiles of temperature at x = 0.5 m, y = 0.00208 m.

AN US

Figure 3 shows the comparison of temporal profiles of temperatures as predicted using the Fourier (τq = τT = 0 s), Hyperbolic (τq = 20 s, τT = 0 s) and DPL (τq = 20 s, τT = 0.05 s) based heat conduction models at spatial location of x = 0.5 m, y = 0.00208 m within the body of the physical domain. It is to be seen from the figure that the Fourier model of heat

M

conduction predicts an almost instantaneous rise in temperature values immediately upon the initiation of the heating process. This instantaneous prediction of temperature changes

ED

without any time delay is to be attributed to the infinite speed of thermal wave propagation associated with the conventional Fourier-based heat conduction model. In contrast, one can clearly see a definite time lag before any significant changes in the temperature values are

PT

predicted on the basis of non-Fourier (hyperbolic and DPL) heat conduction models. In physical terms, this observation may be explained as follows: According to non-Fourier heat

CE

conduction models, the non-homogeneous nature of the biological tissue samples and the resultant micro-structural interactions lead to phase lags or thermal relaxation times

AC

associated with temperature gradients (τT) and the heat flux (τq). These phase lag terms can be interpreted as the mathematical representation of the thermal inertia of the biological medium. The possible impact of these thermal lags on the resultant temperature distribution may be observed from the temperature profiles corresponding to non-Fourier heat conduction models shown in Figure 3 wherein the corresponding temperature profiles reveal that it takes some finite amount of time for the changes in temperature to propagate to any spatial location within the body of the tissue phantom.

25

ACCEPTED MANUSCRIPT

10 Fourier Hyperbolic DPL 8

 (K)

6

4

4

0

0

20

0

50

100

30

150

t (s)

40

50

200

60

250

CR IP T

2

2

300

Figure 3: Comparison of temporal profiles of temperature for different heat conduction models at x = 0.5 m, y = 0.00208 m.

AN US

It is also to be seen from the figure that the hyperbolic heat conduction model predicts relatively higher temperature in comparison with DPL heat conduction model. This comparison has been shown in the inset of Figure 3 for better clarity. Furthermore, the temperature difference between the Fourier and non-Fourier heat conduction models after t =

M

250 s is almost negligible. Therefore, it may be concluded that the temperature distributions predicted by Fourier and non-Fourier heat conduction models differ only during the initial

ED

time instant. Afterwards, the temperature distributions as predicted by any of the three heat conduction models (Fourier and/or non-Fourier) are almost identical and the profiles nearly overlap on each other.

PT

The effect of relaxation time associated with temperature gradients i.e. τT on temporal profiles of temperature at x = 0.5 m, y = 0.00208 m has been shown in Figure 4(a). Values of

CE

τT have been varied as 0, 0.05, 0.5 and 5.0 s while τq = 20 s has been kept fixed. It is to be seen from the figure that the amplitude of oscillations in the thermal profiles keep getting

AC

suppressed with increasing values of τT. For instance, maximum oscillations in the temperature profiles are to be observed for τT = 0 s, while the profiles for higher values of τT are almost smooth. The observed trend is to be attributed to the dominance of thermal diffusion over thermal inertia at relatively higher values of τT. Furthermore, one can also see that rise in temperature takes place relatively early for higher values of τT as compared to its lower values. This is expected due to relatively higher speed of thermal wave propagation in the medium that would be realized for larger values of τT. Figure 4(b) shows the effect of heat flux based relaxation time i.e. τq on the temporal profiles of temperature at x = 0.5 m, y = 0.00208 m. Here the value of τT has been fixed at 0.05 s. It is to be seen from the figure that 26

ACCEPTED MANUSCRIPT as the values of τq increase, the start of temperature rise at the spatial location considered gets relatively delayed due to the fact that the speed of thermal wave decreases with increasing values of τq, which in turn result into the dominance of inertia effects over that of the diffusion process. Moreover, the magnitude of temperature rise is higher for τq = 30 s in comparison with that observed for τq = 10 s. 10

T=0 s T=0.05 s T=0.5 s T=5.0 s

8

8

6

6

 (K)

 (K)

q=10 s q=20 s q=30 s

CR IP T

10

4

4

4

2

2

2

0

20

0

50

25

100

30

35

40

150

45

200

50

55

60

250

t (s)

AN US

0

0

300

0

50

100

150

200

250

300

t (s)

M

(b) (a) Figure 4: (a) Effect of τT on temporal profiles of temperature for τq = 20 s (b) Effect of τq on temporal profiles of temperature (τT = 0.05 s) at x = 0.5 m, y = 0.00208 m. A closer look at the profiles of temperature variations with respect to time shown in

ED

Figures 2, 3 and 4 above reveals the presence of some oscillations that are primarily localized in the region(s) of sharp discontinuity (jumps) in the profiles shown. The observed trend is to be attributed to the fact that, ideally, the analytical solution of non-Fourier heat transfer

PT

equation is in the form of an infinite series (Fourier series). However, under practical conditions, finite number of terms of this infinite series expansion are generally considered.

CE

In the present work, the infinite series has been truncated to a finite number of terms for determining the temperature distribution. Thus, the oscillations that are to be observed in the

AC

region(s) of sharp discontinuity are primarily due to the mathematical limitation of Fourier series in successfully capturing the physical phenomena in regions close to such sharp discontinuities. As has earlier been demonstrated through Figure 2 above, the strength (magnitudes) of these oscillations may be suppressed by increasing the number of terms (eigen functions) in the Fourier series. The other factor that may lead to some oscillations in the solution is the inherent nature of the non-Fourier heat conduction models (dual phase lag as well as hyperbolic) that predicts the wave nature of the thermal front propagating through the physical medium.

27

ACCEPTED MANUSCRIPT 4.2 Time dependent boundary conditions (Sinusoidal surface heating) The present section discusses the thermal response of the physical domain when its bottom horizontal surface has been subjected to time-dependent boundary conditions (sinusoidal surface heating). The discussion has been based on the results obtained using the generalized FIT-based analytical solution as derived in the previous section. With reference to the expression of the sinusoidal form of surface heating, presented in Section 3, the values of and

,

have respectively been fixed as 1000 W/m2, 500 W/m2 and 0.02 rad/s. In order to

three different combinations of eigenvalues i.e. been considered.

,

CR IP T

investigate the sensitivity of the analytical solution with respect to the number of eigenvalues, and

have initially

The effect of number of eigenvalues on temporal profiles of temperature at the spatial

AN US

location of x = 0.5 m, y = 0.00208 m has been shown in Figure 5(a). At any given time instant, the maximum difference in the temperature values predicted using set of eigenvalues is seen to be ≈ 0.7 K, while for

and

and set of

eigenvalues, the maximum temperature difference in the predicted temperatures is ≈ 0.2 K. The observed difference in the temperature values may be attributed to the steady state

M

temperature that has been employed as the initial condition. To elaborate this point, the spatial profiles of the difference in steady state temperature values corresponding to the two

ED

sets of eigenvalues at x = 0.5 m have been plotted in Figure 5(b). As can be seen from the figure, at x = 0.5 m and y = 0.00208 m, the maximum difference (δTSS) in the steady state

for

and

combination of

and

sets of eigenvalues is ≈ 0.7 K, while δTSS ≈ 0.2 K

eigenvalues combination. Based on these observations, a

PT

temperatures for

eigenvalues has been employed for further analysis since the

CE

associated errors with these eigenvalues are relatively low and it also takes lesser

AC

computational time.

28

ACCEPTED MANUSCRIPT 352

0.4

20x26 40x52 60x78

351.5 351

TSS(20x26)-TSS(40x52) TSS(60x78)-TSS(40x52)

0.2

350.5 0

350 349.5

T (K)

TSS (K)

-0.2

349

-0.4

348.5 348

-0.6

347.5 347

-0.8

346

0

100

200

300

400

500

t (s)

-1

0

CR IP T

346.5 0.002

0.004

0.006

0.008

0.01

0.012

y (m)

AN US

(b) (a) Figure 5: (a) Effect of number of eigenvalues on temporal profiles of temperature at x = 0.5 m, y = 0.00208 m; (b) Spatial profile of steady state temperature difference between numbers of eigenvalues and 60 × 78 with respect to along y-direction and x = 0.5 m. The comparison of temporal profiles of temperature predicted using the Fourier (τq = τT = 0 s), Hyperbolic (τq = 20 s, τT = 0 s) and DPL (τq = 20 s, τT = 0.05 s) based heat conduction models at spatial location of x = 0.5 m, y = 0.00208 m has been shown in Figure 6. It is to be seen from the figure that the temperature values predicted using Fourier-based

M

bio-heat transfer equation are comparatively lower than that based on non-Fourier bio-heat transfer equations (zoomed out version has been shown in the inset of the figure for better

ED

clarity). The observed trend may be explained on the basis of the inherent assumption of infinite speed of thermal wave propagation associated with Fourier-based heat conduction model. Under such assumption, the photons propagating through the body of the tissue

PT

phantom spend little time at any given spatial location within the physical domain. Thus, the amount of thermal energy that the medium can gain within this short span of time is not

CE

sufficient to raise the tissue temperature significantly. It is also to be seen from the figure that the temperature waveforms predicted by Fourier and non-Fourier (hyperbolic and DPL)

AC

models are out of phase with respect to each other. The observed shift in the phase of the thermal waveforms can be attributed to the fact that the thermal signal as predicted by nonFourier heat conduction models takes some finite time to reach any given spatial location within the body of the physical domain. On the other hand, Fourier model predicts an almost instantaneous change in temperatures over the whole domain upon the start of the heating process. This effect can clearly be seen at initial time instants wherein the thermal wave form predicted using Fourier model shows an immediate increase in temperature values as soon as the heating process is started, whereas the profiles corresponding to non-Fourier model show almost constant temperature values till about t ≈ 40 s. As a result of this, one sees a shift in 29

ACCEPTED MANUSCRIPT the phase of the thermal waveforms corresponding to Fourier and non-Fourier heat conduction models. 354 Fourier Hyperbolic DPL

349.9

349.85

349.8

352

349.75

82

84

86

88

90

T (K)

349.7 80

348

346

0

200

400

t (s)

600

CR IP T

350

AN US

Figure 6: Comparison of temporal profiles of temperatures predicted using Fourier, hyperbolic and DPL models at x = 0.5 m, y = 0.00208 m. Figure 7(a) shows the effect of varying values of τT on the temporal profiles of temperature at spatial location of x = 0.5 m, y = 0.00208 m. The value of τq has been fixed at 20 s. As expected, when the value of τT is increased from 0 to 5 s, the maximum amplitude of

M

temperature profiles decreases with increasing values of τT due to the dominance of thermal diffusion effects over that of thermal inertia. At any given time instant, minimum amplitude

ED

of temperature values is to be seen for τT = 5.0 s. Thus, in physical terms, the thermal lag associated with the temperature gradients i.e. τT can be interpreted as the dampener that tends

PT

to suppress the amplitude of oscillations of the thermal wave propagating through the medium. The effect of τq on temporal profiles of temperature at x = 0.5 m, y = 0.00208 m has

CE

been shown in Figure 7(b). For the present analysis, the value of τT has been kept constant at 0.05 s. It is to be seen from the figure that with increasing values of τq, the amplitude of

AC

oscillations in the temperature waveforms also increases. Furthermore, the temporal profiles of temperature for different values of τq are not in phase with respect to each other owing to the fact that the thermal wave speed is lower for higher values of τq (due to the dominance of inertia effects over that of diffusion), which in turn results into the relative phase shifts, as observed from the profile shown in the figure above.

30

ACCEPTED MANUSCRIPT 352

352

T=0 s T=0.05 s T=0.5 s T=5.0 s

351.5 351

351 350.5

350

350

349.5

349.5

T (K)

T (K)

350.5

348.5

348.5

348

348

347.5

347.5

347

347

346.5

346.5 0

200

400

346

600

0

CR IP T

349

349

346

q=10 s q=20 s q=30 s

351.5

200

400

600

t (s)

t (s)

(a) (b) Figure 7: (a) Effect of τT (for τq = 20 s); (b) τq (for τT = 0.05 s) on temporal profiles of temperatures predicted at x = 0.5 m, y = 0.00208 m.

AN US

4.3 Laser-irradiated biological tissue phantom

The discussion presented in the above two sections has demonstrated the verification results of the generalized FIT-based analytical solution of DPL-based bio-heat transfer equation for two test cases of constant surface temperature and time-dependent surface heating. In the

M

present section, the analytical solution of DPL model developed in the present work has been coupled with the solution of transient RTE for determining the temperature distribution

ED

within the body of a laser-irradiated biological tissue phantom. The study is important in the context of therapeutic applications wherein a detailed understanding of temperature distribution within the body of the target medium (laser irradiated biological samples) is

PT

important for maximizing the efficiency of the technique of photo thermal therapy and to ensure minimum possible thermal damage to normal healthy cells surrounding the embedded

CE

inhomogeneity (malignant and/or benign cells). The temperature histories with respect to time, as predicted using the present approach, have been recorded at two different spatial

AC

locations (1 and 2, as indicated in Figure 1 in reference [20]). Here Location 1 corresponds to the point of laser irradiation (centre of the top surface of the sample) and Location 2 is within the body of the tissue phantom. The effect of number of eigenvalues on temporal profiles of temperature at Locations 1 and 2 has been shown in Figure 8. It is to be seen from the figure that the temperature profiles for

eigenvalues are in good agreement with those

obtained using numerical results. Therefore, the set of

eigenvalues has been used for

further thermal analysis. The figure also demonstrates the suppression of the oscillation observed in the region of sharp discontinuity with increasing number of eigenvalues employed, an observation discussed earlier in Section 4.1. 31

ACCEPTED MANUSCRIPT

420

318

20x40 40x80 60x120 Numerical

400

316

314

360

312

340

310

320

308

300

306

0

5

10

15

20

304

0

CR IP T

T (K)

T (K)

380

280

20x40 40x80 60x120 Numerical

5

10

15

20

t (s)

t (s)

AN US

(a) (b) Figure 8: Effect of number of eigenvalues on temporal profiles of temperature at Location 1 (a) and 2 (b). The analytical and numerical results for different heat conduction models (Fourier and non-Fourier) have been compared and shown in Figure 9. The temporal profiles of temperature without marker represent the results based on the FIT-based analytical solution

M

and those with marker represent the numerical results. It is to be seen from the figure that the results obtained using numerical simulations predict relatively lower values of temperatures

ED

at Locations 1 and 2 in comparison with those obtained on the basis of the analytical approach. It is also to be seen from the figure that the Fourier heat conduction model predicts relatively lower values of temperature as compared to those calculated using non-Fourier

PT

conduction models, an observation that is to be attributed to the infinite speed of propagation of thermal wave considered in the conventional Fourier model, as has also been discussed

CE

earlier in Section 4.2. Furthermore, the hyperbolic heat conduction model predicts highest temperature values (shown in the inset of Figure 9(a) and (b) for better clarity), while DPL

AC

model predicts temperatures that lie in between Fourier and hyperbolic heat conduction models. The explanation for the observed trend is that the hyperbolic heat conduction model considers the effect of the thermal relaxation time associated with heat flux (τq) only. In contrast, the dual phase lag model takes into account the coupled effects of non-zero values of relaxation times associated with the temperature gradient (τT) as well as heat flux (τq). In physical terms, the phase lag associated with temperature gradients i.e. τT tends to suppress the amplitude of the thermal wave front. Since the value of τT is zero in the case of hyperbolic heat conduction model, the resultant temperature profiles are relatively free of such dampening effects and hence the absolute values of temperatures predicted on the basis of 32

ACCEPTED MANUSCRIPT this form of non-Fourier heat conduction model are expected to be higher than that predicted using the dual phase lag model wherein these two phase lag terms are considered to be nonzero. The profiles shown in Figure 9 clearly support this observation wherein the magnitude of temperatures as predicted using DPL model lie in between those obtained using hyperbolic heat conduction model (maximum values) and Fourier-based heat conduction model (minimum values). 320

Fourier Fourier Hyperbolic Hyperbolic DPL DPL

400

380

Fourier Fourier Hyperbolic Hyperbolic DPL DPL

318

CR IP T

420

316

400

314

395

360

T (K)

T (K)

390

312

385

380

340 375

310 370 0.98

0.985

0.99

0.995

1

1.005

1.01

1.015

316

1.02

320

315.8

AN US

308

315.6

300

306

315.4

315.2 0.98

280

304

0.985

0.99

0.995

1

1.005

1.01

1.015

1.02

10 15 20 (a) (b) t (s) Figure 9: Comparison of temporal profiles of temperature predicted using Fourier and nonFourier heat conduction models at Location 1 (a) and 2 (b) (Without marker: Analytical results; With marker: numerical results). 0

5

10

15

0

5

M

t (s)

20

ED

Figure 9(b) shows the time history of temperature at Location 2, which is situated inside the body of the biological tissue phantom. The profiles corresponding to Fourier heat conduction model show a sudden drop in temperature values immediately after the laser

PT

power is switched off (t> 1.0s). On the other hand, because of the phase lag terms associated with temperature gradients and/or heat flux, the temperature profiles predicted using non-

CE

Fourier heat conduction models (hyperbolic and DPL) show nearly constant values of temperature for a relatively longer period of time (t

11 s) before the drop in temperature

AC

values starts. This, in turn, suggests the prolonged (sustained) effects of thermal energy that has been deposited at a given spatial location within the body of the tissue phantom, according to non-Fourier heat conduction models. The present observation may be physically explained by performing an energy balance in any given control volume within the body of the tissue phantom. In this way, it may be shown that the rate at which the thermal energy accumulates within the control element equals the divergence of the heat flux (for instance, for a one-dimensional case, C v

q T   x ). In the context of biological tissue samples, the t x

relaxation time associated with the heat flux i.e. q is relatively higher in comparison with the 33

ACCEPTED MANUSCRIPT relaxation time involved with the temperature gradients (T) [17, 20]. Thus, one expects that the heat flux would take some finite amount of time before it gets diffused in the surrounding region. As a result, the net heat flux remains unchanged at any given spatial location for the corresponding time duration i.e. q x x  0 . This, in turn, leads to a nearly constant temperature profile (temperature values remaining almost invariant) during this finite time duration, as has also been observed through the profiles shown in Figure 9(b) in the case of

CR IP T

non-Fourier heat conduction models. This phenomenon also explains the observation made earlier wherein the temperature values predicted at any given instant of time using nonFourier heat conduction models were seen to be relatively higher than those obtained based on the Fourier model that works under the assumption of infinite speed of thermal wave propagation in biological medium. It is also to be noted here that no such phenomenon is to

AN US

be observed for temperature profiles at Location 1 (situated on the top horizontal surface of the tissue phantom) wherein one sees an instantaneous drop in temperature values as soon as the laser power is switched off (for t> 1 s) for all the three heat conduction models (Figure 9(a)). This nearly instantaneous drop in temperature values for t > 1 s at Location 1 is to be attributed to the fact that Location 1 is situated at the top boundary (top surface subjected to

M

convective boundary conditions) of the physical domain. Hence the lagging behaviour of the thermal front is not expected to be that significant in this region compared to Location 2,

ED

which lies well within the body of the tissue phantom. In order to develop a clear understanding of the phenomena of propagation of thermal wave front within the body of the tissue phantom as predicted using Fourier and non-Fourier

PT

heat conduction models, the two-dimensional distribution of temperature fields have been plotted in the form of temperature contours. Figure 10 shows the whole-field temperature

CE

distributions at two different time instants of t = 1 s and 5 s as obtained using the FIT-based analytical solution of Fourier (a), hyperbolic (b) and DPL(c) heat conduction models. The

AC

FIT-based solution of non-Fourier heat conduction models clearly brings out the wave nature of the thermal front propagating in the physical domain as it starts from the point of laserirradiation and moves towards the boundaries of the domain. The hyperbolic heat conduction model results into the maximum amplitude of oscillations in comparison with those observed for the DPL model. On the other hand, as expected, no such oscillations are to be seen for Fourier model wherein one observes a nearly uniform and homogenous decay of temperature as one moves from the point of laser-irradiation towards the boundaries of the physical domain. This in turn implies nearly uniform diffusion of thermal energy within the body of 34

ACCEPTED MANUSCRIPT the tissue phantom. Furthermore, for any given non-Fourier heat conduction model (hyperbolic or DPL), the oscillatory nature of the propagating thermal front is more dominant at t = 5 s as compared to that seen at t = 1 s. In physical terms, these differences seen at two different time instants can be explained by the fact that the laser-irradiation is switched off at t = 1 s and hence the thermal energy that has been locally deposited at the point of incident laser beam (centre point of the top surface of the domain) now diffuses throughout the body of the tissue phantom with the passage of time, resulting into the appearance of waviness in

CR IP T

the temperature distribution throughout the entire physical domain.

t = 5s

AN US

t = 1s

M

(a) Fourier model (q= 0, T = 0)

(c) DPL model (q = 16 s, T = 0.05 s)

AC

CE

PT

ED

(b) Hyperbolic model (q = 16 s, T = 0)

Figure 10: Analytically obtained two-dimensional temperature contours as retrieved using Fourier and non-Fourier heat conduction models at two different time instants: t = 1 s (Left Column) and t = 5 s (Right Column). The contours of two-dimensional temperature distribution shown in Figure 10

correspond to the case of laser-irradiated tissue phantom (square enclosure with dimensions 2 mm

2 mm), a physical domain which is identical to that considered in one of our recent

works [20]. The profiles shown in Figure 10 have been obtained using the analytical solution of the governing equations, as developed in the present work, while similar profiles (at the same time instants of t = 1 s and 5 s) generated based on the numerical solutions (based on 35

ACCEPTED MANUSCRIPT finite volume method) of Fourier and non-Fourier heat conduction models were presented in [20]. Thus, it becomes imperative to have a direct comparison of these two approaches, especially in terms of accuracy, computational cost etc. One of the potential advantages of the FIT-based analytical approach pertains to a significant reduction in the computational cost as compared to that generally incurred with numerical simulations. For the results presented in Figure 6 of [20], the numerical simulations were performed on MATLAB® platform using a computer with configuration of Intel® Xeon® CPU [email protected] GHz (dual processor) and

CR IP T

128 GB RAM. The computational time required for the FVM-based numerical solution of the governing equations was estimated to be ≈ 10821 s with the above-mentioned configuration. On the other hand, with the same computational resource, the analytical approach could retrieve the solution of the governing equation at any given location of the physical domain in significantly lesser time span (of the order of few seconds). In addition to the differences in

AN US

the computational costs, in order to assess the accuracy of the analytical approach, a direct comparison of the time histories of temperature distributions at two different spatial locations within the physical domain, as determined through the analytical approach (present work) and that on the basis of numerical simulations [20] has been presented in Figure 9 earlier. A

M

reasonably close agreement in the two sets of the solutions (analytical and numerical) is to be seen from the figure, which in turn demonstrates the ability of the analytical approach in accurately predicting the details of the physical phenomena under study.

ED

The effect of τT on temporal profiles of temperatures at Locations 1 and 2 as predicted using the FIT-based generalized analytical solution have been presented in Figure 11(a) and

PT

(b) respectively. For comparison, the profiles obtained using numerical approach have also been shown. The relaxation time associated with the heat flux (τq) has been kept constant

CE

(=16 s). It is to be seen from the figure that the analytical solution predicts relatively higher temperature values as compared to the numerical results. It can also be seen from the figure

AC

that as the values of relaxation time associated with the temperature gradients (τT) increase, the amplitude of oscillations in the temperature profiles keeps getting suppressed. The temporal profiles of temperature for Location 2, which is situated inside the body of the tissue phantom, have been shown in Figure 11(b). It is to be seen from Figure 11(b) that the temperature profiles obtained on the basis of analytical and numerical approaches show a drop in temperature values at different time instants. The temperature values obtained using analytical solution for τT = 0 s remains almost constant up to t ≈ 13 s before a drop in temperature starts. On the other hand, the temperature profiles corresponding to the numerical

36

ACCEPTED MANUSCRIPT solution for τT = 0 s are almost flat upto t ≈ 8 s beyond which the drop in the profiles is observed. 420

322

Analytical (T=0 s) Numerical (T=0 s) Analytical (T=0.05 s) Numerical (T=0.05 s) Analytical (T=0.5 s) Numerical (T=0.5 s)

400

380

Analytical (T=0 s) Numerical (T=0 s) Analytical (T=0.05 s) Numerical (T=0.05 s) Analytical (T=0.5 s) Numerical (T=0.5 s)

320 318 316

360

312

340

310

320 308

300

280

306

0

5

10

15

20

t (s)

304

0

5

CR IP T

T (K)

T (K)

314

10

15

20

AN US

t (s)

(a) (b) Figure 11: Effect of τT on temporal profiles of temperature at Location 1 (a) and Location 2 (b) (Without marker: Analytical results; With marker: numerical results). The effect of τq on temporal profiles of temperatures predicted at Locations 1 and 2

M

for τT = 0.05 s has been shown in Figure 12. It is to be seen from the figure that analytical solution predicts relatively higher temperature values in comparison with those obtained using numerical simulations. The temperature profiles corresponding to Location 2 (shown in

ED

Figure 12(b)) for τq = 22 s are almost constant for a relatively longer period of time as compared to those obtained for the lower values of for τq. A closer examination of the results

PT

(analytical as well as numerical) presented in Figure 12(b) reveals a rise in temperature values between 15-20 seconds for q = 10 s. However, for other higher values of q, the

CE

corresponding profiles do not reflect this trend. This observation is to be attributed to the fact that the results shown in Figure 12(b) correspond to a total time duration of 20 s. Similar

AC

trends (rise in temperature towards higher instants of time) were also observed for q values greater than 10 s in the simulations that were carried out for a longer time duration. In physical terms, this trend may be associated with the wave nature of the thermal front propagating through the physical domain, as predicted using non-Fourier heat conduction models. With the passage of time, the thermal wave propagating in the physical domain reaches the domain boundaries and gets reflected again back to the medium. This in turn results in the superposition of the backward travelling waves (reflected) with the thermal waves that travel towards the boundaries of the physical domain from the point where the

37

ACCEPTED MANUSCRIPT laser energy has been deposited (point of laser irradiation). This phenomenon of superposition would lead to periodic changes in the temperature values at any given spatial location of the medium under study, as predicted through the temperature profiles shown in Figure 12(b). 420

322 Analytical (q=10 s) Numerical (q=10 s) Analytical (q=16 s) Numerical (q=16 s) Analytical (q=22 s) Numerical (q=22 s)

380

318 316

CR IP T

400

Analytical (q=10 s) Numerical (q=10 s) Analytical (q=16 s) Numerical (q=16 s) Analytical (q=22 s) Numerical (q=22 s)

320

314

T (K)

T (K)

360

312

340

310 308

320

306 300

280

0

5

10

AN US

304 15

20

t (s)

(a)

302

0

5

10

15

20

t (s)

(b)

Figure 12: Effect of τq on temporal profiles of temperature at Location 1 and 2 (Without marker: Analytical results; With marker: numerical results).

M

5. Conclusions

Finite integral transform-based generalized analytical solution of non-Fourier DPL-based

ED

heat conduction model has been presented. The applicability of the developed generalized solutions has been successfully demonstrated in the context of photo-thermal therapy wherein the interest is in understanding the thermal response of biological tissue phantoms that are

PT

subjected to short pulse laser irradiation. For generality, the two-dimensional tissue phantom has been subjected to time-dependent as well as time independent boundary conditions. In

CE

specific terms, boundary conditions in the form of constant surface temperature, sinusoidal surface heating and short pulse laser irradiation have been considered. The test case of

AC

constant surface temperature heating has been considered for verifying the methodology developed in the present work. The temperature profiles obtained using FIT-based analytical solution for this test case are in good agreement with the numerical results reported in the literature. A detailed analysis of the thermal response of biological tissue phantom subjected to short pulse laser irradiation has then been presented using the developed analytical approach. The solution of transient RTE has been coupled with the generalized non-Fourier DPL heat conduction model. In this context, the effect of thermal relaxation times associated with heat flux (τq) and temperature gradients (τT) on the resultant temperature profiles have been discussed. A comparison of the analytical solutions of the Fourier (τq = τT = 0) with 38

ACCEPTED MANUSCRIPT those of hyperbolic (τT = 0) and DPL models revealed that the temperature values predicted using DPL heat conduction equation lie between those obtained using Fourier and hyperbolic heat conduction equations. The hyperbolic heat conduction model resulted in more pronounced wave nature of the predicted temperature profiles as compared to DPL and Fourier heat conduction models. Contours of two-dimensional temperature distributions within the body of the tissue phantom showed maximum amplitude of oscillation for the solution obtained using hyperbolic model. On the other hand, no such oscillations were seen

CR IP T

in the thermal profiles obtained through the solution of conventional Fourier-based heat conduction model. The importance of the study lies in the fact that it presents the generalized analytical solutions of Fourier and non-Fourier heat conduction models for understanding the thermal response of biological tissues phantoms that are subjected to constant as well as timedependent boundary conditions and hence finds its relevance in the area of photo-thermal

AN US

therapy.

References

[1] M. Torabi, S. Saedodin, Analytical and numerical solutions of hyperbolic heat conduction

M

in cylindrical coordinates, J.Thermophy. Heat Transfer 25 (2011) 239-253. [2] Udayraj, P. Talukdar, R. Alagirusamy, A. Das, Heat transfer analysis and second degree

ED

burn prediction in human skin exposed to flame and radiant heating using dual phase lag phenomenon, Int. J. Heat Transfer 78 (2014) 1068-1079. [3] K.J. Baumeister, T.D. Hamill, Hyperbolic heat conduction equation a solution for the

PT

semi-infinite body problem, ASME J. Heat Transfer 91 (1969) 543-548. [4] M.J. Maurer, H.A. Thompson, Non-Fourier effects at high heat flux, ASME J. Heat

CE

Transfer 95 (1973) 284-286.

[5] M. Chester, Second sound in solids, Phys. Rev. 131 (1963) 2013-2015.

AC

[6] M.S. Kazimi, C.A. Erdman, On the interface temperature of two suddenly contacting materials, ASME J. Heat Transfer 97 (1975) 615-617. [7] J.R. Ho, C.P. Kuo, W.S. Jiaung, Study of heat transfer in multilayered structure within the framework of dual-phase-lag heat conduction model using lattice Boltzmann method, Int. J. Heat Mass Transfer 46 (2003) 55-69. [8] K. Kim, Z. Guo, Multi-time-scale heat transfer modeling of turbid tissues exposed to short-pulsed irradiations, Comp. Meth. Prog. Bio. 86 (2007) 112-123. [9] A. Vedavarz, S. Kumar, M.K. Moallemi, Significance of non-Fourier heat waves in conduction, ASME J. Heat Transfer 116 (1994) 221-224. 39

ACCEPTED MANUSCRIPT [10] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, Compt. Rend. 247 (1958) 431-433. [11] P. Vernotte, Some possible complications in the phenomena of thermal conduction, Compt. Rend. 252 (1961) 2190-2191. [12] D.Y. Tzou, A unified field approach for heat conduction from macro- to micro-scales, ASME J. Heat Transfer 117 (1995) 8-16. [13] J. Ghazanfarian, Z. Shomali, Investigation of dual-phase-lag- heat conduction model in

CR IP T

nanoscale metal-oxide-semiconductor filed-effect transistor, Int. J. Heat Mass Transfer 55 (2012) 6231-6237.

[14] H. Liu, M. Bussmann, J. Mostaghimi, A comparison of hyperbolic and parabolic models of phase change of pure metal, Int. J. Heat Mass Transfer 52 (2009) 1177-1184

[15] K. Ramadan, M.A. Al-Nimr, Thermal Wave Reflection and Transmission in a

AN US

Multilayer Slab with Imperfect Contact Using the Dual-Phase-Lag Model, Heat Transfer Eng. 30 (2009) 677–687.

[16] D.Y. Tzou, K.S. Chiu, Temperature-dependent thermal lagging in ultrafast laser heating, Int. J. Heat Mass Transfer 44 (2001) 1725–1734.

M

[17] J. Zhou, J.K. Chen, Y. Zhang, Dual-phase lag effects on thermal damage to biological tissues caused by laser irradiation, Comput. Biol. Med. 39 (2009) 286-293.

ED

[18] J. Zhou, Y. Zhang, J.K. Chen, An axisymmetric dual-phase-lag bioheat model for laser heating of living tissues, Int. J. Thermal Sci. 48 (2009) 1477-1485. [19] A. Narasimhan, S. Sadasivam, Non-Fourier bio heat transfer modelling of thermal

PT

damage during retinal laser irradiation, Int. J. Heat Mass Transfer 60 (2013) 591-597. [20] S. Kumar, A. Srivastava, Thermal analysis of laser-irradiated tissue phantoms using dual

CE

phase lag model coupled with transient radiative transfer equation, Int. J. Heat Mass Transfer 90 (2015) 466-479.

AC

[21] S. Patidar, S. Kumar, A. Srivastava, S. Singh, Dual phase lag model-based thermal analysis of tissue phantoms using lattice Boltzmann method, Int. J. Thermal Sci. 103 (2016) 41–56.

[22] H.T. Chen, J.Y. Lin, Numerical analysis for hyperbolic heat conduction, Int. J. Heat Mass Transfer 36 (1993) 2891-2898. [23] K.C. Liu, Thermal propagation analysis for living tissue with surface heating, Int. J. Thermal Sci. 47 (2008) 507-513.

40

ACCEPTED MANUSCRIPT [24] H.L. Lee, T.H. Lai, W.L. Chen, Y.C. Yang, An inverse hyperbolic conduction problem in estimating surface heat flux of a living skin tissue, Appl. Math. Modelling 37 (2013) 26302643. [25] K.C. Liu, P.J. Cheng, Finite propagation of heat transfer in a multilayer tissue, J. Thermophy. Heat Transfer 22 (2008) 775-782. [26] K.C. Liu, H.T. Chen, Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment, Int. J. Heat Mass Transfer 52 (2009) 1185-1192.

CR IP T

[27] K.C. Liu, J.C. Wang, Analysis of thermal damage to laser irradiated tissue based on the dual-phase-lag model, Int. J. Heat Mass Transfer 70 (2014) 621-628.

[28] J. Liu, Y.X. Zhou, Z.S. Deng, Sinusoidal heating method to noninvasively measure tissue perfusion, IEEE Trans. Biomed. Eng. 49 (2002) 867-877.

[29] T.C. Shih, P. Yuan, W.L. Lin, H.S. Kou, Analytical analysis of the Pennes bioheat

AN US

transfer equation with sinusoidal heat flux condition on the skin surface, Med. Eng. Phys 29 (2007) 946-953.

[30] 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 (1999) 420-428.

M

[31] A. Saleh, M. Al-Nimr, Variational formulation of hyperbolic heat conduction problems applying Laplace transform technique, Int. Commun. Heat Mass Transfer 35 (2008) 204-214.

ED

[32] 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 Transfer 39 (2012) 121-130.

PT

[33] H. Askarizadeh, H. Ahmadikia, Analytical analysis of the dual-phase-lag model of bioheat transfer equation during transient heating of skin tissue, Heat Mass Transfer 50

CE

(2014) 1673-1684.

[34] H. Askarizadeh, H. Ahmadikia, Analytical study on the transient heating of a two-

AC

dimensional skin tissue using parabolic and hyperbolic bioheat transfer equations, Appl. Math. Modelling 39 (2015) 3704-3720. [35] H.Z. Poor, H. Moosavi, A. Moradi, Analysis of the DPL bio-heat transfer equation with constant and time-dependent heat flux conditions on skin surface, Thermal Sci. (2014) 57-57. [36] M. Torabi, K. Zhang, Multi-dimensional dual-phase-lag heat conduction in cylindrical coordinates: analytical and numerical solutions, Int. J. Heat Mass Transfer 78 (2014) 960966. [37] T.T. Lam, A unified solution of several heat conduction models, Int. J. Heat Mass Transfer 56 (2013) 653-666. 41

ACCEPTED MANUSCRIPT [38] T.T. Lam, E. Fong, Thermal dispersion in finite medium under periodic surface disturbance using dual-phase-lag model, ASME J. Heat Transfer 138 (2016) 032401-12. [39] R. Serfaty, R.M. Cotta, Integral transform solutions of diffusion problems with nonlinear equation coefficients, Int. Commun. Heat Mass Transfer 17 (1990) 851-864. [40] R.M. Cotta, M.D. Mikhailov, Integral transform method, Appl. Math. Modelling 17 (1993) 156-161. [41] S. Singh, P.K. Jain, Rizwan-Uddin, Finite integral transform technique to solve

conditions, Nucl. Eng. Des. 241 (2011) 144-154.

CR IP T

asymmetric heat conduction in a multilayer annulus with time dependent boundary

[42] J.I. Frankel, B. Vick, M.N. Ozisik, General formulation and analysis of hyperbolic heat conduction in composite area, Int. J. Heat Mass Transfer 1987 (30) 1293-1305.

[43] B. Abdel-Hamid, Modelling non-Fourier heat conduction with periodic thermal

AN US

oscillation using the finite integral transform, Appl. Math. Modelling 23 (1999) 899-914. [44] E.R. Monteiro, E.N. Macedo, J.N.N. Quaresma, R.M. Cotta, Integral transform solution for hyperbolic heat conduction in a finite slab, Int. Commun. Heat Mass transfer 36 (2009) 297-303.

M

[45] D.W. Tang, N. Araki, Wavy, wavelike, diffusive thermal responses of finite rigid slabs to high-speed heating of laser-pulses, Int. J. Heat Mass Transfer 42 (1999) 855-860.

ED

[46] R.M. Cotta, B.P. Cotta, C.P. Naveira-Cotta, G. Cotta-Pereira, Hybrid integral transforms analysis of the bioheat equation with variable properties, Int. J. Thermal Sci. 49 (2010) 15101516.

PT

[47] D.W. Hahn, M.N. Ozisik, Heat conduction, John Wiley & Sons, Inc., New Jersey, 2012. [48] A. Katzir, Lasers and Optical Fibers in Medicine, Academic Press, California, 1993.

CE

[49] J.L. Boulnois, Photophysical processes in recent medical laser developments: a review, Lasers Med. Sci. 1 (1986) 47-66.

AC

[50] R.R. Anderson and J.A. Parish, Selective photothermolysis: precise microsurgery by selective absorption of pulsed radiation, Science 220 (1983) 524-527. [51] S. Kumar, A. Srivastava, Numerical investigation of thermal response of laser irradiated tissue phantoms embedded with optical inhomogeneities, Int. J. Heat Mass Transfer 77 (2014) 262–277.

42

ACCEPTED MANUSCRIPT Appendix A A.1 Closed form solution Fourier sine or cosine series can be written as ∑ where

(

is Fourier coefficients and (

)

(A-1)

) are the eigenfunctions.

Using the property of orthogonality to eliminate all the terms of the series except for the

(

∫ where ( )

[ (

∫

) (

)]

)

( )

(A-2)

on both sides of Equation (A-1), and using Equation

(

)

AN US

(A-2), the Fourier coefficients becomes

Substituting

{

is norm or normalization integration.

(

Multiplying by operator ∫

)

CR IP T

terms with p = q i.e.

(

∫

)

(A-3)

(Equation (A-3) into Equation(A-1) (

∑

)

M

where

)

(

(

)

ED

∫

(A-4)

(A-5)

The first term on the right hand side of Equation (27) or (28), using the definition of b

PT

(Equation (11)(b) or Equation (17)(b)) can be expressed as ∑ ]

(

)

AC

where

[

CE

∑

∑

(

(

)

)

(∑

[

] (

)

)

(A-6)

Comparing the terms in parenthesis of the right hand side of Equation (A-6) with Equation (A-4), we get

∫

(A-7)

Where ∑ Substituting the values of

(Equation(A-8)) into Equation(A-6), it becomes

43

(A-8)

ACCEPTED MANUSCRIPT

∑ ∑

[

∑

(

(

)

)

[

] (

)

] (

(A-9)

)

A.2 Steady state temperature distribution Equation (33) and (34) are rewritten here and given below; (A-10)

CR IP T

And

(A-11)(a)

(A-11)(b) (A-11)(c)

AN US

(A-11)(d)

The integral transform pair for the function [47] Inversion formula: ∑

)

(

)

̅̅̅̅(

)

(

∫

)

(A-12)

)

(A-13)

ED

̅̅̅̅(

(

M

Integral transform:

with respect to x variable is given below

Using the procedure discussed in [47], governing equation (Equation (A-10)) and boundary conditions (Equation (A-11) (c) and (d)) have been transformed into Equation (A-14) and (A-

PT

15) respectively and are given below )

AC

)̅̅̅̅(

(

CE

̅̅̅̅̅(

̅̅̅̅̅(

)

) (

∫ ̅̅̅̅(

∫

(

)

(A-14)

)

(A-15)(a)

)

(A-15)(b)

The solution of second order ordinary differential equation (Equation(A-14)) is given below ̅̅̅̅(

[

)

]

(A-16)

where The arbitrary constants A and B have been determined by using Equation (A-15). The steady state temperature distribution within the body of the tissue phantom using Equation (A-12) and (A-16) becomes

44

ACCEPTED MANUSCRIPT

∑ (

)

[

]

*

+

(

)

∑

(

)

[

]

*

(A-17)

)

AC

CE

PT

ED

M

AN US

CR IP T

+

(

45