A study of cryosurgery of lung cancer using Modified Legendre wavelet Galerkin method

Author’s Accepted Manuscript A study of cryosurgery of lung cancer using Modified Legendre wavelet Galerkin method

Mukesh Kumar, Subrahamanyam Upadhyay, K.N. Rai www.elsevier.com/locate/jtherbio

PII: DOI: Reference:

S0306-4565(18)30284-5 https://doi.org/10.1016/j.jtherbio.2018.10.012 TB2197

To appear in: Journal of Thermal Biology Received date: 29 June 2018 Revised date: 1 October 2018 Accepted date: 13 October 2018 Cite this article as: Mukesh Kumar, Subrahamanyam Upadhyay and K.N. Rai, A study of cryosurgery of lung cancer using Modified Legendre wavelet Galerkin method, Journal of Thermal Biology, https://doi.org/10.1016/j.jtherbio.2018.10.012 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 study of cryosurgery of lung cancer using Modified Legendre wavelet Galerkin method Mukesh Kumar, Subrahamanyam Upadhyay and K.N. Rai

Abstract In this paper, we have developed a new mathematical model describing bioheat transfer during cryosurgery of lung cancer. The lung tissue cooled by a flat probe whose temperature decreases linearly with time. The freezing process occurs in three stages and the whole region is divided into solid, mushy and liquid region. The heat released in the mushy region is considered as discontinuous heat generation. The model is an initial-boundary value problem of the hyperbolic partial differential equation in stage 1 and moving boundary value problem of parabolic partial differential equations in stage 2 and 3. The method of the solution consists of four-step procedure as transformation of problem in dimensionless form, the problem of hyperbolic partial differential equation converted into ordinary matrix differential equation and the moving boundary problem of parabolic partial differential equations converted into moving boundary problem of ordinary matrix differential equations by using finite differences in space, transferring the problem into the generalized system of Sylvester equations by using Legendre wavelet Galerkin method and the solution of the generalized system of Sylvester equation are solved by using Bartels-Stewart algorithm of generalized inverse. The whole analysis is presented in dimensionless form. The effect of cryoprobe rate on temperature distribution and the effect of Stefan number on moving layer thickness is discussed in detail. Key words: Cryosurgery, DPL Model ,Modified Legendre wavelet Galerkin Method, Generalized inverse. 1. Introduction Cancer has become one of the most dangerous diseases in the world. Every year so many people die due to this disease. It occurs mainly in lungs, stomach, liver, prostate. Lung cancer is the most common cancer worldwide accounting for 1.59 million new cases annually and is responsible for 19.4 percent of all cancer deaths (World Health Organization). Main Preprint submitted to Elsevier

October 19, 2018

causes of lung cancer are cigarette smoking, passive smoking, some organic chemicals and air pollution etc. Cryosurgery, chemotherapy, hyperthermia and radiation therapy are the medical techniques for the treatment of the lung cancer. One of the most important medical techniques to treat the lung cancer is cryosurgery. Nomenclature x ρ c k a ρb cb wb T t L Tb T0 Tc T0 Tl Tz l si q Qm τq τT A fs g

space coordinate density(kg/m3 ) specific heat(J/kg 0 C) thermal conductivity of the tissue(W/m0 C) thermal diffusivity density of blood(kg/m3 ) specific heat of blood(J/kg 0 C) blood perfusion rate(ml/s/ml) temperature(0 C) time(s) latent heat(kJ/kg) arterial temperature(0 C) initial temperature(0 C) cryoprobe temperature(0 C) initial temperature(0 C) liquidus temperature(0 C) solidus temperature(0 C) length of the tissue distance from origin heat flux metabolic heat generation(W/m3 ) phase lag in heat flux(s) phase lag in temperature gradient(s) internal heat generation solid fraction cryoprobe rate

Dimensionless variable X dimensionless space coordinate Fo Fourier number or dimensionless time F oq dimensionless phase lag due to heat flux F oT dimensionless phase lag due to temperature gradient λi dimensionless distance 2

θ θb θu θm θf Pf Pm Ste Pd Subscript u, 1 m, 2 f, 3

dimensionless temperature dimensionless blood temperature dimensionless unfrozen temperature dimensionless mushy temperature dimensionless frozen temperature dimensionless blood perfusion coefficient dimensionless metabolic heat source coefficient Stefan number Predvoditelev number indication for unfrozen indication for mushy indication for frozen

Cryosurgery is the use of intense cold temperature in surgery to destroy diseased tissue [1]. In cryosurgery, Liquid nitrogen is introduced through an instrument called cryoprobe within diseased tissues [2]. Cryosurgery is a procedure that kills cancer cells by freezing them with liquid nitrogen. The main advantages of the cryosurgery are: (i) less cost, (ii) less pain , (iii) shorter hospital stay and (iv) shorter recovery time [3]. Immediate and Delayed effects are the two main mechanisms towards the destructive effect of freezing in cryosurgery. Cryosurgery not only destroying primary tumors but also distant disease [1]. The objective of cryosurgery in lung cancer is to maximize the damage of tumor tissues while preserving healthy lung tissues. The cryosurgery is an important treatment of the lung cancer since it provides the long-term survival. The detailed analysis of cell destruction and temperature distribution in the tumor during cryosurgery has been presented by Chua et al[4]. In their study, they observed that a single cryoprobe of bigger diameter is more effective as compared to multiple cryoprobes in destroying the tumor while preserving healthy tissue. Tarwidi investigated the Cryosurgical simulation of lung cancer based on efficient freezing time [5]. He showed that temperature distribution and phase change interface position in tissue can be used to maximize damage to tumor tissue and minimize the injury to normal tissue. Pennes bio-heat model[6] has been widely used by many authors[3, 4, 7, 8, 9] to solve the phase change heat transfer problem during cryosurgery which is given as cρ

∂T = − ▽ q + ρb wb cb (Tb − T ) + Qm ∂t 3

(1)

where, c and ρ are the specific heat and the density of the tissue, q is the heat flux, cb and ρb are the specific heat and density of blood, wb is the blood perfusion rate, Qm is the metabolic heat generation, Tb is the arterial blood temperature, T is the temperature and t is the time. The above Eqn.(1) is based on Fourier’s law which is defined as q(x, t) = −k ▽ T (x, t)

(2)

Although this model concludes the infinite speed of heat and mass propagation. However, in situations dealing with heat and mass transfer in extremely short periods of time, very high temperature and moisture gradient, very low temperatures and moisture approaching absolute zero, or for micro scale conditions the wave nature of heat and mass propagation becomes dominant. Cattaneo and Vernotte[10] proposed modified Fourier’s model(Single phage model) which are commonly known as non-Fourier’s model. For this type of situations, a constitutive equation which allows a relaxation time in the heat flux vector q. The single phage model is described by where τq ̸= 0 is the relaxation time in flux of heat propagation. Taking first order Taylor’s expansion, the SPL model reduces in the form ) ( ∂q (x, t) = −k ▽ T (x, t). (3) q + τq ∂t Many studies are available in the literature [7, 11, 12, 13, 14], where bioheat model with freezing has been used. Ahmadikia and Moradi[15] used the heat conduction model to describe the non-fourier effect of biological tissue during freezing. Further the physical solutions given by bio-heat model are sometimes unusual [16]. Therefore, to considered the thermal behavior which is not capture by the Fourier’s law, a new model by considering phase lag of temperature and phase lag of heat flux gradient is introduced by Tzou[16]. This model is called dual phase lag model(DPL) and is defined as q(x, t + τq ) = −k ▽ T (x, t + τT ), where τq ̸= 0 and τT ̸= 0 are the relaxation time in heat propagation. Taking first order Taylor’s expansion, the DPL model reduces in the form ( ) ( ) ∂q ∂T q + τq (x, t) = −k ▽ T + τT (x, t). (4) ∂t ∂t The two relaxation time τq (in heat flux) and τT (in temperature gradient) are the main characteristic of DPL model while relaxation time τq (in heat flux) is the main characteristic of SPL model. The study state occurs quickly in 4

the classical Fourier model than SPL model and occurs quickly in SPL model than DPL bio-heat transfer model [7]. In this paper, we have considered the system whose temperature is decreased with time by flat probe. The freezing process has occurred in three stages. (i) In the first stage, the tissue is cooled up to liquidus temperature Tl , (ii) In the second stage when tissue is cooled up to freezing temperature Tz , mushy region is formed. (iii) In this stage when the surface temperature of the flat probe is continuously decreasing, the freezing start and the frozen region, mushy region, and the unfrozen region is formed. The internal heat generation A in the mushy region[17, 18] is defined as: A=ρ∗L

∂fs ∂t

(5)

where fs is solid fraction present in mushy region. The solid fraction in the mushy region depends on the temperature of the mushy region. An exact mathematical expression for the solid fraction cannot be given but the approximate mathematical model is proposed by Gupta [19] for the numerical solution of the problem and defined as fs =

fsu (Tl − Tm ) − fsu1 (Tz − Tm ) Tl − Tz

(6)

where fsu and fsu1 are solid fractions present at liquid-mush and solid-mush boundaries respectively. At the liquid-mush boundary fsu taken as unity and fsu1 = 0, and at the solid-mush boundary fsu1 taken as unity and fsu = 0. So, In this paper we have considered that the solid fraction is also present at the liquid-mush boundary. Then, we have an implicit moving boundary condition which is a difficult task for numerical solution. We read so many problem of moving boundary from this book [20]. We also get the idea of moving boundary problem from these papers [21, 22]. During study of a moving boundary problem, we see that the Stefan number play the key role to determine the thickness of moving layer. The method of the solution consisting of four steps: (i) Transformation of problem in dimensionless form. (ii) Using finite differences in space variable, the initial value problem of hyperbolic partial differential equation and moving boundary problems of parabolic partial differential equations converted into a time boundary value problem of ordinary matrix differential equation and moving boundary problems of ordinary matrix differential equations respectively. 5

(iii) By making use of Legendre wavelet Galerkin method, the problem converted into the generalized systems of Sylvester equations. (iv) Using generalized inverse procedure (Bartels-Stewart algorithm) for the solution of the generalised system of Sylvester equation. The goal of the present study is to determine how the freezing occurs in three different stages. We studied the one-dimensional process of freezing in the different region: unfrozen region, mushy region, and frozen region. The rate of heat transfer determines the rate of tissue cooling. If there is a high temperature gap between cryogenic agent and tissues, rapid heat transfer occurs [1]. This occurs when the temperature of a cryogen like liquid nitrogen (−1960 C) contacts the skin (370 C)[1]. Freezing rate of tissues may be slow (< 100 C/min), moderate (10 − 1000 C/min) and fast (100 − 2600 C/min) . Since the objective is to destroy the tissues, it is essential to reach the fast cooling rate. In the tumor, the freezing interface moves rapidly as it passes from tumortissue into the normal lung tissue [23]. A point of inflection occurs when the frozen interface crosses the tumor-lung boundary into the healthy lung tissue corresponding to the time. A rapid fall of tumor temperature once the freezing interface has propagated slightly beyond the tumor-lung boundary into the lung, with a corresponding rapid rise in the lung temperature beyond the tumor-lung boundary shown in temperature profiles. These studies helps to provide the use of cryosurgery for the treatment of tumor-tissue in the lung and also suggest the technique for controlling the cryosurgery process. 2. Formulation Of Problem Our model is developed for the process of freezing of tumor tissue in the lung shown schematically in fig1. Freezing begins at the cryoprobe boundary and propagates into the tumor then into the lung. The location of the freezing interface is denoted by s1 , the size of the tumor by s2 , and the distance into the tissue by l. At time t = 0, the tissue is cooled by a flat probe whose temperature decreases with time. The cooling start and the whole process of cooling is composed in three stages. In stage I, the probe surface is cooled at a constant rate g, from the initial temperature T0 (370 C) to liquidus temperature Tl (−10 C). Further, in stage II the flat probe is continuously cooled at a constant rate g from the liquidus temperature Tl to freezing temperature Tz (−80 C). During this stage mushy region and liquid region simultaneously co-exist. In stage III, when the temperature of the flat probe is continuously decreased with time the freezing starts from the origin and propagates in the positive x-direction. In this stage frozen region, mushy region and unfrozen 6

Figure 1:

region are formed. The volumetric change in entire freezing process is negligible. Heat source appears due to metabolism and blood perfusion only when tissue is unfrozen [3]. Thermo-physical properties of healthy lung and tumor tissue[24] are taken as   kf , T < Tl k = 12 (kf + ku ), Tl ≤ T ≤ Tz ,   ku , T > Tz   ρf , T < Tl ρ = 21 (ρf + ρu ), Tl ≤ T ≤ Tz ,   ρu , T > Tz   c f , T < T l c = 21 (cf + cu ), Tl ≤ T ≤ Tz ,   cu , T > T z 7

In Stage I, the temperature of the tissue is above −10 C and during freezing, the relaxation time play important role and because of this non-fourier model(DPL) is considered[7]. While in Stage II and III Fourier model is used. From our simulations we observed that, the effect of relaxation time is presented only in unfrozen region and the effect of relaxation time present in mushy and frozen region is negligible. The mathematical model is described as follows: Stage I (

) ∂Tu ∂ 2 Tu ∂ 2 Tu ∂ 3 Tf ρu cu + τq 2 = ku + τ k T u ∂t1 ∂t1 ∂x2 ∂t1 ∂x2 ( ) ∂Tu ρb cb wb Tb − Tu − τq + Qm , 0 < X < l, t1 > 0 ∂t1

(7) (8)

2.1. Initial condition

(9)

Tu (x, 0) = To 2.2. boundary condition

Tu (0, t1 ) = T o − gt1 ,

∂Tu (l, t1 ) = 0 ∂x

(10)

Stage II

ρ m cm

∂Tm ∂ 2 Tm = km + A, ∂t2 ∂x2

ρ u cu

∂ 2 Tu ∂Tu = ku , ∂t2 ∂x2

0 < x < S2 , t 2 > t 1

(11)

S2 < x < l

(12)

2.3. Initial condition

Tm (x, t1 ∗) = Tu (x, t1 ∗) 8

(13)

2.4. boundary condition

Tm (0, t2 ) = Tl − gt2 ,

∂Tu (l, t2 ) = 0 ∂x

(14)

2.5. Interface condition

Tu (s2 , t2 ) = Tm (s2 , t2 ) = Tl ∂Tm ∂Tu km − ku = ρ ∗ Lvn ∂x ∂x

(15) (16)

where t1 ∗ =

T o − Tl , t2 = t1 − t1 ∗ g

Stage III

ρ f cf

∂Tf ∂ 2 Tf = kf , ∂t3 ∂x2

ρ m cm

0 < x < s1 , t 3 > t 2

∂Tm ∂ 2 Tm + A, = km ∂t3 ∂x2

ρ u cu

∂Tu ∂ 2 Tu , = ku ∂t3 ∂x2

s 1 < x < s2

s2 < x < l

(17)

(18)

(19)

2.6. Initial condition

Tf (x, t2 ∗) = Tu (x, t2 ∗) = Tm (x, t2 ∗)

(20)

2.7. boundary condition

Tf (0, t3 ) = Tz − gt3 ,

9

∂Tu (l, t3 ) = 0 ∂x

(21)

2.8. Interface condition

Tf (s1 , t3 ) = Tm (s1 , t3 ) = Tz , Tu (s2 , t3 ) = Tm (s2 , t3 ) = Tl

kf

(22)

∂Tf ∂Tm − km = ρ ∗ Lvn ∂x ∂x

(23)

∂Tm ∂Tu − ku = ρ ∗ Lvn ∂x ∂x

(24)

km where

t2 ∗ =

Tl − Tz , t3 = t2 − t2 ∗ g

3. Dimensionless Analysis Dimensionless variable and similarity criteria defined as follows: x T0 − Tu T0 − Tm T0 − Tf , θu = , θm = , θf = , l T0 − Tl T0 − Tl T0 − Tl a2 t2 a3 t3 T0 − Tb T0 − Tz a1 t1 θb = , θz = , F o1 = 2 , F o2 = 2 , F o3 = 2 , T0 − Tl T0 − Tl l l l 2 2 a1 τ q a1 τ T Qm l ρ b ∗ cb w b l F oq = 2 , F oT = 2 , Pf2 = , Pm = , l l k k(Tl − T0 ) ( ) ( ) ( ) T0 − Tl T0 − Tl Tl − Tz Ste3 = c3 , Ste2 = c2 , Ste1 = c2 L L L 2 2 2 gl gl gl P d1 = , P d2 = , P d3 = a1 (T0 − Tl ) a2 (T0 − Tl ) a3 (Tl − Tz ) ku km ku , kum = , kmf = , a1 = , km kf ρ ∗ c1 kM kf sj (ti ) a2 = , a3 = , λj (F oi ) = , j = 1, 2, i = 2, 3. ρ ∗ c2 ρ ∗ c3 l

X=

The system of equations(7)-(24) reduce in the following form:

10

Stage I ( 2 ) ∂θu ∂θu ∂ 2 θu ∂ 2θ ∂ ∂ θ 2 + F oq + Pf ∗ F oq = + F oT 2 2 ∂F o1 ∂F o1 ∂F o1 ∂X ∂F o1 ∂X 2 +Pf2 (θb − θu ) − Pm

(25)

initial and boundary conditions are:θu (X, 0) = 0, θu (X, F o∗ ) = 1 θu (0, F o1 ) = P d1 ∗ F o1 ∂θu (1, F o1 ) = 0 ∂X

(26) (27) (28)

Stage II ( 1+

1 Ste2

)

∂θm ∂ 2 θm , = ∂F o2 ∂X 2 ∂θu ∂ 2 θu = a12 , ∂F o2 ∂X 2

0 < X < λ2 (F o2 )

(29)

λ2 (F o2 ) < X < 1

(30)

where a12 = aa21 initial and boundary condition are given as:θm (X, F o∗1 ) = θu (X, F o∗1 ), θm (0, F o2 ) = P d2 ∗ F o2 θm (λ2 (F o2 ), F o2 ) = θu (λ2 (F o2 ), F o2 ) = 0 ∂θu (1, F o2 ) = 0. ∂X

(31) (32) (33) (34)

Interface condition is given by: ∂θm ∂θu −1 ∂X − kum = , X = λ2 (F o2 ) ∂X ∂X Ste2 ∂F o2 where F o∗1 =

1 , P d1

F o2 = F o1 − F o∗1

11

(35)

Stage III

∂θf ∂ 2 θf = , ∂F o3 ∂X 2 ( ) 1 ∂θm ∂ 2 θm 1− = a23 , Ste1 ∂F o3 ∂X ∂ 2 θu ∂θu = a13 , ∂F o3 ∂X

0 < X < λ1 (F o3 )

(36)

λ1 (F o3 ) < X < λ2 (F o3 )

(37)

λ2 (F o3 ) < X < 1

(38)

where a23 = aa32 ,a13 = aa31 initial and boundary condition are given as:θm (X, F o∗2 ) = θf (X, F o∗2 ) = θu (X, F o∗2 ), θf (0, F o3 ) = θz + P d3 ∗ F o3 θf (λ1 (F o3 ), F o3 ) = θm (λ1 (F o3 ), F o3 ) = θz ∂θu (1, F o3 ) = 0 ∂X

(39) (40) (41) (42)

Interface condition is given by: ∂θf ∂θm −1 ∂X − kmf = , X = λ1 (F o3 ) ∂X ∂X Ste3 ∂F o3 θm (λ2 (F o3 ), F o3 ) = θu (λ(F o3 ), F o3 ) = 1 ∂θm ∂θu −1 ∂X − kum = , X = λ2 (F o2 ) ∂X ∂X Ste3 ∂F o2 where F o∗2 =

1 , P d2

(43)

(44)

(45)

F o2 = F o2 − F o∗2 .

4. Solution of the Problem Using finite differences in space and after simple algebraic calculation, our problem reduces as follows: For Stage I, F oq

) dΘu ( 2 ) d2 Θu ( 2 + 1 + P ∗ F o − F o U + P I − U Θu q T 1 1 f f dF o21 dF o1 P d1 F oT P d1 F o1 = d+ d + (Pf2 θb + Pm )d1 2 2 h h 12

(46)

Θu (0) = 0, Θu (F o∗ ) = 1,

(47) (48)

where Θu ,U1 , d and d1 are given by: 

     −2, 1, 0, .......0, 0, 0 1 1  1, −2, 1, .......0, 0, 0    0 1  0, 1, −2, .......0, 0, 0      0 1 1  .. .. .. ..        Θu =  ..  , U1 = 2  ., . . . , .; . . . ; ., . . . , .  , d1 =  ..  . , d =  ..   .  . .  h . . . .  .., . . . , ..; . . . ; .., . . . , ..   .  . .    ..   ..   ..  0, 0, 0, ......., 1, −2, 0 N θu N ×1 0 N ×1 1 N ×1 0, 0, 0, ......., 0, 2, −2 N ×N 

 θu1  θu2   3  θu   

For Stage II, ) ( P d2 F o2 ¯ dΘm U2 1 d, 0 < x < λ2 (F o2 ), (49) − 2 Θm = 1+ Ste2 dF o2 h h2 dΘu a12 U3 − Θu = 0, λ2 (F o2 ) < X < 1 (50) dF o2 h2 Θm (F o∗1 ) = Θu (F o∗1 ) (51) ¯ − 1 × 1 column vector in which first entry is one and remaining where d¯ is N are zeros, Θm , Θu , U2 and U3 are given as follows:  ¯   1  θuN +1 θm ¯ +2   2 N  θm  θu   3   N¯ +3   θm  θu   .   , , Θu =  Θm =  .   ...   .     .   ..   ..   .  ¯ −1 N θm θuN N −N¯ ×1 ¯ −1×1 N     −2, 1, 0, .......0, 0, 0 −2, 1, 0, .......0, 0, 0  1, −2, 1, .......0, 0, 0   1, −2, 1, .......0, 0, 0       0, 1, −2, .......0, 0, 0   0, 1, −2, .......0, 0, 0   .  .. .. ..  .. .. ..   ..  . , U3 =  ., . . . , .; . . . ; ., . . . , .  U2 =  ., . . . , .; . . . ; ., . . . , .  .  .   .., . . . , ...; . . . ; ..., . . . , ...   .., . . . , ...; . . . ; ..., . . . , ...      0, 0, 0, ......., 1, −2, 0 0, 0, 0, ......., 1, −2, 0 0, 0, 0, ......., 0, 2, −2 N −N¯ ×N −N¯ 0, 0, 0, ......., 0, 1, −2 N¯ −1×N¯ −1 Interface condition is given as: N −1 1 − θm θN +1 − 1 −1 dλ2 − kum u = . h h Ste1 dF o2 ¯

¯

13

(52)

For Stage III, dΘf A′ θz P d3 F o3 ′ − 21 Θf = 2 B0′ + B1 , dF o3 h h h2 ( ) 1 dΘm a23 A′2 1 1− − Θm = 2 B2′ , 2 Ste1 dF o3 h h ′ dΘu a13 A3 1 − Θu = 2 B3′ , 2 dF o3 h h ∗ ∗ Θm (F o2 ) = Θu (F o2 ) = Θf (F o∗2 )

0 < x < λ1 (F o3 )

λ1 (F o3 ) < X < λ2 (F o3 ) (54) λ2 (F o3 ) < X < 1

where, 

 θf1  θ2   f   θ3   f   Θf =   ...     ..   .  θfp¯−1



 1  1  θm θu 2   θu2   θm  3   3   θm   θu      , Θm =  ..  , Θu =  ..  ,  .   .   .   .   ..   ..  q¯−1 θm θur¯−1 r¯−1×1 q¯−1×1 p¯−1×1     −2, 1, 0, .......0, 0, 0 −2, 1, 0, .......0, 0, 0  1, −2, 1, .......0, 0, 0   1, −2, 1, .......0, 0, 0       0, 1, −2, .......0, 0, 0   0, 1, −2, .......0, 0, 0  . . . . . . . .     A′1 =  .., . . . , ..; . . . ; .., . . . , ..  , A′2 =  .., . . . , ..; . . . ; .., . . . , ..  , .     .., . . . , ...; . . . ; ..., . . . , ...   ..., . . . , ...; . . . ; ..., . . . , ...      0, 0, 0, ......., 1, −2, 0 0, 0, 0, ......., 1, −2, 0 0, 0, 0, ......., 0, 1, −2 p¯−1ׯp−1 0, 0, 0, ......., 0, 1, −2 q×q

14

(53)

(55) (56)



   −2, 1, 0, .......0, 0, 0 1  1, −2, 1, .......0, 0, 0    1  0, 1, −2, .......0, 0, 0    .  1 . . .     A′3 =  .., . . . , ..; . . . ; .., . . . , ..  , B0′ =  ..  , . .   .., . . . , ...; . . . ; ..., . . . , ...  .    ..  0, 0, 0, ......., 1, −2, 0 1 p×1 0, 0, 0, ......., 0, 2, −2 r×r       1 0 1 0 0 0       0 0 0 . .   ′ ′ ′ B1 =  .  , B2 =  .  , B3 =  ..  , . .     . . . .  ..   ..   ..  1 q×1 0 r×1 0 p×1 Interface conditions are given as: −8 − θfp (λ1 ) −1 dλ1 θq (λ1 ) + 8 − kmf m = h h Ste3 dF o3

(57)

q θr (λ2 ) + 8 1 dλ2 −8 − θm (λ2 ) − kum u = h h Ste3 dF o3

(58)

5. Modified Legendre wavelet Galerkin method 5.1. Stage 1 We assume that the unknown function

d2 Θu dF o1 2

d2 Θu ∼ = Cψ(F o1 ), dF o1 2

is approximated by (59)

where C is unknown matrix of order N × 2k−1 M and ψ(F o1 ) is a column vector of order 2k−1 M − 1 × 1 defined by ψ(F o) = [ψ1,0 , ψ1,1 , . . . , ψ1,M −1 , . . . , ψ2k−1 ,0 , ψ2k−1 ,1 , . . . ψ2k−1 ,M −1 ]′ where the elements of ψ(F o) is defined by {√ ( ) n ˆ +1 (m + 1/2)2k/2 Pm 2k F o1 − n ˆ , nˆ2−1 k ≤ F o1 ≤ 2k ψn,m (F o1 ) = 0 , otherwise. where k = 1, 2, 3..., n = 1, 2, . . . , 2k−1 , n ˆ = 2n − 1 and m is the order of Legendre polynomial [25]. Pm (F o1 ) is denoted by Legendre polynomial of 15

order m, m = 0, 1, ..., M − 1 , which are orthogonal with respect to the weight function w(F o1 ) = 1 on the interval [−1, 1], and satisfy the following recursive formula P0 (F o1 ) = 1, P1 (F o1 ) = F o1 , 2m + 1 m Pm+1 (F o1 ) = F o1 Pm (F o1 ) − Pm−1 (F o1 ). m+1 m+1 The operational matrix of integration of ψ defined in [22, 26, 27] and given by ∫F o1 ψ(t)dt = P ψ(F o1 ) (60) 0

where P if is the operational matrix of integration of order 2k−1 M × 1, given by   L O O ··· O  0 L O · · · O    1   0 0 L · · · O P = k  .. .. .. . . ..  2 . . . . .    0 0 0 · · · O 0 0 0 ··· O   2 0 ··· 0 0 0 · · · 0   where F and L are M × M matrices given by O =  .. .. . . ..  . . . . 0 0 ··· 0 and   1 √13 0 0 ··· 0 0 0 −1 1 √  0 0 0  3 0 √15 0 · · ·   0 √−1  1 √ 0 · · · 0 0 0   15 35   0 √−1 0 · · · 0 0 0  0 35 L=  . . . . . . . ...   . . . . . . . . . . . . . .     1 1 √ 0 0 0 0 ··· √  0 (2M −3)(2M −5) (2M −3)(2M −1)   −1 √ 0 0 0 0 ··· 0 0 (2M −1)(2M −3)

respectively.

16

Integrating (59) two time over 0 to F o, we obtained dΘu ∼ dΘu (0) + P ψ(F o1 ), (61) = dF o1 dF o1 dΘu Θu ∼ (0)F o1 + CP 2 ψ(F o1 ) (62) = Θu (0) + dF o1 Substituting F o∗ in the equation and using boundary condition, we get dΘu ∼ d1 − CP 2 ψ(F o∗ ) + CP ψ(F o1 ), = dF o1 F o∗ ( ) d1 − CP 2 ψ(F o∗ ) ∼ Θu = F o1 + CP 2 ψ(F o1 ) ∗ Fo

(63) (64)

( ) ( ) P d1 F F oq Cψ + I + Pf2 ∗ F oq I − F oT U CP ψ + Pf2 I − U CP 2 ψ = ψ h2 P d1 F oT )Gψ (65) +(Pf2 θb + Pm + h2 where d1 = Gψ and F o1 d = F ψ ( ) ( ) F oq C + I + Pf2 ∗ F oq I − F oT U CP + Pf2 I − U CP 2 = P d1 F P d1 F oT + (Pf2 θb + Pm + )G (66) 2 h h2 5.2. Stage 2 dΘm dΘu We assume that the unknown functions dF and dF are approximated o2 o2 by dΘm ∼ (67) = D1 ψ, dF o2 dΘu ∼ (68) = D2 ψ, dF o2 Integrating (67,68) over F o∗ to F o2 , we obtained Θm ∼ (69) = D1 P ψ − D1 P (ψ (F o∗1 ) − ψ (0)) , ∗ ∼ (70) Θu = D2 P ψ − D2 P (ψ (F o ) − ψ (0)) , 1

After solving Eqns(49),(67) and (69), we are getting 1 P d2 1 )D1 I + D1 U2 P ψ(F o∗1 )F − D1 U2 P = 2 dG + 2 dF + Ste1 h h 1 d1 F + U2 d2 F h2 After solving Eqns(50),(68) and (70), we are getting a12 a12 1 D2 I + a12 D2 U3 P ψ(F o∗1 )F − 2 D2 U3 P = 2 U3 d2 F + 2 d1 F h h h (1 +

17

(71)

(72)

5.3. Stage 3 We assume that the unknown functions mated by

dΘf , dΘm dF o3 dF o3

and

dΘu dF o3

are approxi-

dΘf ∼ = E1 ψ, dF o3 dΘm ∼ = E2 ψ, dF o3 dΘu ∼ = E3 ψ, dF o3

(73) (74) (75)

Integrating (73-75) over F o∗ to F o3 , we obtained Θf ∼ = E1 P ψ − E1 P (ψ (F o∗2 ) − ψ (0)) , Θm ∼ = E2 P ψ − E2 P (ψ (F o∗2 ) − ψ (0)) , Θu ∼ = E3 P ψ − E3 P (ψ (F o∗2 ) − ψ (0))

(76) (77) (78)

After solving Eqns(53),(73) and (76), we are getting E1 I +

1 1 θz P d3 8 ′ ′ ∗ ′ E A P ψ(F o )F − E A P = dF + dG − A 1 1 1 2 1 h2 h2 h2 h2 h2 1

(79)

After solving Eqns(54),(74) and (77), we are getting E2 (1 −

1 a23 1 8 a23 )I + 2 E2 A′2 P ψ(F o∗2 )F − 2 E2 A′2 P = 2 d1 F − a23 2 A′2 (80) Ste2 h h h h

After solving Eqns(55),(75) and (78), we are getting E3 I +

a13 a13 1 8 E3 A′3 P ψ(F o∗2 )F − 2 E3 A′3 P = 2 dF − 2 A′3 2 h h h h

(81)

The problem in all the three Stages converted into generalized system of Sylvester equations (66, 71, 72, 79−81) whose solution obtained using BartelsStewart algorithm [28]. Using this result in interface condition, the position of moving layer thickness in Stage 2 and Stage 3 determined. For numerical computation, we have used Matlab Software. 6. Result and Discussion We developed the time relaxation model for freezing of tumor tissue in the lung. The process of freezing is completed in three stages defined in Eq. (7-24). For simplicity of the model we introduced the dimensionless parameters defined in section 3. Therefore, the dimensionless form of the 18

Table 1: Thermal-physical properties of tissues(Bischof et al.[23])

Parameter unit of measurement value Density of unfrozen lung tissue kg/m3 161 3 Density of frozen lung tissue kg/m 149 Density of unfrozen tumor tissue kg/m3 998 3 Density of frozen tumor tissue kg/m 921 3 Density of blood kg/m 1005 Thermal conductivity of unfrozen lung tissue W/m0 C 0.11 0 Thermal conductivity of frozen lung tissue W/m C 0.38 0 Thermal conductivity of unfrozen tumor tissue W/m C 0.552 Thermal conductivity of frozen tumor tissue W/m0 C 2.25 0 Specific heat of unfrozen lung tissue J/kg C 4174 Specific heat of frozen lung tissue J/kg 0 C 1221 0 Specific heat of unfrozen tumor tissue J/kg C 4200 0 Specific heat of frozen tumor tissue J/kg C 1230 Blood perfusion in lung tissue ml/s/ml 0.0005 Blood perfusion in tumor tissue ml/s/ml 0.002 Metabolic heat generation in lung W/m3 42, 000 Metabolic heat generation in tumor W/m3 672 Latent heat kJ/kg 333 0 Arterial blood temperature C 37 model is given in Eq. (25-45). For the solution of a non-dimensional form of model we used Modified Legendre wavelet Galerkin method in each stage.In stage II and stage III, the model is a moving boundary problem of partial differential equations. In stage 2, Eq. (35) represents the interface condition for mushy-liquid region and in stage 3, Eq.(43, 45) represent the interface condition for solid-mushy and mushy-liquid region respectively. For determing the moving layer thickness, we used interface condition in stage 2 and 3 respectively and obtain the values of λ2 (F o2 ) in stage 2 and λ1 (F o3 ), λ2 (F o3 ) in stage 3. The effect of Stefan number on Moving layer thickness in both stages is observed. During the cryosurgical treatment of lung-tumor tissue, temperature distribution and moving layer thickness are important for the prediction of maximum damage to diseased tissue and minimum damage to healthy lung tissue. Therefore, we have examined the temperature distribution in all stages and moving layer thickness in mushy and frozen region. For numerical computation, we have taken the parameter of lung and tumor tissue from table 1[8, 23] given below:

19

Stage I: In stage I, Fig. 2 shows the graph between temperature and time for the different value of g(rate of cryoprobe). In this case we have taken dimensionless relaxation time F oq = 0.001 and F oT = 0.015.As cryoprobe rate g increases the temperature decreases rapidly and temperature also decreases when the time increases. That means the lung tumor is cooled from initial temperature T0 (370 C) to liquidus temperature Tl (−10 C), the temperature (Tu ) decreases as space coordinate x increases. We find the non dimensional temperature Θu by using Modified Legendre wavelet Galerkin Method and then find Tu . In Fig. 4 we have taken different values of relaxation times F oq and F oT . We see that for different values of relaxation time, temperature decreases accordingly. We find the difference between the Fourier and non-Fourier model in our graph. For Fourier model temperature decreasing rate is less than in comparable to non-Fourier model. Stage II: In stage II, the solid fraction in the mushy region depends on the temperature of the mushy region. At the liquid-mush boundary fsu taken as 1 and fsu1 taken as 0. Fig. 5 shows the graph between temperature and time of both mushy(Tm ) and unfrozen(Tu ) region. We see that the temperature decreases as time increases. We also see that temperature of mushy region decreases fastly as comparable to the unfrozen region and it is clear form the figure that the temperature in mushy region decreases at x = 0.3 with the time increases and in the liquid region at x = 0.5. Here the lung tumor is cooled from the liquidus temperature Tl (−10 C) to freezing temperature Tz (−80 C). We find the non-dimensional temperature Θm and Θu by using Modified wavelet Galerkin Method and then find Tm and Tu . We also see that the temperature of the mushy and unfrozen region decreases as space coordinate x increases. Fig. 7 shows the dimensionless moving layer thickness λ2 (F o2 ) in mushy-liquid region and it increases as the dimensionless time F o2 increases. We also see that how Stefan number effects the moving layer thickness in figure 13. We also see that as the cryoprobe rate g increases the temperature decreases rapidly and moving layer thickness increases. Stage III: In this stage, the freezing moving front λ1 and the liquidus moving front λ2 are calculated. Fig. 9 shows that how temperature decreases with the time increases in each region. We also see the effect of cryoprobe rate g in this region. Here we see that the freezing temperature Tf decreases rapidly as the g increases. We also see the effect of Stefan number on moving layer 20

thickness λ1 and λ2 in fig. 14 and 15 respectively. We see that the temperature of each region decreases as space coordinate x increases. Fig.15 shows that the dimensionless moving layer thickness λ1 (F o3 ) in solid-mush region and Fig.16 shows that the moving layer thickness λ2 (F o3 ) in mush-liquid region increases as the dimensionless time F o3 increases. We find the nondimensional temperature Θf , Θm and Θu by using Modified wavelet Galerkin Method and then find Tf , Tm and Tu . We also see that the temperature of frozen, mushy and unfrozen region decreases as space coordinate x increases. Fig.15,16 shows the dimensional moving layer thickness λ1 (F o3 ) in solidmush region and λ2 (F o3 ) in mush-liquid region respectively. Both moving layer thickness increases as the dimensionless time F o3 increases. We also see that how Stefan number effects the moving layer thickness, as Stefan number decreases moving layer thickness increases of both λ1 and λ2 . In this stage, we cooled the tissue from freezing temperature to the lethal temperature. Lethal temperature for tissue destruction in particular the diseased tissue like tumor usually begins around −400 C[13]. We also see that as the cryoprobe rate g increases the temperature decreases. 7. Conclusion In this paper, we have developed the time relaxation model of lung tumor tissue during the freezing process in cryosurgery and then developed a new Modified Legendre wavelets Galerkin method to solve our problem. Firstly, our problem is converted into dimensionless form and with further application of finite difference method, the problem converted into time boundary value problem of ordinary matrix differential equation in stage 1 and moving boundary value problem in stage 2 and 3. Later Legendre wavelet Galerkin method is applied in which we obtained the system of Sylvester equations which are solved by generalised inverse technique(Bartels-Stewart algorithm). In each stage we have used this method to find the temperature. From our simulations, we observed that the effect of relaxation time is present only in unfrozen region and the effect of relaxation time, present in the mushy and frozen region is negligible. Further, we have compared the Fourier model with a non-Fourier model in unfrozen region. We observe from the figure that the cryoprobe rate effect the temperature in all the stages and as the time increases the temperature decreases rapidly whose cryoprobe rate is faster than the other one. In the first Stage, the tissue is cooled up to the liquidus temperature Tl (−10 C) with a fast cooling rate and then in Stage 2, the tissue is cooled up to the freezing temperature Tz (−80 C) where the mushy region is formed and in the last Stage, the tissue is cooled up to the lethal temperature where the frozen region is formed. Stefan number also 21

effects the moving layer thickness in both the stages 2 and 3. Although, we have considered one dimensional dual phase lag bio-heat model in the present investigation, it can be extended for two dimensional bio-heat model to achieve more physically realistic results. Our model will be very useful for the experimental analysis of freezing process of cryosurgical treatment in lung cancer. Acknowledgments The authors expresses their sincere thanks to reviewers for their valuable suggestions in the improvement of this article. Authors are also thankful to DST-Centre for Interdisciplinary Mathematical Sciences, Banaras Hindu University Varanasi, India, for providing necessary facilities. References [1] Paola Pasquali, Cryosurgery: A Practical Manual(2015). [2] H.V. Allington, Liquid nitrogen in the treatment of skin diseases, Calif. Med, 72 (3),153-155(1950). [3] Z.S. Deng, J. Liu, Numerical simulation of selective freezing of target biological tissues following injection of solutions with specific thermal properties, Cryobiology, 50(2),183-192(2005) . [4] K.J. Chua, S.K. Chou, J.C. Ho, An analytical study on the thermal effects of cryosurgery on selective cell destruction, J. Biomech,40(1), 100–116(2007). [5] D. Tarwidi, Godunov method for computerized lung cancer cryosurgery planning with efficient freezing time, In Information and Communication Technology (ICoICT), 2015 3rd International Conference on IEEE, pp. 494–499(2015). [6] H.H Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm. J. Appl. Physiol. 1 (2), 93–122 (1948). [7] D. Kumar,S. Singh,K.N. Rai, Analysis of classical fourier, SPL and DPL heat transfer model in biological tissues in presence of metabolic and external heat source. Heat Mass Transf. 52 (6), 1089–1107(2016). [8] Ajay Kumar, Sushil Kumar, V.K. Katiyar, Shirley Telles, Dual phase lag bio-heat transfer during cryosurgery of lung cancer: Comparison of three heat transfer models, Journal of Thermal Biology,69(2017). 22

[9] P kumar, D kumar, K.N. Rai, A numerical study on dual phase lag model of bio heat transfer during hyperthermia treatment, Journal of thermal biology 49,98-105(2015). [10] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, Comptes Rendus, vol. 247, pp. 431–433, 1958. [11] Xin Zeng and Amir Faghri, Temperature transforming model for binary solid-liquid phase change problem part I: Mathematical modelling and numerical methodology, Numerical Heat Transfer B, vol. 25, pp. 467– 480, 1994a. [12] Xin Zeng and Amir Faghri, Temperature transforming model for binary solid-liquid phase change problem part II: Numerical simulation, Numerical Heat Transfer B, vol. 25, pp. 481–500,1994b. [13] Xin Zhang, S.M. Chapal Hossain, Gang Zhao, beisheng Qiu, Xiaoming He, Two-phase heat transfer model for multiprobe cryosurgery, Applied Thermal Engineering,113,47-57(2017). [14] P.K. Gupta, J Singh, K.N.Rai, Numerical simulation for heat transfer in tissues during thermal therapy, Journal of Thermal Biology,35(6),295301(2017). [15] H. Ahmadikia, A. Moradi, Non-Fourier phase change heat transfer in biological tissues during solidification. Heat. Mass Transf. 48 (9), 1559– 1568 (2012). [16] D.Y. Tzou, A unified field approach for heat conduction from macro to microscales Transc,117,8-16(1995). [17] S. Yadav, S. Upadhyay, K.N. Rai, A Mathematical Model for solidification of Binary Eutectic System including relaxation time, Computational Thermal Sciences,8,583-589(2016). [18] K.N. Rai, S.K. Rai, Boundary fixation approach to one dimensional inverse Stefan problem in non-ideal biological tissues, Fourth National Conference on Thermal system ,290-297(2003). [19] SC Gupta, The Classical Stefan Problems: Basic Concepts,Modelling and Analysis, Elsevier, 2003. [20] J. Crank, Free and Moving Boundary Problems, Clarendon Press, 1987. 23

[21] S. Yadav, S Upadhyay, K.N. Rai, Wavelet Galerkin and Wavelet Collocation method in moving boundary problem with temperature dependant thermal physical properties, Preceeding of Conv-14 Int,Symp.on Convective hheat and mass Transfer,790,1061-1079(2014). [22] S Upadhyay, S Yadav, K.N. Rai, Modelling and Simulation of a Moving Boundary Problem arising during immersion frying of foods, National Academic Science Letters(Accepted-2017). [23] J.C. Bischof,J.Bastacky,B. Rubinsky, An analytical study of cryosurgery in the lung. J. Biomech. Eng. 114 (4), 467–472 (1992). [24] S. Kumar, V.K. Katiyar, Numerical study on phase change heat transfer during combined hyperthermia and cryosurgical treatment of lung cancer, Int. J. Appl. Math. Mech. 3 (3) 1–17(2007). [25] S Yadav, D Kumar, K.N. Rai, Finite Element Legendre Wavelet Galerkin Approach to inward solidification in simple body under most generalised boundary condition, Z. Naturforsch. 69a, 501-510(2014). [26] M. Razzaghi and S.Yousefi, Legendre wavelets operational matrix of integration, International Journal of System Science,vol. 32, pp. 495– 502, (2001). [27] S Upadhyay, K.N. Rai, Collocation method applied to unsteady flow of gas through a porous medium, Int.J. App Math Research,3(3) 251259(2014). [28] R. Bartels and G. Stewart, Solution of the Matrix Equation AX + XB = C: Algorithm 432, Comm. ACM, 15 ,pp.820–826(1972).

24

Figure 2:

25

Figure 3:

26

Figure 4:

27

Figure 5:

28

Figure 6:

29

Figure 7:

30

Figure 8:

31

Figure 9:

Figure 10:

32

Figure 11:

33

Figure 12:

34

Figure 13:

35

Figure 14:

Figure 15:

36

Figure 16:

37