Author’s Accepted Manuscript A numerical study on dual-phase-lag model of bioheat transfer during hyperthermia treatment P. Kumar, Dinesh Kumar, K.N. Rai
www.elsevier.com/locate/jtherbio
PII: DOI: Reference:
S0306-4565(15)00026-1 http://dx.doi.org/10.1016/j.jtherbio.2015.02.008 TB1617
To appear in: Journal of Thermal Biology Received date: 26 November 2014 Revised date: 13 February 2015 Accepted date: 13 February 2015 Cite this article as: P. Kumar, Dinesh Kumar and K.N. Rai, A numerical study on dual-phase-lag model of bio-heat transfer during hyperthermia treatment, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2015.02.008 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 galley proof before it is published in its final citable 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.
A numerical study on dual-phase-lag model of bio-heat transfer during hyperthermia treatment P.Kumar∗, Dinesh Kumar †, K.N.Rai‡
Abstract The success of hyperthermia in the treatment of cancer depends on the precise prediction and control of temperature. It was absolutely a necessity for hyperthermia treatment planning to understand the temperature distribution within living biological tissues. In this article,dualphase-lag model of bio-heat transfer has been studied using Gaussian distribution source term under most generalised boundary condition during hyperthermia treatment. An approximate analytical solution of the present problem has been done by Finite element wavelet Galerkin method which uses Legendre wavelet as a basis function. Multi-resolution analysis of Legendre wavelet in the present case localizes small scale variations of solution and fast switching of functional bases.The whole analysis is presented in dimensionless form. The dual-phase-lag model of bio-heat transfer has compared with Pennes and Thermal wave model of bio-heat transfer and it has been found that large differences in the temperature at the hyperthermia position and time to achieve the hyperthermia temperature exist,when we increase the value of τT . Particular cases when surface subjected to boundary condition of 1st, 2nd and 3rd kind are discussed in detail. The use of dual-phase-lag model of bio-heat transfer and Finite element wavelet Galerkin method as a solution method helps in precise prediction of temperature. Gaussian distribution source term helps in control of temperature during hyperthermia treatment. So,it makes this study more useful for clinical applications.
Keywords : DPL bio-heat conduction model, Finite element wavelet Galerkin method, Boundary condition, Gaussian distribution, Hyperthermia.
1
Introduction
Cancer is a life threatening disease. It is leading cause of death in developed countries and second leading cause of death in developing countries (Global Cancer Facts and Figures , 2011). Cancer is a group of disease characterized by unlimited growth and spread of abnormal cells. It is caused by both internal (hormones, inherited mutations)and external (Tobacco, infectious disease etc.) factors. The development of most cancers requires multiple steps that occur over many years, so their early detection and treatment are possible. Hyperthermia is one of the most commonly used thermal therapies for cancer treatment. It is an adjuvant therapy means it is used along with some other treatment modalities like chemotherapy, radiotherapy and surgery to enhance the effectiveness ∗ Corresponding author: Department of Mathematical Sciences, IIT-BHU, Email:[email protected], Mob.+918419079330 † DST-CIMS, Faculty of Science, BHU, Varanasi, India,Email:[email protected] ‡ Dept. of Mathematical Sciences,IIT-BHU,Varanasi,India,Email:[email protected]
1
Varanasi,
India,
of the treatment. In hyperthermia treatment,the temperature at the site of cancerous cell (tumor) is increased by using some external means (ultrasound, radio-frequency, micro-waves, infrared radiation, magnetically excitable thermoseeds, tube with hot water etc. (Xu et al., 2009) resulting in changing the physiology of diseased cells which leads to apoptosis (cell death). Depending on the degree of temperature raise,it may be classified into thermal ablation (46o C < T < 56o C), moderate hyperthermia (41o C < T < 46o C) and diathermia (T < 41o C).Further,on the basis of location of the disease, it is classified into local, regional and whole body hyperthermia. Local hyperthermia involves subjection of heat only to a small area such as tumor. Regional hyperthermia involves larger areas such as whole tissue and organ. Whole body hyperthermia is applied to treat metastatic cancerous cells when it spreads throughout the body (Kumar and Mohammad , 2011). For the success of any type of hyperthermia treatment precise prediction and control of temperature are always needed (Das et al., 2013; Gupta et al., 2013). There are a number of bio-heat transfer equations (Bhowmik et al., 2013) for living biological tissues given by Pennes (1948), Weinbaum and Jiji (1985), Nakayama and Kuwahara (2008) and many others. But due to simplicity Pennes bio-heat transfer equation is used most commonly for interpretation of thermal data. The conduction term in the Pennes bio-heat transfer equation is based on the classical Fourier’s law q(r, t) = −k∇T (r, t).
(1)
It assumes that the heat flux vector q(r, t) and temperature gradient ∇T (r, t) appear at the same instant of time i.e. thermal signal propagates with infinite speed. It means that any thermal disturbance produced at a certain location will be felt throughout the medium at that instant. In fact, heat is always found to propagate with a finite speed within living biological tissues as they have highly non-homogeneous inner structure. To solve the paradox occurred in the Pennes bio-heat equation, Thermal wave model of bio-heat transfer has been proposed based on singlephase-lagging(SPL) (Cattaneo , 1958; Vernotte , 1958) constitutive relation q(r, t + τq ) = −k∇T (r, t).
(2)
where a relaxation time τq has added to capture the micro-scale response in time. The relaxation time τq actually represents the time needed to establish the heat flux when a temperature gradient is suddenly imposed.Although the Thermal wave model of bio-heat transfer taken into account of micro-scale responses in time, it doesn’t capture the micro-scale responses in space. Due to this it produces some unusual thermal behavior. In order to consider the effect of micro-structural interactions along with the fast transient effects,a phase lag for temperature gradient τT has introduced in SPL constitutive relation (3) q(r, t + τq ) = −k∇T (r, t + τT ). According to this relation,the temperature gradient at a point r at time t + τT corresponds to the heat flux at r at time t + τq . The corresponding model is called the DPL model of bio-heat transfer (Tzou , 1996). Recently, Rukolaine (2014) has explained some unphysical effects of DPL model of heat conduction. Since Pennes bio-heat equation is a modified form of energy equation, so different numerical methods (Wang et al., 2008) are available in literature for solving them like Finite Difference Method (Shen et al., 2005; Yuan , 2009; Pletcher , 2009), Finite Difference-Decomposition Method (Gupta et al., 2013), Homotopy Perturbation Method (Gupta et al., 2010) and Finite Volume Method (Versteeg and Malalasekera, 1995). Shih et al. (2007) solved the Pennes bio-heat transfer 2
equation with sinusoidal surface heat flux on the skin as semi-infinite domain. Ahmadikia et al. (2012) did analytic solution of both parabolic and hyperbolic bio-heat transfer equation. Liu and Chen (2009, 2010) solved the DPL bio-heat model using hybrid numerical scheme. In the present study, we have to obtain the solution of the DPL model of bio-heat transfer equation under most generalised boundary conditions. Discretizing in space co-ordinate,the problem is converted in to a system of O.D.E’s with initial conditions.This system of O.D.E’s in unknown variables has been solved by Wavelet Galerkin approach. This reduces our problem into Sylvester equation.Solution of this Sylvester equation gives dimensionless temperature. The analytical solution using Laplace transform technique and approximate analytical solution obtained using FEWGM show good agreement. This proves the rationality and reliability of our solution scheme.
Nomenclature c cb k L qm qr t T Tb Wb q y τq τT
thermal wave propagation speed, m/s specific heat of blood,J/KgK thermal conductivity of tissue,W/mK length of tissue, m metabolic heat generation,W/m3 spatial heating source,W/m3 time,s temperature of tissue, K arterial temperature, K perfusion rate of blood,m3 /s/m3 heat flux,W/m2 spatial coordinate,m phase lag of heat flux , s phase lag of temperature gradient, s
Dimensionless variable and similarity criteria x Fo Foq FoT θ θb θw θf Pf Pr Pm Ki Bi a
Dimensionless space coordinate Fourier number or dimensionless time Dimensionless phase lag of heat flux Dimensionless phase lag of temperature gradient Dimensionless local tissue temperature Dimensionless arterial blood temperature Dimensionless surface temperature Dimensionless ambient temperature Dimensionless blood perfusion coefficient Dimensionless external heat source coefficient Dimensionless metabolic heat source coefficient Kirchhoff number Biot number dilation parameter
3
x∗
2
location of tumor in tissue
Mathematical formulation
During Hyperthermia treatment, the body tissue which is initially at a constant temperature T0 = 37o C is heated by some external heat source. In order to consider the micro-scale responses in both time and space DPL constitutive relation is used to derive the DPL model of bio-heat transfer. q(y, t + τq ) = −k∇T (y, t + τT ).
(4)
Where q is the heat flux, T the temperature,∇T the temperature gradient, k the thermal conductivity, τq the phase lag of heat flux and τT the phase lag of temperature gradient. τq has been introduced to take account of fast transient effects which induces the behavior of thermal wave whereas τT to represent micro-structural interactions. The heat flux precedes the temperature gradient for τq < τT and temperature gradient precedes the heat flux for τT < τq . In a local energy balance,the one dimensional energy equation of the present problem is ρc
∂q ∂T =− + ωb cb (Tb − T ) + qm + qr , ∂t ∂y
(5)
where ρ,c denotes density and specific heat of tissue,ωb ,cb ,Tb is the perfusion rate, specific heat and arterial temperature of blood. qm is the metabolic heat generation (Rai and Rai, 1999) and qr is the Gaussian distribution heating source term (Liu, 2011). Now eliminating q from Eq.(4) and Eq.(5) gives rise to ρcτq
∂2T ∂(∇2 T ) ∂T + ωb cb T = k(∇2 T + τT ) + (ωb cb Tb + qm + qr ), + (ρc + ωb cb τq ) 2 ∂t ∂t ∂t
(6)
which is DPL model of bio-heat transfer in non-dimensional form with initial condition T (y, 0) = Tb ,
∂T (y, 0) = 0, ∂t
(7)
generalised boundary condition A where in B.C. of First kind Second kind Third kind
∂T + B T = g (t). ∂y
(8)
A = 0, B = 1, g (t) = Tw , A = −k, B = 0, g (t) = qw , A = −k, B = h, g (t) = hTf ,
and symmetric condition as ∂T (0, t) = 0. ∂y
4
(9)
3
Solution of the Problem
Introducing the dimensionless variable and similarity criteria Tb − T0 Tw − T0 Tf − T0 T − T0 ; θb = ; θw = ; θf = ; T0 T0 T0 T0 wb cb kt kτq kτT y L; ; F = ; F = ; P = x = ; Fo = oq oT f L ρcL2 ρcL2 ρcL2 k L2 qm L2 qr0 L2 ∗ a1 Pm = ; Pr0 = ;x = ;a = 2 . kT0 kT0 L b
θ=
(10)
The system of Eqs.(6)-(9) reduces to the following form Foq
∂2θ ∂3θ ∂2θ 2 ∂θ 2 + (1 + F P ) + P θ = + F + Pf2 θb + Pm + Pr , oq oT f f ∂Fo2 ∂Fo ∂x2 ∂x2 ∂Fo
with initial condition θ(x, 0) = 0,
∂θ(x, 0) = 0, ∂Fo
(11)
(12)
generalised boundary condition A
∂θ(1, F0 ) + Bθ(1, Fo ) = g(Fo ), ∂x
where A = A , B = B L & g(Fo ) = the symmetric condition is
g (t)L T0
(13)
− B L and
∂θ(0, Fo ) = 0. (14) ∂x The Eqn.(11) reduces to Thermal Wave Model of bio-heat transfer by setting τT = 0 and reduces to Pennes model of bio-heat transfer by setting τq = τT . Discrtizing in space co-ordinate using central difference after simple algebraic computation the system of Eqs.(11)-(14) reduces in vector matrix form as follows Foq
d2 θ d(M1 θ) dθ + (1 + Pf2 Foq ) + Pf2 θ = M1 θ + FoT + N, dFo2 dFo dFo
(15)
with initial condition θ(0) = [0
0
...
dθ(0) = [0 dFo
0
...
0]T
(16)
and
5
0]T ,
(17)
where M1 is k × k matrix, θ and N are k × 1 matrix i.e. ⎡ −29
21h2 ⎢ h12 ⎢
⎢ ⎢ M1 = ⎢ ⎢ ⎢ ⎣
N =
3.1
0 .. .
38 21h2 −2 h2 1 h2
.. . 0 0
0 0
−9 21h2 1 h2 −2 h2
.. . 0 0
−a(x1 −x∗ )2 2θ + P Pf m + Pr0 e b
... ... ... .. . ... ...
0 0 0 .. .
0 0 0 .. .
1 h2 −9A (21A+20hB)h2
−2 h2 17A (21A+20hB)h2
θ = [θ1
...
θ2
−a(x2 −x∗ )2 2θ + P Pf m + Pr0 e b
⎤
0 0 0
+
1 h2
1 h2 13A (21A+20hB)h2
−
2 h2
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
θk ]T . ...
(18)
(19)
T dg(Fo ) −a(xk −x∗ )2 + 2θ + P 20 . + g(Fo )) ( Pf m + Pr0 e b dFo h(21A+20hB) (20)
Finite Element Wavelet Galerkin Method
To solve the O.D.E’s Eqs.(15)-(17), let us assume that d2 θ = C T Ψ(Fo ), dFo2
(21)
where C T is unknown of 2k−1 M × 2k−1 M matrix and Ψ(Fo ) is a column vector of 2k−1 M × 1, which is defined as Ψ(Fo ) =
Ψ10
Ψ11
...
Ψ1M −1
Ψ20
...
Ψ2M −1
...
Ψ2(k−1) 0
...
Ψ2(k−1) M −1
Legendre wavelet is defined on the interval (0, 1) by (Razzaghi and Yousefi, 2000)
(m + 1/2)2k/2 Pm 2k t − n ˆ , nˆ2−1 ≤ t ≤ nˆ2+1 k k Ψnm (t) = , 0 , otherwise
. (22)
(23)
where m = 0, 1, . . . , M − 1 and n = 1, 2, . . . , 2k−1 . Here Pm (Fo ) is the well known Legendre polynomials of order m which are given as follows: P0 (Fo ) = 1, P1 (Fo ) = Fo , 2m + 1 m Fo Pm (Fo ) − Pm−1 (Fo ), m = 1, 2, 3, . . . , M − 1; Pm+1 (Fo ) = m+1 m+1
(24)
where k = 1, 2, 3, . . . , n = 1, 2, ..., 2k − 1, n ˆ = (2n − 1), m is the order of Legendre polynomial and Fo is the normalized time. They are defined on the interval [0, 1]. Pm (Fo ) is denoted by Legendre polynomials of order m which are orthogonal with respect to the weight function w(Fo ) = 1, Fo ∈ [0, 1].
6
Integrating Eq.(21) with respect to Fo , we have dθ = C T P Ψ(Fo ), dFo
(25)
where P is an operational matrix of order 2k−1 M × 2k−1 M. The operational matrix of integration (Razzaghi and Yousefi, 2001) is defined as
Fo
Ψ(s)ds = P Ψ(Fo ), Fo ∈ [0, 1],
0
⎛ ⎜ ⎜ ⎜ 1 ⎜ P = k⎜ 2 ⎜ ⎜ ⎜ ⎝
L 0 0 ... ... 0 0
O L 0 ... ... 0 0
O O L ... ... 0 0
where O and L are M − 1 × M − 1 matrices given ⎛ 2 0 ⎜ 0 0 O=⎜ ⎝ ... ... 0 0 ⎛
and
1
⎜ √ −1 ⎜ 3 ⎜ ⎜ 0 ⎜ ⎜ 0 L=⎜ ⎜ . ⎜ .. ⎜ ⎜ 0 ⎜ ⎝ 0
√1 3
0
−1 √ 15
0 .. . 0 0
0
0 0
.. . 0
0 .. . 0
... ... ... ... .. . ...
0
0
...
√1 15
√1 35
0
−1 √ 35
... ... ... ... ... ... ...
O O O ... ... O L
⎞ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠
(27)
by ⎞ 0 0 ⎟ ⎟ ... ⎠ 0
... ... ... ... 0 0 0 0 .. . 0
√
(26)
−1 (2M −1)(2M −3)
√
0 0 0 0 .. .
1 (2M −1)(2M −3)
0
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
Integrating Eq.(25) with respect to Fo , we get θ(x, Fo ) = C T P 2 Ψ(Fo ). dθ and Substituting the value θ, dF o
d2 θ dFo2
(28)
in Eq.(15), we have
M1 C T X T − C T + Z = 0, 2
Foq Pf2 )P
Pf2 P 2 ]−1
(29) T
Foq Pf2 )P
where X = (P + FoT P )[Foq I + (1 + + and Z = N d [Foq I + (1 + + 2 2 −1 Pf P ] . Eq.(29) is known as Sylvester equation. Solving this equation we obtain the value of unknown vector C T . Putting the value of C T in Eq.(28) we get the required dimensionless temperature θ in biological tissue. 7
3.2
Exact solution
In order to justify the validity of the present solution scheme(FEWGM),an exact solution is required.The exact solution of the present problem is not possible,so reducing it in more simplified form by setting Foq = FoT = 0 in Eq.(11) under boundary condition of first kind with symmetry condition, the Eqs.(11)-(14) reduce to the following form ∂θ ∂2θ + Pf2 θ = + Pf2 θ + Pm + Pr , ∂Fo ∂x2
(30)
θ(x, 0) = 0,
(31)
θ(1, Fo ) = θw ,
(32)
with initial condition first kind boundary condition and symmetric condition is ∂θ(0, Fo ) = 0. (33) ∂x Applying Laplace transform technique, the solution of Eq.(30) under the initial condition Eq.(31), boundary condition Eq.(32) and symmetric condition Eq.(33) comes out to be ∞
cosh(Bs0 )x 2θw esn1Fo Bs1 cosh(Bs1 )x + cosh(Bs0 ) s1 sinh(Bs1 ) n=1 ∞ 2esn1Fo Bsn1 cosh(Bs1 )x cosh(Pf x) 2 + −(Pf θb + Pm + Pr ) 2 sinh(B Pf2 cosh(Pf ) n=1 sn1Bsn1 sn1 ) 2 1 e−Pf Fo +(Pf2 θb + Pm + Pr ) − Pf2 Pf2
θ(x, Fo ) = θw
(34)
2
where s0 = 0, s1 = −Pf2 , sn1 = −(P f 2 + (2n + 1)2 π4 ), Bs0 = P f 2 , Bs1 = 0 and Bsn1 = sn1 + Pf2 .
4
Numerical Analysis and Discussion
As in the case of moderate hyperthermia,the temperature at the site of tumor(hyperthermia position) has been kept in between 41o C − 46o C for up to 30 minutes so that tumor becomes dead,without disaffecting the neighbouring healthy tissues.The typical thermophysical properties of living biological tissues require different boundary conditions and internal heat source for hyperthermia treatment.The selected reference values including the material properties and parameters are taken as ρ = 1000kgm−3 , C = 4.18 × 103 Jkg −1 K −1 , Wb = 8kgm−3 s−1 , Cb = 3.344 × 103 Jkg −1 K −1 , k = 0.5W m−1 , qm = 10910W m−1 , T0 = 37o C, Tf = 37o C, Tw = 37o C, L = 0.05m, qw = 10W m−2 K −1 , H = 1W m−2 K −1 , a = 0.001 and x∗ = 0.025m. MATLAB 7.5 software has been used for all computational work.
8
These graphical representation for temperature distribution at the site of tumor(hyperthermia position) shows the effect of different boundary conditions. By suitably choosing the value of a and x∗ , we can produce required heat (41o C to 46o C) at the site of tumor (x = 0.5). It is independent of B. Cs. as suggested by Fig.1a. The alteration in temperature distribution at the surface is due to different type of B.Cs, because the boundary condition of first, second and third kind represents constant heating (cooling), constant heat flux and heat convection respectively. Time to achieve the hyperthermia position in our problem is 15-30 min. As suggested by Fig.1b. for this time interval it is independent of B.Cs. Differences due to different boundary condition appear after long time interval. At the tumor location, temperature profiles by different bio-heat models (Classical, Thermal wave and DPL model of bio-heat transfer) are different in Fig.2a. Under the same physical properties and parameter, differences in temperature among Classical, Thermal wave and DPL model of bio-heat transfer are of the order of 2o C, which is of great concern as far as clinical applications are concerned. The time to achieve the hyperthermia position in Classical and Thermal wave model of bio-heat transfer is almost same,but significantly more in the case of DPL model of bio-heat transfer which is shown in Fig.2b that strongly favours the theory of DPL model of bio-heat transfer. Exact solution of Pennes model under the B.C. of first kind has been obtained by using Laplace transform technique and it is compared with approximate solution obtained by FEWGM.It can be seen from Fig.3 that the result obtained by FEWGM is approximately the same as that obtained by Laplace transform technique. As suggeted by Fig.4a,keeping the value of τT constant at 20 sec and increasing the value of τq doesn’t affect the graph significantly at the location of tumor,so τq has no significant effect as far as hyperthermia psition is concerned.Similar behaviour has been observed in Fig.4b ,where we draw θ Vs.Fo for different values of τq ,so τq has no significant effect as far as time to achieve hyperthermia position is concerned. As it is clear from Figs.5a & 5b that τT has significant effect on temperature profile in living biological tissues.In Fig.5a as we increase the value of τT by keeping the value of τq constant at 28 sec,value of dimensionless temperature decreases at the hyperthermia position. It means that temperature predictied by Thermal wave model of bio-heat transfer will always be greater than that of DPL bio-heat transfer model.In the case of hyperthermia treatment where precise prediction and control of temperature are very much required , we must take care about the value of τT . It is clear from Fig.5b that the time to achieve hyperthermia position in a living biological tissue is varying with value of τT .As we increase the value of τT , it achieves hyperthermia position more quickly. It is clear from Fig.6 that when Tw is kept constant at body temperature and the value of Qro is increased gradually from 1.85×104 to 4.85×104 , the temperature at the location of tumor increased from 41o C to 46o C, which is needed for hyperthermia treatment. If Qro is kept constant at 4.85 × 104 and θw is decreased then temperature at the site of tumor remains unaffected,but it appears on the boundary,so peak temperature of hyperthermia treatment (46o C) always keeps maintained at 46o C(Fig.7). From Figs.8 & 9 it is clear that Ki and Bi do not affect the temperature at the hyperthermia
9
position in the B.C. of second kind and third kind respectively. Although it’s effect appears on the boundary.
5
Conclusion
We have presented an approximate analytical solution of the dual-phase-lag bio-heat transfer equation using the FEWGM. It is a hybrid numerical scheme based on Legendre wavelet that combines the multi-resolution and multi-scale computational property of Legendre wavelet with the specific finite element wise analysis.Multi-resolution analysis of Legendre wavelet in the case of present problem localizes small scale variations of solution and fast switching of functional bases.FEWGM has the capability to minimize error and produce higher degree of accuracy. Another results from Figs.2a & 2b is that temperature prediction at the site of hyperthermia position and time to achieve the hyperthermia position are different by different bio-heat models.DPL model of bio-heat transfer differs significantly from Thermal wave model and Pennes model of bioheat transfer. The phase lag of heat flux τq and temperature gradient τT that represents short time thermal inertia and micro-structural interactions respectively are the main characteristics of dual-phase-lag model of bio-heat transfer.In the present study,we found that τq has no significant effect but τT affects significantly. The selection of Gaussian distribution source term as spatial heating source helps to make the hyperthermia position independent of B.C.(Figs.1a & 1b) Although effect of B.C. appears on the surface(Figs.7,8 & 9) which doesn’t affect the treatment so make the treatment independent of B.C. The use of DPL model of bio-heat transfer and FEWGM as a solution method helps in precise prediction of temperature. Gaussian distribution source term helps in control of temperature during hyperthermia treatment. So,it makes this study more useful for medical doctors in the clinical field for prediction and control of temperature.The whole analysis is presented in dimensionless form.Therefore,these results provide more comprehensive and better insight into the temperature distribution inside biological tissues.
Acknowledgements The authors of this article expresses their sincere thanks to the editor and reviewers for their valuable suggestions in the improvement of this article. The first author is also thankful to the UGC, New Delhi, India for the financial support under the JRF(17-06/2012/(i)EU-V)scheme. The second and third authors are also thankful to the DST-CIMS,BHU for providing the necessary facility during the preparation of the manuscript.
10
References Ahmadikia, H., Fazlali, R., Moradi, A., 2012. Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat conditions on skin tissue. Int. Commun. in Heat and Mass Transfer. 39, 121-130. Bhowmik, A., Singh, R., Repaka, R., Mishra, S.C., 2013. Conventional and newly developed bioheat transport models in vascularized tissues: A review. Journal of Thermal Biology. 38, 107-125. Cattaneo, C., 1958. A form of heat conduction equation which eleminates the paradox of intantaeous propagation. C.R. Acad. Sci. 247, 431-433. Das K., Singh, R., Mishra, S.C., 2013. Numerical analysis for determination of the presence of a tumor and estimation of its size and location in a tissue. Journal of Thermal Biology. 38, 32-40. Global Global Cancer Facts and Figures, 2011, Secondn edition, American Cancer Society. Gupta, P.K., Singh, J., Rai, K.N., 2010. Numerical simulaton for heat transfer in tissues during thermal therapy , J. of Themal Biology. 35,295-301. Gupta, P.K., Singh, J., Rai, K.N., 2013. A numerical study on heat transfer in tissues during hyperthermia. Mathematical and Computer Modelling. 57, 1018-1037. Gupta, P.K., Singh, J., Rai, K.N., 2013. Solution of the heat transfer problem in tissues during hyperthermia by finite difference - decomposition method, J. of Appl. Mathematics and Computation. 219, 6882-6892. Kumar, S. S., Challa, Mohammad F., 2011. Magnetic nanomaterials for hyperthermia-based therapy and controlled drug delivery. Advanced Drug Delivery Reviews. 63, 789-808. Liu, K.C., Chen, H.T., 2009. Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment.International Journal of Heat and Mass Transfer. 52, 1185-1192. Liu, K.C., Chen, H.T.,2010. Investigation for the dual phase lag behavior of bio-heat transfer. International Journal of Thermal Sciences. 49, 1138-1146. Liu, K.C., 2011. Nonlinear behavior of thermal lagging in concentric living tissues with Gaussian distribution source. International Journal of Heat and Mass Transfer. 54, 2829-2836. Nakayama., , Kuwahara, F., 2008. A general bioheat transfer model based on the theory of porous media, Int.J.Heat Mass Transfer. 51, 3190-3199. Pennes, H.H., 1948. Analysis of tissue and arterial blood temperature in the resting forearm. Journal of Applied Physiology. 1, 93-122. Pletcher, R.H., 2009. Finite Difference Method, Handbook of Numerical Heat Transfer. in: W.J. Minkowycz,E.M.Sparrow,J.Y.Murthy,J.Y,Second Edition John Wiley and Sons,Inc. Hoboken,NJ,USA. Rai, K.N., Rai, S.K., 1999. Effect of metabolic heat generation and blood perfusion on the heat transfer in the tissues with a blood vessel. Heat and Mass Transfer. 35, 74-79.
11
Razzaghi, M., Yousefi, S., 2000. Legendre wavelets direct method for variational problems. Math. Comput. Simul. 53 (3), 185-192. Razzaghi, M., Yousefi, S., 2001. The Legendre wavelets operational matrix of integration. International Journal of System Sciences. 32, 495-502. Rukolaine, S.A., 2014. Unphysical effects of the dual-phase-lag model of heat conduction. International Journal of Heat and Mass Transfer. 78, 5863. Shen, W., Zhang, J., Yang, F., 2005. Modeling and numerical simulation of bioheat transfer and biomechanics in soft tissue. Math.Comput.Modelling. 41, 1251-1265. Shih, T.C., Yuan, P., Lin, W.L., Kao, H.S., 2007. Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface. Medical Engineering & Physics. 29, 946-953. Tzou, D.Y., 1996. Macro- to Microscale Heat Transfer: The Lagging Behavior. Taylor and Francis, Washington. Vernotte, P., 1958. Les paradoxes de la theorie continue de I equation de la chleur. C.R. Acad. Sci. 246, 3154-3155. Versteeg, H.K., Malalasekera,W., 1995. An Introduction to Computational Fluid Dynamics, The Finite Volume Method. Longman Group Limited, New York. Wang, L.Q., Zhou, X., Wei, X., 2008. Heat Conduction: Mathematical Models and Analytical solutions, Springer Verlag, Berlin. Weinbaum, S., Jiji, L.M., 1985. A new simplified bioheat equation for the effect of blood flow on local average tissue temperature. ASME J.Biomech.Eng. 107, 131-139. Xu, F., Lu, T.J., Seffen, K.A., Ng, E.Y.K., 2009. Mathematical modeling of skin bioheat transfer. Appl.Mech.Rev. 62, 050801-05080135. Yuan, P., 2009. Numerical analysis of an equivalent heat transfer coefficient in a porous model for simulating a biological tissue in a hyperthermia therapy. Int.J. Heat Mass Transfer. 52, 1734-1740.
12
0.3 B.C. of first kind B.C. of second kind B.C. of third kind
Dimensionless temperature (θ)
0.25 0.2 0.15 0.1 0.05 0 −0.05
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
0.9
Fig.1a : plot of θ vs. x for different kinds of B.Cs.
1
1
Dimensionless temperature (θ)
2.5 B.C. of first kind B.C. of second kind B.C. of third kind
2
1.5 2.25 1 2.2 0.5
0
2.15 0.94 0
0.2
0.4 0.6 Dimensionless time (F )
0.96
0.98
0.8
o
Fig.1b : plot of θ vs. F0 for different kinds of B.Cs.
1
1 1
Classical model Thermal wave DPL model
Dimensionless temperature (θ)
0.25
0.255 0.25
0.2
0.245 0.48
0.5
0.52
0.15
0.1
0.05
0
0
0.2
0.4 0.6 Dimensionless distance (x)
0.8
Fig.2a : plot of θ vs. x for different models.
1
1
Dimensionless temperature (θ)
2.5 Classical model Thermal wave model DPL model
2
1.5 2.3 1 2.2 2.1
0.5
0.9
0.85 0
0
0.2
0.4 0.6 Dimensionless time (F )
0.95
0.8
o
Fig.2b : plot of θ vs. F0 for for different models.
1
1
Dimensionless temperature (θ)
0.3 0.25 Exact solution FEWGM solution
0.2 0.15 0.1 0.05 0
0
0.2
0.6 0.4 Dimensionless time (F )
0.8
1
o
Fig.3 : Comprison of Exact and FEWGM solution for B.C. of first kind.
1
0.25 τ =0, τ =20
0.247
τ =14, τ =20
0.246
τ =28, τ =20
0.245
Dimensionless temperature (θ)
q q
0.2
q
T
T T
0.244 0.15 0.495
0.5
0.505
0.1
0.05
0
0
0.2
0.4 0.6 Dimensionless distance (x)
0.8
Fig.4a : plot of θ vs. x for different values of τq .
1
1
Dimensionless temperature (θ)
0.35 τq = 0,τT = 20 0.3
τq = 14,τT = 20
0.25
τq = 28,τT = 20
0.32 0.315
0.2
0.31 0.15
0.305 0.3 0.94
0.1
0.96
0.98
1
0.05 0
0
0.2
0.4
0.6
Dimensionless time (Fo)
0.8
Fig.4b : plot of θ vs. Fo for different values of τq .
1
1
Fig.5b : plot of θ vs. Fo for different values of τT .
1
Dimensionless temperature (θ)
τT=0, τq=28
0.26
τT=10, τq=28
0.25
0.25
τ =20, τ =28 T
q
0.24
0.2
0.23 0.15
0.45
0.5
0.4 0.6 Dimensionless distance (x)
0.8
0.55
0.1
0.05
0
0
0.2
Fig.5a : plot of θ vs. x for different values of τT .
1
1
0.25 4
Dimensionless temperature (θ)
Qro=1.85*10
Qro=2.85*104
0.2
4
Qro=3.85*10
Qro=4.85*104 0.15
0.1
0.05
0
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
0.9
1
Fig.6 : plot of θ vs. x for different values of Qr0 in B.C. of first kind.
1
0.25 θ = − 0.1892 w
Dimensionless temperature (θ)
0.2
θ = − 0.1081 w
0.15
θ =0 w
0.1 0.05 0 −0.05 −0.1 −0.15 −0.2
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.7 : plot of θ vs. x for B.C. of first kind.
1
0.9
1
0.3 K = 0.0270 Dimensionless temperature (θ)
i
Ki= 0.1351
0.2
K = 0.2703 i
0.1
0
−0.1
−0.2
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.8 : plot of θ vs. x for B.C. of second kind..
1
0.9
1
0.3
Dimensionless temperature (θ)
Bi= 0.1 Bi= 0.5
0.2
B = 1.0 i
0.1
0
−0.1
−0.2
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.9 : plot of θ vs. x for B.C. of third kind.
1
0.9
1
www.elsevier.com/locate/jtherbio
PII: DOI: Reference:
S0306-4565(15)00026-1 http://dx.doi.org/10.1016/j.jtherbio.2015.02.008 TB1617
To appear in: Journal of Thermal Biology Received date: 26 November 2014 Revised date: 13 February 2015 Accepted date: 13 February 2015 Cite this article as: P. Kumar, Dinesh Kumar and K.N. Rai, A numerical study on dual-phase-lag model of bio-heat transfer during hyperthermia treatment, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2015.02.008 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 galley proof before it is published in its final citable 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.
A numerical study on dual-phase-lag model of bio-heat transfer during hyperthermia treatment P.Kumar∗, Dinesh Kumar †, K.N.Rai‡
Abstract The success of hyperthermia in the treatment of cancer depends on the precise prediction and control of temperature. It was absolutely a necessity for hyperthermia treatment planning to understand the temperature distribution within living biological tissues. In this article,dualphase-lag model of bio-heat transfer has been studied using Gaussian distribution source term under most generalised boundary condition during hyperthermia treatment. An approximate analytical solution of the present problem has been done by Finite element wavelet Galerkin method which uses Legendre wavelet as a basis function. Multi-resolution analysis of Legendre wavelet in the present case localizes small scale variations of solution and fast switching of functional bases.The whole analysis is presented in dimensionless form. The dual-phase-lag model of bio-heat transfer has compared with Pennes and Thermal wave model of bio-heat transfer and it has been found that large differences in the temperature at the hyperthermia position and time to achieve the hyperthermia temperature exist,when we increase the value of τT . Particular cases when surface subjected to boundary condition of 1st, 2nd and 3rd kind are discussed in detail. The use of dual-phase-lag model of bio-heat transfer and Finite element wavelet Galerkin method as a solution method helps in precise prediction of temperature. Gaussian distribution source term helps in control of temperature during hyperthermia treatment. So,it makes this study more useful for clinical applications.
Keywords : DPL bio-heat conduction model, Finite element wavelet Galerkin method, Boundary condition, Gaussian distribution, Hyperthermia.
1
Introduction
Cancer is a life threatening disease. It is leading cause of death in developed countries and second leading cause of death in developing countries (Global Cancer Facts and Figures , 2011). Cancer is a group of disease characterized by unlimited growth and spread of abnormal cells. It is caused by both internal (hormones, inherited mutations)and external (Tobacco, infectious disease etc.) factors. The development of most cancers requires multiple steps that occur over many years, so their early detection and treatment are possible. Hyperthermia is one of the most commonly used thermal therapies for cancer treatment. It is an adjuvant therapy means it is used along with some other treatment modalities like chemotherapy, radiotherapy and surgery to enhance the effectiveness ∗ Corresponding author: Department of Mathematical Sciences, IIT-BHU, Email:[email protected], Mob.+918419079330 † DST-CIMS, Faculty of Science, BHU, Varanasi, India,Email:[email protected] ‡ Dept. of Mathematical Sciences,IIT-BHU,Varanasi,India,Email:[email protected]
1
Varanasi,
India,
of the treatment. In hyperthermia treatment,the temperature at the site of cancerous cell (tumor) is increased by using some external means (ultrasound, radio-frequency, micro-waves, infrared radiation, magnetically excitable thermoseeds, tube with hot water etc. (Xu et al., 2009) resulting in changing the physiology of diseased cells which leads to apoptosis (cell death). Depending on the degree of temperature raise,it may be classified into thermal ablation (46o C < T < 56o C), moderate hyperthermia (41o C < T < 46o C) and diathermia (T < 41o C).Further,on the basis of location of the disease, it is classified into local, regional and whole body hyperthermia. Local hyperthermia involves subjection of heat only to a small area such as tumor. Regional hyperthermia involves larger areas such as whole tissue and organ. Whole body hyperthermia is applied to treat metastatic cancerous cells when it spreads throughout the body (Kumar and Mohammad , 2011). For the success of any type of hyperthermia treatment precise prediction and control of temperature are always needed (Das et al., 2013; Gupta et al., 2013). There are a number of bio-heat transfer equations (Bhowmik et al., 2013) for living biological tissues given by Pennes (1948), Weinbaum and Jiji (1985), Nakayama and Kuwahara (2008) and many others. But due to simplicity Pennes bio-heat transfer equation is used most commonly for interpretation of thermal data. The conduction term in the Pennes bio-heat transfer equation is based on the classical Fourier’s law q(r, t) = −k∇T (r, t).
(1)
It assumes that the heat flux vector q(r, t) and temperature gradient ∇T (r, t) appear at the same instant of time i.e. thermal signal propagates with infinite speed. It means that any thermal disturbance produced at a certain location will be felt throughout the medium at that instant. In fact, heat is always found to propagate with a finite speed within living biological tissues as they have highly non-homogeneous inner structure. To solve the paradox occurred in the Pennes bio-heat equation, Thermal wave model of bio-heat transfer has been proposed based on singlephase-lagging(SPL) (Cattaneo , 1958; Vernotte , 1958) constitutive relation q(r, t + τq ) = −k∇T (r, t).
(2)
where a relaxation time τq has added to capture the micro-scale response in time. The relaxation time τq actually represents the time needed to establish the heat flux when a temperature gradient is suddenly imposed.Although the Thermal wave model of bio-heat transfer taken into account of micro-scale responses in time, it doesn’t capture the micro-scale responses in space. Due to this it produces some unusual thermal behavior. In order to consider the effect of micro-structural interactions along with the fast transient effects,a phase lag for temperature gradient τT has introduced in SPL constitutive relation (3) q(r, t + τq ) = −k∇T (r, t + τT ). According to this relation,the temperature gradient at a point r at time t + τT corresponds to the heat flux at r at time t + τq . The corresponding model is called the DPL model of bio-heat transfer (Tzou , 1996). Recently, Rukolaine (2014) has explained some unphysical effects of DPL model of heat conduction. Since Pennes bio-heat equation is a modified form of energy equation, so different numerical methods (Wang et al., 2008) are available in literature for solving them like Finite Difference Method (Shen et al., 2005; Yuan , 2009; Pletcher , 2009), Finite Difference-Decomposition Method (Gupta et al., 2013), Homotopy Perturbation Method (Gupta et al., 2010) and Finite Volume Method (Versteeg and Malalasekera, 1995). Shih et al. (2007) solved the Pennes bio-heat transfer 2
equation with sinusoidal surface heat flux on the skin as semi-infinite domain. Ahmadikia et al. (2012) did analytic solution of both parabolic and hyperbolic bio-heat transfer equation. Liu and Chen (2009, 2010) solved the DPL bio-heat model using hybrid numerical scheme. In the present study, we have to obtain the solution of the DPL model of bio-heat transfer equation under most generalised boundary conditions. Discretizing in space co-ordinate,the problem is converted in to a system of O.D.E’s with initial conditions.This system of O.D.E’s in unknown variables has been solved by Wavelet Galerkin approach. This reduces our problem into Sylvester equation.Solution of this Sylvester equation gives dimensionless temperature. The analytical solution using Laplace transform technique and approximate analytical solution obtained using FEWGM show good agreement. This proves the rationality and reliability of our solution scheme.
Nomenclature c cb k L qm qr t T Tb Wb q y τq τT
thermal wave propagation speed, m/s specific heat of blood,J/KgK thermal conductivity of tissue,W/mK length of tissue, m metabolic heat generation,W/m3 spatial heating source,W/m3 time,s temperature of tissue, K arterial temperature, K perfusion rate of blood,m3 /s/m3 heat flux,W/m2 spatial coordinate,m phase lag of heat flux , s phase lag of temperature gradient, s
Dimensionless variable and similarity criteria x Fo Foq FoT θ θb θw θf Pf Pr Pm Ki Bi a
Dimensionless space coordinate Fourier number or dimensionless time Dimensionless phase lag of heat flux Dimensionless phase lag of temperature gradient Dimensionless local tissue temperature Dimensionless arterial blood temperature Dimensionless surface temperature Dimensionless ambient temperature Dimensionless blood perfusion coefficient Dimensionless external heat source coefficient Dimensionless metabolic heat source coefficient Kirchhoff number Biot number dilation parameter
3
x∗
2
location of tumor in tissue
Mathematical formulation
During Hyperthermia treatment, the body tissue which is initially at a constant temperature T0 = 37o C is heated by some external heat source. In order to consider the micro-scale responses in both time and space DPL constitutive relation is used to derive the DPL model of bio-heat transfer. q(y, t + τq ) = −k∇T (y, t + τT ).
(4)
Where q is the heat flux, T the temperature,∇T the temperature gradient, k the thermal conductivity, τq the phase lag of heat flux and τT the phase lag of temperature gradient. τq has been introduced to take account of fast transient effects which induces the behavior of thermal wave whereas τT to represent micro-structural interactions. The heat flux precedes the temperature gradient for τq < τT and temperature gradient precedes the heat flux for τT < τq . In a local energy balance,the one dimensional energy equation of the present problem is ρc
∂q ∂T =− + ωb cb (Tb − T ) + qm + qr , ∂t ∂y
(5)
where ρ,c denotes density and specific heat of tissue,ωb ,cb ,Tb is the perfusion rate, specific heat and arterial temperature of blood. qm is the metabolic heat generation (Rai and Rai, 1999) and qr is the Gaussian distribution heating source term (Liu, 2011). Now eliminating q from Eq.(4) and Eq.(5) gives rise to ρcτq
∂2T ∂(∇2 T ) ∂T + ωb cb T = k(∇2 T + τT ) + (ωb cb Tb + qm + qr ), + (ρc + ωb cb τq ) 2 ∂t ∂t ∂t
(6)
which is DPL model of bio-heat transfer in non-dimensional form with initial condition T (y, 0) = Tb ,
∂T (y, 0) = 0, ∂t
(7)
generalised boundary condition A where in B.C. of First kind Second kind Third kind
∂T + B T = g (t). ∂y
(8)
A = 0, B = 1, g (t) = Tw , A = −k, B = 0, g (t) = qw , A = −k, B = h, g (t) = hTf ,
and symmetric condition as ∂T (0, t) = 0. ∂y
4
(9)
3
Solution of the Problem
Introducing the dimensionless variable and similarity criteria Tb − T0 Tw − T0 Tf − T0 T − T0 ; θb = ; θw = ; θf = ; T0 T0 T0 T0 wb cb kt kτq kτT y L; ; F = ; F = ; P = x = ; Fo = oq oT f L ρcL2 ρcL2 ρcL2 k L2 qm L2 qr0 L2 ∗ a1 Pm = ; Pr0 = ;x = ;a = 2 . kT0 kT0 L b
θ=
(10)
The system of Eqs.(6)-(9) reduces to the following form Foq
∂2θ ∂3θ ∂2θ 2 ∂θ 2 + (1 + F P ) + P θ = + F + Pf2 θb + Pm + Pr , oq oT f f ∂Fo2 ∂Fo ∂x2 ∂x2 ∂Fo
with initial condition θ(x, 0) = 0,
∂θ(x, 0) = 0, ∂Fo
(11)
(12)
generalised boundary condition A
∂θ(1, F0 ) + Bθ(1, Fo ) = g(Fo ), ∂x
where A = A , B = B L & g(Fo ) = the symmetric condition is
g (t)L T0
(13)
− B L and
∂θ(0, Fo ) = 0. (14) ∂x The Eqn.(11) reduces to Thermal Wave Model of bio-heat transfer by setting τT = 0 and reduces to Pennes model of bio-heat transfer by setting τq = τT . Discrtizing in space co-ordinate using central difference after simple algebraic computation the system of Eqs.(11)-(14) reduces in vector matrix form as follows Foq
d2 θ d(M1 θ) dθ + (1 + Pf2 Foq ) + Pf2 θ = M1 θ + FoT + N, dFo2 dFo dFo
(15)
with initial condition θ(0) = [0
0
...
dθ(0) = [0 dFo
0
...
0]T
(16)
and
5
0]T ,
(17)
where M1 is k × k matrix, θ and N are k × 1 matrix i.e. ⎡ −29
21h2 ⎢ h12 ⎢
⎢ ⎢ M1 = ⎢ ⎢ ⎢ ⎣
N =
3.1
0 .. .
38 21h2 −2 h2 1 h2
.. . 0 0
0 0
−9 21h2 1 h2 −2 h2
.. . 0 0
−a(x1 −x∗ )2 2θ + P Pf m + Pr0 e b
... ... ... .. . ... ...
0 0 0 .. .
0 0 0 .. .
1 h2 −9A (21A+20hB)h2
−2 h2 17A (21A+20hB)h2
θ = [θ1
...
θ2
−a(x2 −x∗ )2 2θ + P Pf m + Pr0 e b
⎤
0 0 0
+
1 h2
1 h2 13A (21A+20hB)h2
−
2 h2
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
θk ]T . ...
(18)
(19)
T dg(Fo ) −a(xk −x∗ )2 + 2θ + P 20 . + g(Fo )) ( Pf m + Pr0 e b dFo h(21A+20hB) (20)
Finite Element Wavelet Galerkin Method
To solve the O.D.E’s Eqs.(15)-(17), let us assume that d2 θ = C T Ψ(Fo ), dFo2
(21)
where C T is unknown of 2k−1 M × 2k−1 M matrix and Ψ(Fo ) is a column vector of 2k−1 M × 1, which is defined as Ψ(Fo ) =
Ψ10
Ψ11
...
Ψ1M −1
Ψ20
...
Ψ2M −1
...
Ψ2(k−1) 0
...
Ψ2(k−1) M −1
Legendre wavelet is defined on the interval (0, 1) by (Razzaghi and Yousefi, 2000)
(m + 1/2)2k/2 Pm 2k t − n ˆ , nˆ2−1 ≤ t ≤ nˆ2+1 k k Ψnm (t) = , 0 , otherwise
. (22)
(23)
where m = 0, 1, . . . , M − 1 and n = 1, 2, . . . , 2k−1 . Here Pm (Fo ) is the well known Legendre polynomials of order m which are given as follows: P0 (Fo ) = 1, P1 (Fo ) = Fo , 2m + 1 m Fo Pm (Fo ) − Pm−1 (Fo ), m = 1, 2, 3, . . . , M − 1; Pm+1 (Fo ) = m+1 m+1
(24)
where k = 1, 2, 3, . . . , n = 1, 2, ..., 2k − 1, n ˆ = (2n − 1), m is the order of Legendre polynomial and Fo is the normalized time. They are defined on the interval [0, 1]. Pm (Fo ) is denoted by Legendre polynomials of order m which are orthogonal with respect to the weight function w(Fo ) = 1, Fo ∈ [0, 1].
6
Integrating Eq.(21) with respect to Fo , we have dθ = C T P Ψ(Fo ), dFo
(25)
where P is an operational matrix of order 2k−1 M × 2k−1 M. The operational matrix of integration (Razzaghi and Yousefi, 2001) is defined as
Fo
Ψ(s)ds = P Ψ(Fo ), Fo ∈ [0, 1],
0
⎛ ⎜ ⎜ ⎜ 1 ⎜ P = k⎜ 2 ⎜ ⎜ ⎜ ⎝
L 0 0 ... ... 0 0
O L 0 ... ... 0 0
O O L ... ... 0 0
where O and L are M − 1 × M − 1 matrices given ⎛ 2 0 ⎜ 0 0 O=⎜ ⎝ ... ... 0 0 ⎛
and
1
⎜ √ −1 ⎜ 3 ⎜ ⎜ 0 ⎜ ⎜ 0 L=⎜ ⎜ . ⎜ .. ⎜ ⎜ 0 ⎜ ⎝ 0
√1 3
0
−1 √ 15
0 .. . 0 0
0
0 0
.. . 0
0 .. . 0
... ... ... ... .. . ...
0
0
...
√1 15
√1 35
0
−1 √ 35
... ... ... ... ... ... ...
O O O ... ... O L
⎞ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠
(27)
by ⎞ 0 0 ⎟ ⎟ ... ⎠ 0
... ... ... ... 0 0 0 0 .. . 0
√
(26)
−1 (2M −1)(2M −3)
√
0 0 0 0 .. .
1 (2M −1)(2M −3)
0
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
Integrating Eq.(25) with respect to Fo , we get θ(x, Fo ) = C T P 2 Ψ(Fo ). dθ and Substituting the value θ, dF o
d2 θ dFo2
(28)
in Eq.(15), we have
M1 C T X T − C T + Z = 0, 2
Foq Pf2 )P
Pf2 P 2 ]−1
(29) T
Foq Pf2 )P
where X = (P + FoT P )[Foq I + (1 + + and Z = N d [Foq I + (1 + + 2 2 −1 Pf P ] . Eq.(29) is known as Sylvester equation. Solving this equation we obtain the value of unknown vector C T . Putting the value of C T in Eq.(28) we get the required dimensionless temperature θ in biological tissue. 7
3.2
Exact solution
In order to justify the validity of the present solution scheme(FEWGM),an exact solution is required.The exact solution of the present problem is not possible,so reducing it in more simplified form by setting Foq = FoT = 0 in Eq.(11) under boundary condition of first kind with symmetry condition, the Eqs.(11)-(14) reduce to the following form ∂θ ∂2θ + Pf2 θ = + Pf2 θ + Pm + Pr , ∂Fo ∂x2
(30)
θ(x, 0) = 0,
(31)
θ(1, Fo ) = θw ,
(32)
with initial condition first kind boundary condition and symmetric condition is ∂θ(0, Fo ) = 0. (33) ∂x Applying Laplace transform technique, the solution of Eq.(30) under the initial condition Eq.(31), boundary condition Eq.(32) and symmetric condition Eq.(33) comes out to be ∞
cosh(Bs0 )x 2θw esn1Fo Bs1 cosh(Bs1 )x + cosh(Bs0 ) s1 sinh(Bs1 ) n=1 ∞ 2esn1Fo Bsn1 cosh(Bs1 )x cosh(Pf x) 2 + −(Pf θb + Pm + Pr ) 2 sinh(B Pf2 cosh(Pf ) n=1 sn1Bsn1 sn1 ) 2 1 e−Pf Fo +(Pf2 θb + Pm + Pr ) − Pf2 Pf2
θ(x, Fo ) = θw
(34)
2
where s0 = 0, s1 = −Pf2 , sn1 = −(P f 2 + (2n + 1)2 π4 ), Bs0 = P f 2 , Bs1 = 0 and Bsn1 = sn1 + Pf2 .
4
Numerical Analysis and Discussion
As in the case of moderate hyperthermia,the temperature at the site of tumor(hyperthermia position) has been kept in between 41o C − 46o C for up to 30 minutes so that tumor becomes dead,without disaffecting the neighbouring healthy tissues.The typical thermophysical properties of living biological tissues require different boundary conditions and internal heat source for hyperthermia treatment.The selected reference values including the material properties and parameters are taken as ρ = 1000kgm−3 , C = 4.18 × 103 Jkg −1 K −1 , Wb = 8kgm−3 s−1 , Cb = 3.344 × 103 Jkg −1 K −1 , k = 0.5W m−1 , qm = 10910W m−1 , T0 = 37o C, Tf = 37o C, Tw = 37o C, L = 0.05m, qw = 10W m−2 K −1 , H = 1W m−2 K −1 , a = 0.001 and x∗ = 0.025m. MATLAB 7.5 software has been used for all computational work.
8
These graphical representation for temperature distribution at the site of tumor(hyperthermia position) shows the effect of different boundary conditions. By suitably choosing the value of a and x∗ , we can produce required heat (41o C to 46o C) at the site of tumor (x = 0.5). It is independent of B. Cs. as suggested by Fig.1a. The alteration in temperature distribution at the surface is due to different type of B.Cs, because the boundary condition of first, second and third kind represents constant heating (cooling), constant heat flux and heat convection respectively. Time to achieve the hyperthermia position in our problem is 15-30 min. As suggested by Fig.1b. for this time interval it is independent of B.Cs. Differences due to different boundary condition appear after long time interval. At the tumor location, temperature profiles by different bio-heat models (Classical, Thermal wave and DPL model of bio-heat transfer) are different in Fig.2a. Under the same physical properties and parameter, differences in temperature among Classical, Thermal wave and DPL model of bio-heat transfer are of the order of 2o C, which is of great concern as far as clinical applications are concerned. The time to achieve the hyperthermia position in Classical and Thermal wave model of bio-heat transfer is almost same,but significantly more in the case of DPL model of bio-heat transfer which is shown in Fig.2b that strongly favours the theory of DPL model of bio-heat transfer. Exact solution of Pennes model under the B.C. of first kind has been obtained by using Laplace transform technique and it is compared with approximate solution obtained by FEWGM.It can be seen from Fig.3 that the result obtained by FEWGM is approximately the same as that obtained by Laplace transform technique. As suggeted by Fig.4a,keeping the value of τT constant at 20 sec and increasing the value of τq doesn’t affect the graph significantly at the location of tumor,so τq has no significant effect as far as hyperthermia psition is concerned.Similar behaviour has been observed in Fig.4b ,where we draw θ Vs.Fo for different values of τq ,so τq has no significant effect as far as time to achieve hyperthermia position is concerned. As it is clear from Figs.5a & 5b that τT has significant effect on temperature profile in living biological tissues.In Fig.5a as we increase the value of τT by keeping the value of τq constant at 28 sec,value of dimensionless temperature decreases at the hyperthermia position. It means that temperature predictied by Thermal wave model of bio-heat transfer will always be greater than that of DPL bio-heat transfer model.In the case of hyperthermia treatment where precise prediction and control of temperature are very much required , we must take care about the value of τT . It is clear from Fig.5b that the time to achieve hyperthermia position in a living biological tissue is varying with value of τT .As we increase the value of τT , it achieves hyperthermia position more quickly. It is clear from Fig.6 that when Tw is kept constant at body temperature and the value of Qro is increased gradually from 1.85×104 to 4.85×104 , the temperature at the location of tumor increased from 41o C to 46o C, which is needed for hyperthermia treatment. If Qro is kept constant at 4.85 × 104 and θw is decreased then temperature at the site of tumor remains unaffected,but it appears on the boundary,so peak temperature of hyperthermia treatment (46o C) always keeps maintained at 46o C(Fig.7). From Figs.8 & 9 it is clear that Ki and Bi do not affect the temperature at the hyperthermia
9
position in the B.C. of second kind and third kind respectively. Although it’s effect appears on the boundary.
5
Conclusion
We have presented an approximate analytical solution of the dual-phase-lag bio-heat transfer equation using the FEWGM. It is a hybrid numerical scheme based on Legendre wavelet that combines the multi-resolution and multi-scale computational property of Legendre wavelet with the specific finite element wise analysis.Multi-resolution analysis of Legendre wavelet in the case of present problem localizes small scale variations of solution and fast switching of functional bases.FEWGM has the capability to minimize error and produce higher degree of accuracy. Another results from Figs.2a & 2b is that temperature prediction at the site of hyperthermia position and time to achieve the hyperthermia position are different by different bio-heat models.DPL model of bio-heat transfer differs significantly from Thermal wave model and Pennes model of bioheat transfer. The phase lag of heat flux τq and temperature gradient τT that represents short time thermal inertia and micro-structural interactions respectively are the main characteristics of dual-phase-lag model of bio-heat transfer.In the present study,we found that τq has no significant effect but τT affects significantly. The selection of Gaussian distribution source term as spatial heating source helps to make the hyperthermia position independent of B.C.(Figs.1a & 1b) Although effect of B.C. appears on the surface(Figs.7,8 & 9) which doesn’t affect the treatment so make the treatment independent of B.C. The use of DPL model of bio-heat transfer and FEWGM as a solution method helps in precise prediction of temperature. Gaussian distribution source term helps in control of temperature during hyperthermia treatment. So,it makes this study more useful for medical doctors in the clinical field for prediction and control of temperature.The whole analysis is presented in dimensionless form.Therefore,these results provide more comprehensive and better insight into the temperature distribution inside biological tissues.
Acknowledgements The authors of this article expresses their sincere thanks to the editor and reviewers for their valuable suggestions in the improvement of this article. The first author is also thankful to the UGC, New Delhi, India for the financial support under the JRF(17-06/2012/(i)EU-V)scheme. The second and third authors are also thankful to the DST-CIMS,BHU for providing the necessary facility during the preparation of the manuscript.
10
References Ahmadikia, H., Fazlali, R., Moradi, A., 2012. Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat conditions on skin tissue. Int. Commun. in Heat and Mass Transfer. 39, 121-130. Bhowmik, A., Singh, R., Repaka, R., Mishra, S.C., 2013. Conventional and newly developed bioheat transport models in vascularized tissues: A review. Journal of Thermal Biology. 38, 107-125. Cattaneo, C., 1958. A form of heat conduction equation which eleminates the paradox of intantaeous propagation. C.R. Acad. Sci. 247, 431-433. Das K., Singh, R., Mishra, S.C., 2013. Numerical analysis for determination of the presence of a tumor and estimation of its size and location in a tissue. Journal of Thermal Biology. 38, 32-40. Global Global Cancer Facts and Figures, 2011, Secondn edition, American Cancer Society. Gupta, P.K., Singh, J., Rai, K.N., 2010. Numerical simulaton for heat transfer in tissues during thermal therapy , J. of Themal Biology. 35,295-301. Gupta, P.K., Singh, J., Rai, K.N., 2013. A numerical study on heat transfer in tissues during hyperthermia. Mathematical and Computer Modelling. 57, 1018-1037. Gupta, P.K., Singh, J., Rai, K.N., 2013. Solution of the heat transfer problem in tissues during hyperthermia by finite difference - decomposition method, J. of Appl. Mathematics and Computation. 219, 6882-6892. Kumar, S. S., Challa, Mohammad F., 2011. Magnetic nanomaterials for hyperthermia-based therapy and controlled drug delivery. Advanced Drug Delivery Reviews. 63, 789-808. Liu, K.C., Chen, H.T., 2009. Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment.International Journal of Heat and Mass Transfer. 52, 1185-1192. Liu, K.C., Chen, H.T.,2010. Investigation for the dual phase lag behavior of bio-heat transfer. International Journal of Thermal Sciences. 49, 1138-1146. Liu, K.C., 2011. Nonlinear behavior of thermal lagging in concentric living tissues with Gaussian distribution source. International Journal of Heat and Mass Transfer. 54, 2829-2836. Nakayama., , Kuwahara, F., 2008. A general bioheat transfer model based on the theory of porous media, Int.J.Heat Mass Transfer. 51, 3190-3199. Pennes, H.H., 1948. Analysis of tissue and arterial blood temperature in the resting forearm. Journal of Applied Physiology. 1, 93-122. Pletcher, R.H., 2009. Finite Difference Method, Handbook of Numerical Heat Transfer. in: W.J. Minkowycz,E.M.Sparrow,J.Y.Murthy,J.Y,Second Edition John Wiley and Sons,Inc. Hoboken,NJ,USA. Rai, K.N., Rai, S.K., 1999. Effect of metabolic heat generation and blood perfusion on the heat transfer in the tissues with a blood vessel. Heat and Mass Transfer. 35, 74-79.
11
Razzaghi, M., Yousefi, S., 2000. Legendre wavelets direct method for variational problems. Math. Comput. Simul. 53 (3), 185-192. Razzaghi, M., Yousefi, S., 2001. The Legendre wavelets operational matrix of integration. International Journal of System Sciences. 32, 495-502. Rukolaine, S.A., 2014. Unphysical effects of the dual-phase-lag model of heat conduction. International Journal of Heat and Mass Transfer. 78, 5863. Shen, W., Zhang, J., Yang, F., 2005. Modeling and numerical simulation of bioheat transfer and biomechanics in soft tissue. Math.Comput.Modelling. 41, 1251-1265. Shih, T.C., Yuan, P., Lin, W.L., Kao, H.S., 2007. Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface. Medical Engineering & Physics. 29, 946-953. Tzou, D.Y., 1996. Macro- to Microscale Heat Transfer: The Lagging Behavior. Taylor and Francis, Washington. Vernotte, P., 1958. Les paradoxes de la theorie continue de I equation de la chleur. C.R. Acad. Sci. 246, 3154-3155. Versteeg, H.K., Malalasekera,W., 1995. An Introduction to Computational Fluid Dynamics, The Finite Volume Method. Longman Group Limited, New York. Wang, L.Q., Zhou, X., Wei, X., 2008. Heat Conduction: Mathematical Models and Analytical solutions, Springer Verlag, Berlin. Weinbaum, S., Jiji, L.M., 1985. A new simplified bioheat equation for the effect of blood flow on local average tissue temperature. ASME J.Biomech.Eng. 107, 131-139. Xu, F., Lu, T.J., Seffen, K.A., Ng, E.Y.K., 2009. Mathematical modeling of skin bioheat transfer. Appl.Mech.Rev. 62, 050801-05080135. Yuan, P., 2009. Numerical analysis of an equivalent heat transfer coefficient in a porous model for simulating a biological tissue in a hyperthermia therapy. Int.J. Heat Mass Transfer. 52, 1734-1740.
12
0.3 B.C. of first kind B.C. of second kind B.C. of third kind
Dimensionless temperature (θ)
0.25 0.2 0.15 0.1 0.05 0 −0.05
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
0.9
Fig.1a : plot of θ vs. x for different kinds of B.Cs.
1
1
Dimensionless temperature (θ)
2.5 B.C. of first kind B.C. of second kind B.C. of third kind
2
1.5 2.25 1 2.2 0.5
0
2.15 0.94 0
0.2
0.4 0.6 Dimensionless time (F )
0.96
0.98
0.8
o
Fig.1b : plot of θ vs. F0 for different kinds of B.Cs.
1
1 1
Classical model Thermal wave DPL model
Dimensionless temperature (θ)
0.25
0.255 0.25
0.2
0.245 0.48
0.5
0.52
0.15
0.1
0.05
0
0
0.2
0.4 0.6 Dimensionless distance (x)
0.8
Fig.2a : plot of θ vs. x for different models.
1
1
Dimensionless temperature (θ)
2.5 Classical model Thermal wave model DPL model
2
1.5 2.3 1 2.2 2.1
0.5
0.9
0.85 0
0
0.2
0.4 0.6 Dimensionless time (F )
0.95
0.8
o
Fig.2b : plot of θ vs. F0 for for different models.
1
1
Dimensionless temperature (θ)
0.3 0.25 Exact solution FEWGM solution
0.2 0.15 0.1 0.05 0
0
0.2
0.6 0.4 Dimensionless time (F )
0.8
1
o
Fig.3 : Comprison of Exact and FEWGM solution for B.C. of first kind.
1
0.25 τ =0, τ =20
0.247
τ =14, τ =20
0.246
τ =28, τ =20
0.245
Dimensionless temperature (θ)
q q
0.2
q
T
T T
0.244 0.15 0.495
0.5
0.505
0.1
0.05
0
0
0.2
0.4 0.6 Dimensionless distance (x)
0.8
Fig.4a : plot of θ vs. x for different values of τq .
1
1
Dimensionless temperature (θ)
0.35 τq = 0,τT = 20 0.3
τq = 14,τT = 20
0.25
τq = 28,τT = 20
0.32 0.315
0.2
0.31 0.15
0.305 0.3 0.94
0.1
0.96
0.98
1
0.05 0
0
0.2
0.4
0.6
Dimensionless time (Fo)
0.8
Fig.4b : plot of θ vs. Fo for different values of τq .
1
1
Fig.5b : plot of θ vs. Fo for different values of τT .
1
Dimensionless temperature (θ)
τT=0, τq=28
0.26
τT=10, τq=28
0.25
0.25
τ =20, τ =28 T
q
0.24
0.2
0.23 0.15
0.45
0.5
0.4 0.6 Dimensionless distance (x)
0.8
0.55
0.1
0.05
0
0
0.2
Fig.5a : plot of θ vs. x for different values of τT .
1
1
0.25 4
Dimensionless temperature (θ)
Qro=1.85*10
Qro=2.85*104
0.2
4
Qro=3.85*10
Qro=4.85*104 0.15
0.1
0.05
0
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
0.9
1
Fig.6 : plot of θ vs. x for different values of Qr0 in B.C. of first kind.
1
0.25 θ = − 0.1892 w
Dimensionless temperature (θ)
0.2
θ = − 0.1081 w
0.15
θ =0 w
0.1 0.05 0 −0.05 −0.1 −0.15 −0.2
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.7 : plot of θ vs. x for B.C. of first kind.
1
0.9
1
0.3 K = 0.0270 Dimensionless temperature (θ)
i
Ki= 0.1351
0.2
K = 0.2703 i
0.1
0
−0.1
−0.2
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.8 : plot of θ vs. x for B.C. of second kind..
1
0.9
1
0.3
Dimensionless temperature (θ)
Bi= 0.1 Bi= 0.5
0.2
B = 1.0 i
0.1
0
−0.1
−0.2
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 Dimensionless distance (x)
0.8
Fig.9 : plot of θ vs. x for B.C. of third kind.
1
0.9
1
