A study of heat transfer during cryosurgery of lung cancer

Accepted Manuscript A study of heat transfer during cryosurgery of lung cancer Mukesh Kumar, Subrahamanyam Upadhyay, K.N. Rai PII:

S0306-4565(19)30031-2

DOI:

https://doi.org/10.1016/j.jtherbio.2019.05.023

Reference:

TB 2341

To appear in:

Journal of Thermal Biology

Received Date: 18 January 2019 Revised Date:

4 May 2019

Accepted Date: 22 May 2019

Please cite this article as: Kumar, M., Upadhyay, S., Rai, K.N., A study of heat transfer during cryosurgery of lung cancer, Journal of Thermal Biology (2019), doi: https://doi.org/10.1016/ j.jtherbio.2019.05.023. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A study of Heat Transfer during cryosurgery of lung cancer

RI PT

Mukesh Kumar, Subrahamanyam Upadhyay and K.N. Rai

SC

Abstract

In this study, a mathematical model describing two-dimensional bio-heat transfer during cryosurgery of lung cancer is developed. The lung tissue is

M AN U

cooled by a cryoprobe by imposing its surface at a constant temperature or a constant heat flux or a constant heat transfer coefficient. The freezing starts and the domain is distributed into three stages namely: unfrozen, mushy and frozen regions. In stage I where the only unfrozen region is formed, our problem is an initial-boundary value problem of the hyperbolic partial differ-

TE D

ential equation. In stage II where mushy and unfrozen regions are formed, our problem is a moving boundary value problem of parabolic partial differential equations and in stage III where frozen, mushy, and unfrozen regions are formed, our problem is a moving boundary value problem of parabolic par-

EP

tial differential equations. The solution consists of the three-step procedure: (i) transformation of problem in non-dimensional form, (ii) by using finite

AC C

differences, the problem converted into ordinary matrix differential equation and moving boundary problem of ordinary matrix differential equations, (iii) applying Legendre wavelet Galerkin method the problem is transferred into the generalized system of Sylvester equations which are solved by applying Bartels-Stewart algorithm of generalized inverse. The complete analysis is presented in the non-dimensional form. The consequence of the imposition of boundary conditions on moving layer thickness and temperature distribution

Preprint submitted to Elsevier

May 23, 2019

ACCEPTED MANUSCRIPT

are studied in detail. The consequence of Stefan number, Kirchoff number and Biot number on moving layer thickness are also studied in specific.

RI PT

Key words: Cryosurgery, DPL model, FEM, Generalized boundary condition, Generalized inverse. 1. Introduction

SC

Cancer has turned out to be one of the most threatening diseases all over the globe. Presently 90.5 million population is suffering from cancer and

M AN U

as estimated from studies by 2020 around 20 million fresh cancer cases and 10 million cancer deaths will occur annually(World Health Organization). Diet and obesity (30-35 percent), Tobacco(25-30 percent), infections(15-20 percent), radiation (up to 10 percent), mental pressure, less body movement, and pollution are common contributes to cancer death (Islami et al.,2018).

TE D

However, it is still difficult to specify an exact cause behind a certain type of cancer due to the existence of a large number of factors in environment. for e.g. tobacco is not only a factor for lung cancer, air pollution or radiation also contributes to its growth in human body (Tolar and Neglia, 2003). Cancer

EP

mainly effects the lungs, stomach, liver, prostate, and colorectal among all other body parts.

AC C

Nomenclature x

space coordinate

y

space coordinate

ρ

density(kg/m3 )

c

specific heat(J/kg 0 C)

k

thermal conductivity of the tissue(W/m0 C)

a

thermal diffusivity 2

ACCEPTED MANUSCRIPT

density of blood(kg/m3 )

cb

specific heat of blood(J/kg 0 C)

wb

blood perfusion rate(ml/s/ml)

T

temperature(0 C)

t

time(s)

L

latent heat(kJ/kg)

Tb

arterial temperature(0 C)

T0

initial temperature(0 C)

Tc

cryoprobe temperature(0 C)

Tl

liquidus temperature(0 C)

Ts

solidus temperature(0 C)

l

length of the tissue

si

distance from origin

q

heat flux

Qm

metabolic heat generation(W/m3 )

τq

phase lag in heat flux(s)

τT

phase lag in temperature gradient(s)

Tw

surface temperature

h

SC

M AN U

TE D

surrounding temperature heat transfer coefficient error function

AC C

erf

EP

Tr

RI PT

ρb

erf c

complementary error function

Dimensionless variable X

dimensionless space coordinate

Y

dimensionless space coordinate

Fo

Fourier number or dimensionless time

3

ACCEPTED MANUSCRIPT

dimensionless phase lag due to heat flux

F oT

dimensionless phase lag due to temperature gradient

λi

dimensionless distance

θ

dimensionless temperature

θb

dimensionless blood temperature

θu

dimensionless unfrozen temperature

θm

dimensionless mushy temperature

θf

dimensionless frozen temperature

θw

dimensionless surface temperature

θr

dimensionless surrounding temperature

θs

dimensionless solidus temperature

Pf

dimensionless blood perfusion coefficient

Pm

dimensionless metabolic heat source coefficient

Ste

Stefan number

Ki

Kirchoff number

Bi

Biot number

m, 2

TE D

M AN U

SC

u, 1

EP

Subscript

indication for unfrozen

indication for mushy indication for frozen

AC C

f, 3

RI PT

F oq

Lung cancer is one amongst the commonly found cancer globally com-

prising 1.69 million new patients yearly and is cause for 22 percent of total deaths (WHO) due to cancer. Cigarette smoking, second-hand smoke, exposure to radon and air pollution etc are some of its main factors. There is no assured method to intercept lung cancer but danger can be decreased if we: (i) Stop smoking, (ii) Avoid second-hand smoke, (iii) Check our home 4

ACCEPTED MANUSCRIPT

for radon, (iv) Include fresh fruits and vegetables in our diet, (v) make exercise on the daily routine. Although it can be cured by some therapies too

RI PT

commonly known as cryosurgery, chemotherapy, radiation therapy, and hyperthermia. But amongst these ones of the most important is cryosurgery (Kumar et al., 2018).

Cryosurgery can be stated as the destruction of diseased tissue applying

SC

intense cold temperature in surgery (Paola, 2015). Here, diseased tissues are

injected by an instrument called cryoprobe with liquid nitrogen (Allington, 1950). Cancer cells are frozen with liquid nitrogen or argon gas and by

M AN U

the time killed in the process of cryosurgery. There are several benefits of the cryosurgery such as (i) minimum cost incurred, (ii) minimum pain suffered, (iii) shorter stay in the hospital and (iv) shorter recovery time (Deng and Liu, 2005). In lung cancer, the main objective of cryosurgery is to increase the destruction of tumor tissues while protecting enclosed healthy

TE D

lung tissues. Cryosurgery is considered to be a beneficial treatment of lung cancer as it increases the possibilities of longer survival (Kumar et al., 2018). For successful cryosurgical treatment, temperature distribution and positions of phase change interface in lung-tumor tissue are required (Katiyar et al.,

EP

2017). The comprehensive study has been presented regarding the rate of cell destruction and temperature distribution in the tumor during cryosurgery

AC C

(Chua et al., 2007).

Cryosurgical simulation of lung cancer build on systematic freezing time

is explored (Tarwidi, 2015). His studies show that temperature distribution and phase change interface position in tissue can be applied to increase destruction to tumor tissue and decrease the injury to normal tissue. Pennes bio-heat model (Pennes, 1948) has been globally applied by several authors (Kumar et al.,2018; Deng and Liu, 2005; Katiyar et al.,2017; Chua

5

ACCEPTED MANUSCRIPT

et al.,2007; Kumar et al. 2016) to answer the phase change heat problem in cryosurgery which is defined as ∂T = − 5 q + ρb wb cb (Tb − T ) + Qm ∂t

(1)

RI PT

cρ

where, c and ρ are the specific heat and the density of the tissue, cb and ρb

are the specific heat and density of blood, wb is the blood perfusion rate, q

is the heat flux, Qm is the metabolic heat generation, T is the temperature,

SC

Tb is the arterial blood temperature and t is the time.

The above Eqn.(1) is based on Fourier’s law which is defined as

M AN U

q(S, t) = −k 5 T (S, t)

(2)

where k is the thermal conductivity of the tissue, q(S, t) is the heat flux and T(S, t) is the temperature at position S = (x, y, z).

However, this model concludes the infinite speed of heat and mass propagation. Although in situations dealing with heat and mass transfer in ex-

TE D

tremely short durations of time, extremely high temperature and moisture gradient, extremely low temperatures and moisture towards absolute zero, or for microscale situations the wave nature of heat and mass propagation be-

EP

come dominant (Kumar et al.,2018). Cattaneo and Vernotte (1958) proposed modified Fourier’s and Fick’s models (Single phage model) generally called

AC C

as non-Fourier’s and non-Fick’s models. For such conditions, a constitutive equation which premises a relaxation time in the heat flux vector q. The single phage model is expressed by q(S, t + τq ) = −k 5 T (S, t),

where τq 6= 0 is the relaxation time due to heat flux. Using first order Taylor’s expansion, the SPL model reduces in the form   ∂q q + τq (S, t) = −k 5 T (S, t). ∂t 6

(3)

ACCEPTED MANUSCRIPT

There are several studies accessible in the literature (Deng and Liu, 2005; Katiyar et al., 2017; Kumar et al.,2016; Kumar et al., 2015; Ahmadikia and

RI PT

Moradi, 2012) where the bio-heat model with freezing is applied. The nonFourier effect of biological tissue of heat conduction during freezing has been used (Ahmadikia and Moradi, 2012). Further, the physical results provided by bio-heat model are sometimes different (Tzou, 1995). Thus to consider

SC

the thermal behavior that is not captured by the Fourier’s law, a new model by observing phase lag of temperature and phase lag of heat flux gradient is

and is expressed as

M AN U

established (Tzou, 1995). This model is known as dual phase lag model(DPL)

q(S, t + τq ) = −k 5 T (S, t + τT ),

where τq 6= 0 and τT 6= 0 are the relaxation time in heat propagation. Using

TE D

first order Taylor’s expansion, the reduced form of DPL model is     ∂q ∂T q + τq (S, t) = −k 5 T + τT (S, t). ∂t ∂t

(4)

Eq. (4) is called the dual phase lag constitutive relation. From Eqs. (1) and (4) we obtain the below equation:

(5)

AC C

EP

∂ 2T ∂T τq cρ 2 + (cρ + ρb wb cb τq ) = ∂t  ∂t  ∂T (S, t) + ρb wb cb (Tb − T ) + Qm k 52 T + τ T ∂t

The above equation is known as a dual phase lag(DPL) bio-heat equation. Many researchers (Majchrzak, 2010; Liu and Chen, 2009; Zhang, 2009;

Zhou et al., 2009) have applied Dual-phase lag bio-heat conduction model without phase change. DPL model is suggested to analyse the non-Fourier heat conduction (Antaki, 2005). Zhou et al.(2009) concluded that the thermal behavior of the DPL bio-heat model varies from the other bio-heat models. Similarly, Kumar et al.(2015) investigated the temperature behavior in 7

ACCEPTED MANUSCRIPT

biological tissues during hyperthermia treatment by using the DPL model. Ahmadikia and Moradi (2012) considered the freezing process in biologi-

RI PT

cal tissue with metabolic heat generation and blood flow by applying onedimensional DPL model. Recently, Kumar et al.(2018) have also used onedimensional DPL model to study the temperature distribution and moving layer thickness during cryosurgery of lung cancer.

SC

For a real study of lung cancer using cryosurgery, a two-dimensional math-

ematical bio-heat transfer model is required. In a series of papers Katiyar et al.(2007, 2016, 2017) studied two-dimensional bio-heat transfer models

M AN U

[Fourier model, SPL model, and DPL model] during cryosurgery of lung cancer when tumor tissue is cooled by a cryoprobe with constant temperature. They solved these problems by enthalpy based finite difference method. A mathematical model during cryosurgery of lung cancer at stage1 is a twodimensional boundary value problem in presence of lagging time (relaxation

TE D

time τq and τT ) and a DPL bio-heat transfer model is appropriate. In stage 2 and 3, it will be a moving boundary problem in absence of lagging time and a Fourier bio-heat transfer model is appropriate. To the best of our knowledge, this type of mathematical model is not considered yet.

EP

The aim of our study is to determine how the freezing appears in three stages when lung tissue cooled by a cryoprobe by imposing on it a constant

AC C

temperature or a constant heat flux or a constant heat transfer coefficient, i.e. by imposing on it the boundary condition of first kind or second kind or third kind. We see how the moving layer thickness and temperature distribution changes as the boundary condition changes. We studied the two-dimensional process of freezing in three various regions: unfrozen region, mushy region, and frozen region. The frequency of heat transfer decides the frequency of tissue cooling. As the aim is to damage the tissues, it is necessary to obtain

8

ACCEPTED MANUSCRIPT

the rapid cooling rate (Paola, 2015). In the tumor, the freezing interface accelerates as it moves from tumor-tissue into the normal lung tissue (Bischof et

RI PT

al., 1992). Cryosurgical process is simulated mathematically as the equations of heat transfer in liquid and solid phases, where the mushy region (interface region) between two phases is subject to Interface condition (Kotova et

AC C

EP

TE D

M AN U

SC

al.,2016).

Figure 1: Schematic of freezing a tumor within the lung

9

ACCEPTED MANUSCRIPT

2. Formulation of the problem We observe a lung tissue occupied in the domain ω = [(x, y) : 0 ≤ x ≤ l, 0 ≤ y ≤ l]

RI PT

in which a tumor tissue of dimension ωt = [(x, y) : 0 ≤ x ≤ a, 0 ≤ y ≤ a] is

fixed. A cryoprobe is positioned at x = 0, dcm ≤ y ≤ ecm. The lung tissue freezed by a cryoprobe by keeping it at a persistent temperature or a persis-

tent heat flux or a persistent heat transfer coefficient. The freezing begins

SC

and the complete procedure of freezing is controlled in three stages. In stage

I, the cryoprobe surface is cooled from the initial temperature T0 (370 C) to

M AN U

liquidus temperature Tl (−10 C). In this stage, unfrozen region is formed. Further, in stage II the cryoprobe is perpetually cooled from the liquidus temperature Tl to solidus temperature Ts (−80 C). Throughout this stage, unfrozen region and mushy region co-exist at the same time. In stage III, when the temperature of the cryoprobe is perpetually reduced with time the freezing begins from the origin and spreads in the positive x-direction. In

TE D

this stage, unfrozen region, mushy region, and frozen region are set up. The present problem can be described in the given Fig. 1. The temperature distribution of the tissue is above −10 C in stage I and

EP

during freezing, the relaxation time plays the crucial role and due to this nonFourier model(DPL) is considered (Kumar et al., 2016). While the Fourier

AC C

model is used in Stage II and III. From our calculations, we see that the effect of relaxation time is negligible in the mushy and frozen region and present only in unfrozen region. The following presumptions are examined to solve the two-dimensional mathematical model:-

(i) Non-Fourier heat conduction law followed by Heat conduction (Budman et al.,1995). (ii) When tissue is unfrozen heat source seems due to blood perfusion and 10

ACCEPTED MANUSCRIPT

metabolism (Deng and Liu, 2005; Frayn, 1996). (iii) Flawless property of tissues is manipulated with liquidus(upper phase)

RI PT

and solidus(lower phase) temperature as −10 C and −80 C respectively (Rabin and Shitzer, 1995).

(iv) The primary(initial) temperature of the tissue is examined as an arterial temperature (370 C).

from liquid to solid (Lee and Bastacky, 1995).

SC

(v) Thermophysical properties vary at the point of complete phase change

(vi) As the solid continuously forms within the intermediate zone, the latent

M AN U

heat of fusion is released accordingly. The heat of fusion can be treated as internal heat generation A within the mushy region (Kumar et al., 2018; Chua et al, 2007; Yadav et al.,2016) is defined as: A=ρ∗L

∂fs ∂t

(6)

TE D

where fs is solid fraction present in mushy region.

(vii) The rate of change of solid fraction with respect to time in the mushy region provides a heat generation effect. An exact mathematical expression for the solid fraction cannot be given but the approximate mathematical

EP

model is proposed by Gupta (2003) for the numerical solution of the problem

AC C

and defined as

fs =

fsu (Tl − Tm ) − fsu1 (Ts − Tm ) Tl − Ts

(7)

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. 11

ACCEPTED MANUSCRIPT

(viii) Perfect thermal contact is assumed between the probe surface and the adjacent medium.

RI PT

The mathematical model governing the action of heat transfer comes out in three stages are as follows: Stage I

  3  ∂ 2 Tu ∂ 2 Tu ∂ Tu ∂ 3 Tu = ku + τT ku (8) + + ∂x2 ∂y 2 ∂t1 ∂x2 ∂t1 ∂y 2   ∂Tu + Qm , 0 < x, y < l, t1 > 0 ρb cb wb Tb − Tu − τq ∂t1





M AN U

ρ u cu

∂Tu ∂ 2 Tu + τq 2 ∂t1 ∂t1

SC



Initial and boundary condition are given as:-

EP

Stage II

TE D

Tu (x, y, 0) = T0 , ∂Tu (x, y, 0) = 0, ∂t ∂Tu (l, y, t1 ) = 0, ∂x ∂Tu (x, l, t1 ) = 0. ∂y

 ∂ 2 Tm ∂ 2 Tm + + A, ∂x2 ∂y 2   2 ∂Tu ∂ Tu ∂ 2 Tu , ρ u cu = ku + ∂t2 ∂x2 ∂y 2

AC C

∂Tm ρm cm = km ∂t2

(9) (10) (11) (12)



0 < x, y < s2 , t2 > t1

(13)

s2 < x, y < l

(14)

Initial and boundary condition are given as:Tm (x, y, t1 ∗) = Tu (x, y, t1 ∗), ∂Tm (l, y, t1 ) = 0, ∂x ∂Tm (x, l, t1 ) = 0, ∂y Tu (s2 , t2 ) = Tm (s2 , t2 ) = Tl 12

(15) (16) (17) (18)

ACCEPTED MANUSCRIPT

Interface condition is given by: ∂Tm ∂Tu − ku = ρ ∗ Lvn . ∂x ∂x

where t2 = t1 − t1 ∗ Stage III

∂Tm = km ρm cm ∂t3

∂ 2 Tf ∂ 2 Tf + ∂x2 ∂x2 

∂Tu = ku ρ u cu ∂t3

 ,

∂ 2 Tm ∂ 2 Tm + ∂x2 ∂y 2 

SC



0 < x, y < s1 , t3 > t2 

M AN U

∂Tf ρ f cf = kf ∂t3

(19)

RI PT

km

∂ 2 Tu ∂ 2 Tu + ∂x2 ∂y 2

+ A,

s1 < x, y < s2

(20)

(21)



,

s2 < x, y < l

(22)

Initial and boundary condition are given as:-

AC C

EP

TE D

Tf (x, y, t2 ∗) = Tm (x, y, t2 ∗) = Tu (x, y, t2 ∗), ∂Tf (l, y, t1 ) = 0, ∂x ∂Tf (x, l, t1 ) = 0, ∂y Tf (s1 , t3 ) = Tm (s1 , t3 ) = Ts , Tu (s2 , t3 ) = Tm (s2 , t3 ) = Tl

(23) (24) (25) (26) (27)

Interface condition are given by: ∂Tf ∂Tm − km = ρ ∗ Lvn ∂x ∂x ∂Tm ∂Tu km − ku = ρ ∗ Lvn . ∂x ∂x kf

(28) (29)

where t3 = t2 − t2 ∗ and subscripts u, m, f and b denotes unfrozen, mushy, frozen and blood, T is the temperature, t is the time, Qm is the metabolic heat generation, τq and τT are the relaxation times. 13

ACCEPTED MANUSCRIPT

2.1. General Boundary conditions: Due to existence of cryoprobe at surface x = 0 the tissue is cooled under

Ao

∂T (0, y, t) + Bo T (0, y, t) = f (y, t) ∂x

SC

I kind:

RI PT

generalised boundary condition is as follows:

Ao = 0, Bo = 1, f (t) = Tw

M AN U

II kind:

Ao = −k, Bo = 0, f (t) = qw

TE D

III kind:

Ao = −k, Bo = −h, f (t) = −hTs

(30)

(31)

(32)

(33)

as

EP

Continuity of heat flux and temperature at lung-tumor boundary is given

AC C

kt

∂Tt ∂Tl (x, y, t) = kl (x, y, t) ∂x ∂x Tt = Tl

(34)

where, subscripts t and l stands for tumor and normal lung tissue, respectively.

Thermal properties of tumor and healthy lung tissues are different (Bischof

14

ACCEPTED MANUSCRIPT

et al. 1992; Deng and Liu, 2005). Thermophysical properties of healthy lung

RI PT

and tumor tissue (Kumar and Katiayar, 2016) are taken as     kf , T < T l    k = 1 (kf + ku ), Tl ≤ T ≤ Ts , 2      ku , T > Ts

M AN U

SC

    ρf , T < T l    ρ = 1 (ρf + ρu ), Tl ≤ T ≤ Ts , 2      ρ u , T > T s

TE D

    cf , T < T l    c = 1 (cf + cu ), Tl ≤ T ≤ Ts . 2      c u , T > T s

3. Dimensionless Analysis

EP

We define some dimensionless variable and similarity criteria as follows y s s1 s2 To − Ts x To − Tw , Y = , S = , S1 = , S2 = , θs = , θw = , l l l l l To − Tl To − Tl To − Tr To − Tb kt kτq kτT θr = , θb = , F0 = , F oq = , F oT = , 2 2 To − Tl To − Tl ρ ∗ cl ρ ∗ cl ρ ∗ cl2 cb w b l 2 Qm l2 ku km ku km Pf2 = , Pm = , kum = , kmf = , a1 = , a2 = , k k(To − Tl ) km kf ρ ∗ c1 ρ ∗ c2     kf To − Tl Tl − Ts qw l a3 = , Ste2 = cm , Ste3 = cf , Ki = , ρ ∗ c3 L L (To − Tl )ku hl sj (ti ) Bi = , λj (F oi ) = , j = 1, 2, i = 2, 3, ku l Ao (To − Tl ) A1 = , B1 = Bo (To − Tl ), g(Fo ) = f (t) − Bo (To − Tl ). l

AC C

X=

15

ACCEPTED MANUSCRIPT

Using these dimensionless variable our problem can be reduce in the following

Stage I

0 < X, Y < 1(35)

M AN U

Initial and boundary condition are given as:-

SC

 2 
∂θu ∂ 2 θu ∂ θu ∂ 2 θu 2 F oq = + 1 + Pf F o q + ∂F o21 ∂F o1 ∂X 2 ∂Y 2  2  ∂ θu ∂ 2 θu ∂ + Pf2 (θb − θu ) + Pm , +F oT + 2 2 ∂F0 ∂X ∂Y

RI PT

form:-

(36) (37) (38) (39)

TE D

θu (X, Y, 0) = 0 ∂θu (X, Y, 0) = 0 ∂F0 ∂θu (1, Y, F o1 ) = 0 ∂X ∂θu (X, 1, F o1 ) = 0 ∂Y

where

AC C

Stage II

To − Tu To − Tl

EP

θu =

 1+

 1 ∂θm ∂ 2 θm ∂ 2 θm + , 0 < X, Y < λ2 (F o2 ) = Ste2 ∂F o2 ∂X 2 ∂Y 2  2  ∂θu ∂ θu ∂ 2 θu + , λ2 (F o2 ) < X, Y < 1 = a12 ∂F o2 ∂X 2 ∂Y 2

16

(40) (41)

ACCEPTED MANUSCRIPT

Initial and boundary condition are given as:-

SC

Interface condition is given by:

(43)

(44)

(45)

(46)

M AN U

∂θm ∂θu −1 ∂X − kum = , X = λ2 (F o2 ) ∂X ∂X Ste2 ∂F o2 where θm =

(42)

RI PT

θm (X, Y, F o1 ∗) = θu (X, Y, F o1 ∗), ∂θm (1, Y, F o1 ) = 0, ∂X ∂θm (X, 1, F o1 ) = 0, ∂Y θm (λ2 (F o2 ), F o2 ) = θu (λ2 (F o2 ), F o2 ) = 0

Tl − Tm Tl − Tu a1 , θu = , F o∗1 = F o2 − F o1 , a12 = . To − Tl To − Tl a2

TE D

Stage III

AC C

EP

∂ 2 θf ∂ 2 θf ∂θf = + , 0 < X, Y < λ1 (F o3 ) (47) ∂F o ∂X 2 ∂Y 2   2 3 2   1 ∂θm ∂ θm ∂ θm 1− + = a23 , λ1 (F o3 ) < X, Y < λ2 (F o3 ) (48) Ste3 ∂F o3 ∂X 2 ∂Y 2   2 ∂θu ∂ θu ∂ 2 θu = a13 + , λ2 (F o3 ) < X, Y < 1 (49) ∂F o3 ∂X 2 ∂Y 2 Initial and boundary condition are given as:θf (X, Y, F o2 ∗) = θm (X, Y, F o2 ∗) = θu (X, Y, F o2 ∗), ∂θf (1, Y, F o1 ) = 0, ∂X ∂θf (X, 1, F o1 ) = 0, ∂Y θf (λ1 (F o3 ), F o3 ) = θm (λ1 (F o3 ), F o3 ) = θs , θm (λ2 (F o3 ), F o3 ) = θu (λ2 (F o3 ), F o3 ) = 1 17

(50) (51) (52) (53) (54)

ACCEPTED MANUSCRIPT

Interface condition are given by:

where

(56)

Tm − Ts Tu − Ts a2 a1 Tf − Ts , θm = , θu = , F o∗2 = F o3 − F o2 , a23 = , a13 = . Tl − Ts Tl − Ts Tl − Ts a3 a3

SC

θf =

(55)

RI PT

∂θf ∂θm −1 ∂X − kmf = , X = λ1 (F o3 ), ∂X ∂X Ste3 ∂F o3 ∂θm ∂θu −1 ∂X − kum = , X = λ2 (F o3 ). ∂X ∂X Ste3 ∂F o3

3.1. Boundary condition in dimensionless form:

M AN U

Now also reducing boundary conditions in dimensionless form, Eqns.(11),(12),(13) becomes:A1

∂θ (0, Y, Fo ) + B1 θ(0, Y, Fo ) = g(Fo ) ∂X

I kind:

(58)

∂θ (0, Y, Fo ) = −Ki ∂X

(59)

∂θ (0, Y, Fo ) = −Bi(θw − θs ). ∂X

(60)

TE D

AC C

III kind:

θ(0, Y, Fo ) = θw

EP

II kind:

(57)

4. Solution of the problem: Using finite differences and after simple algebraic computation, our prob-

lem reduces as follows:

18

ACCEPTED MANUSCRIPT

Boundary condition of I kind: For Stage I,
dΘu
d2 Θu 2 2 + I + P F o − F o U + P I − U Θu q T 1 1 f f dF o21 dF o1 θw θr θr = (Pf2 θb + Pm )d + 2 d1 + 2 d2 + 2 d3 h h k

RI PT

F oq

(61)

EP

TE D

M AN U

SC

where Θu , U1 , d, d1 , d2 and d3 are given by:     −2, 1, 0, .......0, 0, 0   θu1      1, −2, 1, .......0, 0, 0   2    θu        0, 1, −2, .......0, 0, 0  3    θu  1 .. .. ..   ..   , , U1 = 2 Θu =  .   ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   ..    . . . .    .., . . . , ..; . . . ; .., . . . , ..   ..     .      0, 0, 0, ......., 1, −2, 0 N   θu N ×1 0, 0, 0, ......., 0, 2, −2        N ×N 1 1 0 1                 1 0 0 0                 1 0 0 0        d = . , d1 =  .  , d2 =  .  , d3 =  .  ..  . . . . . . .          ..   ..   ..   ..  . . . .         1 0 1 1

AC C

N ×1

N ×1

N ×1

For Stage II,   1 dΘm θw θr θr 1+ − U2 Θm = 2 d1 + 2 d2 + 2 d3 , Ste2 dF o2 h h k dΘu − a12 U3 Θu = 0, dF o2 Θm (F o∗1 ) = Θu (F o∗1 )

19

N ×1

0 < X, Y < λ2 (F o(62) 2 ), λ2 (F o2 ) < X, Y <(63) 1 (64)

ACCEPTED MANUSCRIPT

For Stage III,

TE D

M AN U

SC

RI PT

where Θm , Θu , U2 and U3 are given as follows:     −2, 1, 0, .......0, 0, 0   θ1    m   1, −2, 1, .......0, 0, 0   2     θm       0, 1, −2, .......0, 0, 0   3     θm  1 .. .. ..   ..  , U = , Θm =   ., . . . , .; . . . ; ., . . . , .  2  ..   (h2 + k 2 )   .   .. . . .    ., . . . , ..; . . . ; .., . . . , ..   ..     .        0, 0, 0, ......., 1, −2, 0 ¯ −1 N   θm ¯ −1×1 N 0, 0, 0, ......., 0, 1, −2 ¯ ¯ N −1×N −1    −2, 1, 0, .......0, 0, 0 ¯   θuN +1      1, −2, 1, .......0, 0, 0   N¯ +2    θu        0, 1, −2, .......0, 0, 0  N¯ +3    θu  1 .. .. ..   ..  , U = Θu =   ., . . . , .; . . . ; ., . . . , .  3  ..   (h2 + k 2 )   .   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −2, 0 N   θu ¯ ×1 N −N 0, 0, 0, ......., 0, 2, −2 ¯ ¯

N −N ×N −N

AC C

EP

dΘf θw θs − U4 Θf = 2 d0 + 2 d01 , 0 < X, Y < λ1 (F o3 )(65) 2 2 dF o3 (h + k ) (h + k )   1 dΘm 1 1− − a23 U5 Θm = 2 d0 , λ1 (F o3 ) < X, Y < λ2 (F o3 )(66) Ste3 dF o3 (h + k 2 ) 2 dΘu 1 − a13 U6 Θu = 2 d0 , λ2 (F o3 ) < X, Y < 1(67) dF o3 (h + k 2 ) 3 Θf (F o∗2 ) = Θm (F o∗2 ) = Θu (F o∗2 )(68)

20

ACCEPTED MANUSCRIPT

where,  

θf1



−2, 1, 0, .......0, 0, 0



     1, −2, 1, .......0, 0, 0  2   θf    0, 1, −2, .......0, 0, 0    θf3 1 .. .. ..   .. , U = ,  ., . . . , .; . . . ; ., . . . , .  4 ..  (h2 + k 2 )  .  .. . . .  ., . . . , ..; . . . ; .., . . . , ..  ..   .   0, 0, 0, ......., 1, −2, 0   θfp¯−1 p¯−1×1 0, 0, 0, ......., 0, 1, −2 p¯−1ׯ p−1     −2, 1, 0, .......0, 0, 0 1   θm      1, −2, 1, .......0, 0, 0   2     θm       0, 1, −2, .......0, 0, 0   3     θm  1 .. .. ..   ..   , U5 = 2 Θm =  .   ., . . . , .; . . . ; ., . . . , .  ,  (h + k 2 )   ..   .. .. .. ..      .; . . . ; ., . . . , . ., . . . ,  ..     .      0, 0, 0, ......., 1, −2, 0 q¯−1   θm q¯−1×1 0, 0, 0, ......., 0, 1, −2 q×q     −2, 1, 0, .......0, 0, 0 1   θu      1, −2, 1, .......0, 0, 0   2     θu       0, 1, −2, .......0, 0, 0   3     θu  1 .. .. ..   ..  Θu =  , U =   , ., . . . , .; . . . ; ., . . . , . 6  ..   (h2 + k 2 )   .   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −2, 0   θur¯−1 r¯−1×1 0, 0, 0, ......., 0, 2, −2

SC

RI PT

            

AC C

EP

TE D

M AN U

      Θf =       

r×r

(69)

21

ACCEPTED MANUSCRIPT

p×1

  0     0     0 0  , d2 =   ..  .    ..  .   1

q×1

  1     0     0 0  , d3 =   ..  .    ..  .   0

.

RI PT

p×1

  1     0     0 0  , d1 =   ..  .    ..  .   0

r×1

SC

  1     1     1 0  d =  ..  .    ..  .   1

Boundary condition of II kind:

M AN U

For Stage I,


dΘu
d2 Θu 2 0 2 0 + 1 + P F o − F o U + P I − U Θu F oq q T f 1 f 1 dF o21 dF o1 Ki Ki Ki Ki d1 − d2 + d1 − d2 = (Pf2 θb + Pm )d + h h k k where Θu and U10 are given by:





,

AC C

EP

θu1



TE D

−1, 1, 0, .......0, 0, 0        1, −1, 1, .......0, 0, 0   2    θu        0, 1, −1, .......0, 0, 0  3    θu  1 .. .. ..   .. 0   , U1 = 2 Θu =  .   ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   ..   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −1, 0 N   θu N ×1 0, 0, 0, ......., 0, 2, −1 

(70)

N ×N

For Stage II,   1 dΘm Ki Ki Ki Ki 1+ − U20 Θm = d1 − d2 + d1 − d2 , 0 < X, Y < λ2 (F o2 ),(71) Ste2 dF o2 h h k k dΘu − a12 U30 Θu = 0, λ2 (F o2 ) < X, Y < 1(72) dF o2 Θm (F o∗1 ) = Θu (F o∗1 )(73)

22

ACCEPTED MANUSCRIPT

For Stage III,

TE D

M AN U

SC

RI PT

where Θm , Θu , U20 and U30 are given as follows:     −1, 1, 0, .......0, 0, 0   θ1    m   1, −1, 1, .......0, 0, 0   2     θm        0, 1, −1, .......0, 0, 0  3     θm  1 .. .. ..   .. 0  = , , U Θm =   ., . . . , .; . . . ; ., . . . , .  2  ..   (h2 + k 2 )   .   .. . . .    ., . . . , ..; . . . ; .., . . . , ..   ..     .        0, 0, 0, ......., 1, −1, 0 ¯ −1 N   θm ¯ −1×1 N 0, 0, 0, ......., 0, 1, −1 ¯ ¯  N −1×N −1   −1, 1, 0, .......0, 0, 0 ¯   θuN +1      1, −1, 1, .......0, 0, 0   N¯ +2    θu       0, 1, −1, .......0, 0, 0   N¯ +3    θu  1 .. .. ..   .. 0   , U3 = 2 Θu =  .   ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   ..   .. .. .. ..      ., . . . , .; . . . ; ., . . . , .  ..     .      0, 0, 0, ......., 1, −1, 0 N   θu ¯ ×1 N −N 0, 0, 0, ......., 0, 2, −1 ¯ ¯

N −N ×N −N

1 1− Ste3



dΘm 1 − a23 U50 Θm = 2 d0 , dF o3 (h + k 2 ) 2 1 dΘu − a13 U60 Θu = 2 d0 , dF o3 (h + k 2 ) 3 Θm (F o∗2 ) = Θu (F o∗2 ) = Θf (F o∗2 )

AC C



EP

Ki Ki Ki Ki dΘf − U40 Θf = d1 − d2 + d1 − d2 , 0 < X, Y < λ1 (F o3 ) (74) dF o3 h h k k

23

λ1 (F o3 ) < X, Y < λ2 (F(75) o3 ) λ2 (F o3 ) < X, Y < 1

(76) (77)

ACCEPTED MANUSCRIPT

where,  θf2

            

θf3 .. . .. . θfp¯−1



     1, −1, 1, .......0, 0, 0       0, 1, −1, .......0, 0, 0    1 .. .. ..   .. , U40 = 2  ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   .. .. .. ..   ., . . . , .; . . . ; ., . . . , .      0, 0, 0, ......., 1, −1, 0   p¯−1×1 0, 0, 0, ......., 0, 1, −1

RI PT

θf1

−1, 1, 0, .......0, 0, 0

,

SC

      Θf =       



M AN U





p¯−1ׯ p−1

(78)



−1, 1, 0, .......0, 0, 0         1, −1, 1, .......0, 0, 0  2     θm       0, 1, −1, .......0, 0, 0   3     θm  1 .. .. ..   .. 0   , U5 = 2 Θm =  .   ., . . . , .; . . . ; ., . . . , .  ,  (h + k 2 )   ..   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −1, 0 q¯−1   θm q¯−1×1 0, 0, 0, ......., 0, 1, −1 q×q     −1, 1, 0, .......0, 0, 0   θu1      1, −1, 1, .......0, 0, 0   2     θu       0, 1, −1, .......0, 0, 0   3     θu  1 .. .. ..   .. 0  Θu =  , U =  ., . . . , .; . . . ; ., . . . , .  6  ..   (h2 + k 2 )   .   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −1, 0   θur¯−1 r¯−1×1 0, 0, 0, ......., 0, 2, −1 1 θm



AC C

EP

TE D



r×r

24

ACCEPTED MANUSCRIPT

Boundary condition of III kind: For Stage I,

RI PT



d2 Θu 2 00 dΘu 2 00 + 1 + P F o − F o U + P I − U Θu = (Pf2 θb + Pm )d q T f 1 f 1 dF o21 dF o1 Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) + d1 − d2 + d1 − d2(79) h h k k

where Θu and U100 are given by: 

SC

F oq



−1, 1, 0, .......0, 0, 0         1, −1, 1, .......0, 0, 0  2    θu       0, 1, −1, .......0, 0, 0   3    θu  1 .. .. ..   .. 00   Θu =  .  , U1 = 2  ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   ..   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −1, 0   θuN N ×1 0, 0, 0, ......., 0, 2, −1 θu1



,

N ×N

TE D

M AN U



EP

For Stage II,   1 dΘm Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) 1+ − U200 Θm = d1 − d2 + d1 Ste2 dF o2 h h k Bi (θw − θr ) − d2 , 0 < X, Y < λ2 (F o2 ),(80) k

AC C

dΘu − a12 U300 Θu = 0, dF o2 Θm (F o∗1 ) = Θu (F o∗1 )

25

λ2 (F o2 ) < X, Y < 1

(81) (82)

ACCEPTED MANUSCRIPT

For Stage III,

TE D

M AN U

SC

RI PT

where Θm , Θu , U200 and U300 are given as follows:     −1, 1, 0, .......0, 0, 0   θ1    m   1, −1, 1, .......0, 0, 0   2     θm        0, 1, −1, .......0, 0, 0  3     θm  1 .. .. ..   .. 00  = , , U Θm =   ., . . . , .; . . . ; ., . . . , .  2  ..   (h2 + k 2 )   .   .. . . .    ., . . . , ..; . . . ; .., . . . , ..   ..     .        0, 0, 0, ......., 1, −1, 0 ¯ −1 N   θm ¯ −1×1 N 0, 0, 0, ......., 0, 1, −1 ¯ ¯  N −1×N −1   −1, 1, 0, .......0, 0, 0 ¯   θuN +1      1, −1, 1, .......0, 0, 0   N¯ +2    θu       0, 1, −1, .......0, 0, 0   N¯ +3    θu  1 .. .. ..   .. 00   , U3 = 2 Θu =  .   ., . . . , .; . . . ; ., . . . , .   (h + k 2 )   ..   .. .. .. ..      ., . . . , .; . . . ; ., . . . , .  ..     .      0, 0, 0, ......., 1, −1, 0 N   θu ¯ ×1 N −N 0, 0, 0, ......., 0, 2, −1 ¯ ¯

N −N ×N −N

AC C

EP

Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) dΘf − U400 Θf = d1 − d2 + d1 dF o3 h h k Bi (θw − θr ) − d2 , 0 < X, Y < λ1 (F o3 )(83) k   1 dΘm 1 1− − a23 U500 Θm = 2 d0 , λ1 (F o3 ) < X, Y < λ2 (F o3 )(84) Ste3 dF o3 (h + k 2 ) 2 dΘu 1 − a13 U600 Θu = 2 d0 , λ2 (F o3 ) < X, Y < 1(85) dF o3 (h + k 2 ) 3 Θm (F o∗2 ) = Θu (F o∗2 ) = Θf (F o∗2 )(86)

26

ACCEPTED MANUSCRIPT

where,  

θf1



−1, 1, 0, .......0, 0, 0



     1, −1, 1, .......0, 0, 0  2   θf    0, 1, −1, .......0, 0, 0    θf3 1 .. .. ..   .. 00 = , U ,  ., . . . , .; . . . ; ., . . . , .  4 ..  (h2 + k 2 )  .  .. . . .  ., . . . , ..; . . . ; .., . . . , ..  ..   .   0, 0, 0, ......., 1, −1, 0   θfp¯−1 p¯−1×1 0, 0, 0, ......., 0, 1, −1 p¯−1ׯ p−1     −1, 1, 0, .......0, 0, 0 1   θm      1, −1, 1, .......0, 0, 0   2     θm       0, 1, −1, .......0, 0, 0   3     θm  1 .. .. ..   .. 00   , U5 = 2 Θm =  .   ., . . . , .; . . . ; ., . . . , .  ,  (h + k 2 )   ..   .. .. .. ..      .; . . . ; ., . . . , . ., . . . ,  ..     .      0, 0, 0, ......., 1, −1, 0 q¯−1   θm q¯−1×1 0, 0, 0, ......., 0, 1, −1 q×q     −1, 1, 0, .......0, 0, 0 1   θu      1, −1, 1, .......0, 0, 0   2     θu       0, 1, −1, .......0, 0, 0   3     θu  1 .. .. ..   .. 00  Θu =  , U =   . ., . . . , .; . . . ; ., . . . , . 6  ..   (h2 + k 2 )   .   .. .. .. ..     ., . . . , .; . . . ; ., . . . , .   ..     .      0, 0, 0, ......., 1, −1, 0   θur¯−1 r¯−1×1 0, 0, 0, ......., 0, 2, −1

SC

RI PT

            

AC C

EP

TE D

M AN U

      Θf =       

r×r

5. Modified Legendre wavelet Galerkin Method : We are using Legendre wavelet Galerkin method to solve our problem in all the stages with generalised boundary condition are as follows:

27

ACCEPTED MANUSCRIPT

Boundary condition of I kind 5.1. Stage 1 d2 Θu dF o1 2

is approximated by

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

RI PT

Let us assume that the unknown function

(87)

0 ··· 0



  ψ · · · 0  .. . . ..  . . .  0 ··· ψ

M AN U

of order 22(k−1) M 2 × J defined by  ψ   0 ¯ ψ=  .. .  0

SC

¯ o1 ) is a matrix where C is unknown matrix of order N × 22(k−1) M 2 and ψ(F

where ψ(F o) = [ψ1,0 , ψ1,1 , . . . , ψ1,M −1 , . . . , ψ2k−1 ,0 , ψ2k−1 ,1 , . . . ψ2k−1 ,M −1 ]0 and

TE D

the elements of ψ(F o) is defined by  p
  (m + 1/2)2k/2 Pm 2k F o1 − n ˆ , ψn,m (F o1 ) =  0 ,

n ˆ −1 2k

≤ F o1 ≤

n ˆ +1 2k

otherwise.

EP

where k = 1, 2, 3..., n = 1, 2, . . . , 2k−1 , n ˆ = 2n − 1 and m is the order of Legendre polynomial (Yadav et al., 2014). Pm (F o1 ) is denoted by Legendre

AC C

polynomial of 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

28

ACCEPTED MANUSCRIPT

The operational matrix of integration of ψ defined in (Kumar et al., 2018;

ZF o1 ψ(t)dt = P ψ(F o1 )

RI PT

Razzaghi and Yousefi, 2001; Kumar et al., 2016) and given by

(88)

0

where P if is the operational matrix of integration of order 2k−1 M × 1, given by O O ··· O

  L O · · · O   0 L · · · O  .. .. . . ..  . . . .   0 0 · · · O  0 0 ··· O 

M AN U

L   0   0 1  P = k  2  ...   0  0



SC



AC C

EP

TE D

 2 0 ··· 0     0 0 · · · 0   where O and L are M × M matrices given by O =  . . . ..  . . . . . . .   0 0 ··· 0 and   0 0 ··· 0 0 0 1 √13    −1  1 √ √  3 0  0 ··· 0 0 0 15     √1  0 √−1  0 · · · 0 0 0 15 35     −1 √  0 0 0 ··· 0 0 0 L= 35     .. .. .. .. .. .. .. ...   . . . . . . .     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.

29

ACCEPTED MANUSCRIPT

dΘu ∼ dΘu (0) + P¯ ψ(F o1 ) = dF o1 dF o1 dΘu ∼ ¯ = C P ψ(F o1 ) dF o1 where P¯ is defined as P

0 ··· 0



also

M AN U

    0 P ··· 0 ¯   P =. . . ..  . . . . . . .   0 0 ··· P

(89)

SC



RI PT

Integrating (85) over 0 to F o1 , we obtained

ZF o ¯ o)dF o = P¯ ψ(F ¯ o) ψ(F 0

Integrating again (85) over 0 to F o1 , we obtained

TE D

¯ o1 ) Θu ∼ = Θu (0) + C P¯ 2 ψ(F ¯ o1 ) Θu ∼ = C P¯ 2 ψ(F

(90)

EP

After solving Eqns(61),(85),(87) and (88), we obtained

AC C



F oq C ψ¯ + I + Pf2 F oq I − F oT U C P¯ ψ¯ + Pf2 I − U C P¯ 2 ψ¯ θr θr θw = (Pf2 θb + Pm )H ψ¯ + 2 H1 ψ¯ + 2 H2 ψ¯ + 2 H3 ψ¯ h h k

(91)

¯ d1 = H1 ψ, ¯ d2 = H2 ψ, ¯ d3 = H3 ψ¯ and F od = F ψ¯ where d = H ψ,

F oq C + I + Pf2 F oq I − F oT U C P¯ + Pf2 I − U C P¯ 2 θw θr θr = (Pf2 θb + Pm )H + 2 H1 + 2 H2 + 2 H3 h h k

30

(92)

ACCEPTED MANUSCRIPT

5.2. Stage 2 Let us assume that the unknown functions

dΘm dF o2

and

dΘu dF o2

are approxi-

RI PT

mated by dΘm ∼ ¯ = D1 ψ, dF o2 dΘu ∼ ¯ = D2 ψ, dF o2

(93)

(94)

SC

Integrating (73,74) over F o∗1 to F o2 , we obtained

(95)

¯ o1 ∗) Θu ∼ = D2 P¯ ψ¯ − D2 P¯ ψ(F

(96)

M AN U

¯ o1 ∗), Θm ∼ = D1 P¯ ψ¯ − D1 P¯ ψ(F

After solving Eqns (62),(91) and (93), we obtained   1 ¯ o1 ∗)F − D1 U2 P¯ = θw H1 + θr H2 + θr H3 + U2 d1 F(97) D1 I + D1 U2 P¯ ψ(F 1+ Ste2 h2 h2 k2 After solving Eqns (63),(92) and (94), we obtained

5.3. Stage 3

TE D

¯ o1 ∗)F − a12 D2 U3 P¯ = D2 I + a12 D2 U3 P¯ ψ(F

Let us assume that the unknown functions

AC C

EP

proximated by

(h2

1 dF + k2)

dΘf , dΘm dF o3 dF o3

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

and

(98)

dΘu dF o3

are ap-

(99) (100) (101)

Integrating (97-99) over F o∗2 to F o3 , we obtained Θf ∼ = E1 P¯ ψ − E1 P¯ ψ¯ (F o∗2 ) ,

(102)

Θm ∼ = E2 P¯ ψ − E2 P¯ ψ¯ (F o∗2 ) ,

(103)

Θu ∼ = E3 P¯ ψ − E3 P¯ ψ¯ (F o∗2 ) .

(104)

31

ACCEPTED MANUSCRIPT

After solving Eqns (65), (97) and (100), we obtained ¯ o∗ )F − E1 U4 P¯ = E1 I + E1 U4 P¯ ψ(F 2

θw θs 0 H + H0 (h2 + k 2 ) (h2 + k 2 ) 1

(105)

RI PT

After solving Eqns (66), (98) and (101), we obtained   1 θw ¯ o∗ )F − a23 E2 U5 P¯ = 1− E2 I + a23 E2 U5 P¯ ψ(F H20 (106) 2 2 2 Ste3 (h + k ) After solving Eqns (67), (99) and (102), we obtained

Boundary condition of II kind

θw H0 + k2) 3

M AN U

5.4. Stage 1

(h2

After solving Eqns (69), (85), (87) and (88), we obtained

F oq C + I + Pf2 F oq I − F oT U 0 C P¯ + Pf2 I − U 0 C P¯ 2 Ki Ki Ki Ki H1 − H2 + H1 − H2 = (Pf2 θb + Pm )H + h h k k 5.5. Stage 2

(108)

After solving Eqns (70),(91) and (93), we obtained  1 ¯ o1 ∗)F − D1 U 0 P¯ = Ki H1 − Ki H2 + Ki H1 − Ki H(109) D1 I + D1 U20 P¯ ψ(F 1+ 2 2 Ste2 h h k k

TE D



(107)

SC

¯ o∗ )F − a13 E3 U6 P¯ = E3 I + a13 E3 U6 P¯ ψ(F 2

After solving Eqns (71),(92) and (94), we obtained

EP

¯ o1 ∗)F − a12 D2 U 0 P¯ = D2 I + a12 D2 U30 P¯ ψ(F 3

(h2

1 dF + k2)

(110)

5.6. Stage 3

AC C

After solving Eqns (73),(97) and (100), we obtained

¯ o∗ )F − E1 U 0 P¯ = Ki H1 − Ki H2 + Ki H1 − Ki H2 (111) E1 I + E1 U40 P¯ ψ(F 2 4 h h k k

After solving Eqns (74),(98) and (101), we obtained   1 1 ¯ o∗ )F − a23 E2 U 0 P¯ = 1− E2 I + a23 E2 U50 P¯ ψ(F H 0 (112) 2 5 2 Ste3 (h + k 2 ) 2

After solving Eqns (75),(99) and (102), we obtained ¯ o∗ )F − a13 E3 U 0 P¯ = E3 I + a13 E3 U60 P¯ ψ(F 2 6 32

1 H0 (h2 + k 2 ) 3

(113)

ACCEPTED MANUSCRIPT

Boundary condition of III kind 5.7. Stage 1

RI PT

After solving Eqns (77),(85), (87) and (88), we obtained

F oq C + I + Pf2 F oq I − F oT U100 C P¯ + Pf2 I − U100 C P¯ 2 = (Pf2 θb + Pm )H Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) Bi (θw − θr ) H1 − H2 + H1 − H2(114) + h h k k

SC

5.8. Stage 2

M AN U

After solving Eqns (78),(91) and (93), we obtained   1 ¯ o1 ∗)F − D1 U 00 P¯ = Bi (θw − θr ) H1 − Bi (θw − θr ) H2 1+ D1 I + D1 U200 P¯ ψ(F 2 Ste2 h h Bi (θw − θr ) Bi (θw − θr ) H1 − H2(115) + k k After solving Eqns (79),(92) and (94), we obtained

5.9. Stage 3

(h2

1 dF + k2)

(116)

TE D

¯ o1 ∗)F − a12 D2 U 00 P¯ = D2 I + a12 D2 U300 P¯ ψ(F 3

After solving Eqns (81),(97) and (100), we obtained

EP

¯ o∗ )F − E1 U 00 P¯ = Bi (θw − θr ) H1 − Bi (θw − θr ) H2 E1 I + E1 U400 P¯ ψ(F 2 4 h h Bi (θw − θr ) Bi (θw − θr ) + H1 − H2 (117) k k

AC C

After solving Eqns (82),(98) and (101), we obtained   1 1 ¯ o∗ )F − a23 E2 U 00 P¯ = 1− E2 I + a23 E2 U500 P¯ ψ(F H 0 (118) 2 5 2 Ste3 (h + k 2 ) 2 After solving Eqns (83),(99) and (102), we obtained ¯ o∗ )F − a13 E3 U 00 P¯ = E3 I + a13 E3 U600 P¯ ψ(F 2 6

(h2

1 H0 + k2) 3

(119)

The problem in all the three Stages with different boundary condition converted into generalized system of Sylvester equations (90, 95, 96, 103 − 117) 33

ACCEPTED MANUSCRIPT

whose solution is obtained by applying Bartels-Stewart algorithm (1972). And applying these results in interface condition, the position of moving

RI PT

layer thickness in Stage 2 and Stage 3 with generalized boundary condition determined. For numerical computation, Matlab Software is used. 6. Model verification for one-dimensional case

SC

The exact solutions of the equations (47)-(56) are available for semi-

infinite one-dimension case. We have described the structure of one-dimensional problem in cryosurgery as a semi-infinite domain, which is initially in liquid

M AN U

phase at temperature T0 . The surface at x = 0 cooled by a cryoprobe at temperature Tc . The freezing starts and the domain is divided into three regions. In Stage 1 only unfrozen region is formed,in Stage 2 mushy and unfrozen regions are formed while in Stage 3 all the three regions : frozen, mushy and unfrozen are formed. We used Boundary fixation method to find

TE D

the exact solution in Stage 3. Exact solution for this problem is as follows:θf = −(C + 4)i2 erf c(e1 z) + Ci2 erf c(−e1 z)

EP

where

4i2 erf c(e1 ) i2 erf c(−e1 ) − i2 erf c(−e1 )

AC C

C=

q q ¯3 ¯3 1+Ste erf (e1 z 1+aSte ) − erf (e ) 1 a23 23 q q θm = ¯3 ¯3 ) − erf (e1 1+aSte ) erf (e2 1+aSte 23 23

and

√ erf c(e1 z a13 ) θu = 2 − . √ erf c(e2 a13 )

34

ACCEPTED MANUSCRIPT

where L a2 a1 X S1 S2 , a23 = , a13 = , z = , e1 = √ , e2 = √ , cm (Ts − Tl ) a3 a3 S1 2 F0 2 F0 Z∞ in erf c(x) = in−1 erf c(t)dt; n = 0, 1, 2... x

RI PT

¯3= Ste

Modified Legendre wavelet Galerkin Method is presented in Section 5 of each

SC

Stage. The temperature was calculated by this method in all the three region. And with the help of temperature, we find the interface position.

M AN U

Figures 2 and 3 show the temperature distribution and interface position calculated by Modified Legendre wavelet Galerkin Method, and their comparison with exact solutions. These figures 2 and 3 clearly demonstrate that numerical Modified Legendre wavelet Galerkin Method agrees completely with exact solution.

TE D

7. Result and Discussion

In this paper, a two-dimensional mathematical model for the freezing of tumor tissue in lung is developed. The procedure of freezing is accom-

EP

plished in three stages by imposing on it the boundary condition of I kind or II kind or III kind. For easy understanding of the model, we established the non-dimensional parameters described in section 3. Consequently, the

AC C

dimensionless form of the model is described in Eqs. (26-42). We applied Modified Legendre wavelet Galerkin method (Kumar et al., 2018) in each phase of boundary condition I, II and III kind for the solution of a dimensionless model. In stage 2 and 3, the model is a moving boundary problem of partial differential equations. Eq. (39) shows the interface condition in stage 2 and Eq.(41, 42) represent the interface condition in stage 3.

We used interface condition of stage 2 and 3 to determine the moving layer 35

ACCEPTED MANUSCRIPT

thickness and obtain the values of λ2 (F o2 ) in stage 2 and λ1 (F o3 ), λ2 (F o3 ) in stage 3. In both stages, the effect of Stefan number on Moving layer thickness

RI PT

is observed. Moving layer thickness increases as the Stefan number decreases. Also, the effect of Kirchoff number on moving layer thickness has been seen in Stage 3 of boundary condition II kind and the effect of Biot number on moving layer thickness has been seen in stage3 of boundary condition III

SC

kind. Moving layer thickness increases as the Kirchoff number increases and moving layer thickness decreases as the Biot number increases.

Moving layer thickness and temperature distribution are two important

M AN U

factors during the cryosurgical treatment of lung-tumor tissue for the prophecy of extreme damage to diseased tissue and the least damage to healthy lung tissue (Kumar et al., 2018). Consequently, we have analyzed the temperature distribution in all stages with generalized boundary condition and moving layer thickness in the mushy and frozen region of boundary condition I, II

TE D

and III kind. We have extracted the parameter of lung and tumor tissue for numerical computation from table 1 (Kumar et al., 2018 Katiyar et al., 2007; Bischof et al., 1992) given below:

EP

Stage I:

Fig .4 exhibits the graph between temperature distribution and time for

AC C

the generalized boundary condition in Stage I. In this stage, we see the temperature distribution at 200 seconds. In this figure(Fig.4) we see the difference in the temperature distribution by keeping it at a constant temperature(I kind), a constant heat flux(II kind) and a constant heat transfer coefficient(III kind). We see how much temperature distribution vary by imposing boundary condition of I, II and III kind. In this stage, temperature decreases rapidly as the time increases. Here the lung tumor is cooled from initial temperature T0 (370 C) to liquidus temperature Tl (−10 C), and also the 36

RI PT

ACCEPTED MANUSCRIPT

SC

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

AC C

EP

TE D

M AN U

Parameter unit of measurement Density of unfrozen lung tissue kg/m3 Density of frozen lung tissue kg/m3 Density of unfrozen tumor tissue kg/m3 Density of frozen tumor tissue kg/m3 Density of blood kg/m3 Thermal conductivity of unfrozen lung tissue W/m0 C Thermal conductivity of frozen lung tissue W/m0 C Thermal conductivity of unfrozen tumor tissue W/m0 C Thermal conductivity of frozen tumor tissue W/m0 C Specific heat of unfrozen lung tissue J/kg 0 C Specific heat of frozen lung tissue J/kg 0 C Specific heat of unfrozen tumor tissue J/kg 0 C Specific heat of frozen tumor tissue J/kg 0 C Blood perfusion in lung tissue ml/s/ml Blood perfusion in tumor tissue ml/s/ml Metabolic heat generation in lung W/m3 Metabolic heat generation in tumor W/m3 Latent heat kJ/kg 0 Arterial blood temperature C

37

value 161 149 998 921 1005 0.11 0.38 0.552 2.25 4174 1221 4200 1230 0.0005 0.002 42, 000 672 333 37

ACCEPTED MANUSCRIPT

temperature (Tu ) decreases as space coordinate x, y increases. We obtained the non-dimensional temperature Θu by applying Modified Legendre wavelet

RI PT

Galerkin Method and then obtained Tu . Stage II:

In this stage, Fig. 5, 6 exhibits the graph between temperature distri-

SC

bution and time for the generalized boundary condition. In this stage, we see the temperature distribution at 500 seconds. In these figures (Fig. 5,6) we see the difference in the temperature distribution by keeping it at a con-

M AN U

stant temperature(I kind), a constant heat flux(II kind) and a constant heat transfer coefficient(III kind). We see how much temperature distribution vary by imposing boundary condition of I, II and III kind. As observed the temperature of unfrozen region diminishes slowly as compared to the mushy region. In this stage, the lung tumor is cooled from the liquidus temperature

TE D

Tl (−10 C) to solidus temperature Ts (−80 C). We obtained the dimensionless temperature Θm and Θu by applying Modified wavelet Galerkin Method and then obtained Tm and Tu . We see that the temperature distribution in unfrozen(Tu ) and mushy region(Tm ) decreases as space coordinate x, y

EP

increases. We also see that there is slight variation not too much in temperature distribution by imposing generalized boundary condition (I, II and

AC C

III kind). Fig. 7 shows the comparison of dimensionless moving layer thickness λ2 (F o2 ) in the mushy-liquid region of boundary condition I, II and III kind and moving layer thickness λ2 (F o2 ) increases as the dimensionless time F o2 increases. The effect of Stefan number on moving layer thickness is also observed in Fig. 13.

38

ACCEPTED MANUSCRIPT

Stage III: Here, the solidus moving front λ1 and the liquidus moving front λ2 are

RI PT

calculated. Fig.8 exhibits the dimensionless moving layer thickness λ1 (F o3 )

in solid-mush region and Fig. 9 exhibits the moving layer thickness λ2 (F o3 )

in the mush-liquid region by imposing on it boundary condition of I, II and III kind. We see that λ1 (F o3 ) and λ2 (F o3 ) increases as the dimensionless

SC

time F o3 increases. Also, the effect of Stefan number is observed on moving

layer thickness λ1 and λ2 in Fig. 14 and 15 respectively. We obtained the

M AN U

non-dimensional temperature Θf , Θm and Θu by applying Modified wavelet Galerkin Method and then obtained Tf , Tm and Tu . The temperature distribution in frozen region(Tf ), mushy region(Tm ) and unfrozen region(Tu ) decreases as space coordinate x, y increases. Fig. 10, 11 and 12 shows the graph between temperature distribution and time by keeping it at a constant temperature(I kind), a constant heat flux(II kind) and a constant heat trans-

TE D

fer coefficient(III kind). In this stage, we see the temperature distribution at 900 seconds. In these figures, we see the difference in the temperature distribution of boundary condition I, II and III kind. We see how much tem-

EP

perature distribution vary with different boundary condition. As observed the temperature of the frozen region decreases rapidly as compared to the

AC C

mushy region and the unfrozen region. Fig. 8, 9 exhibits the dimensionless moving layer thickness λ1 (F o3 ) in the solid-mush region and λ2 (F o3 ) in mush-liquid region respectively. As the time F o3 increases, moving layer thickness λ1 (F o3 ) and λ2 (F o3 ) increases. The effect of Stefan number on moving layer thickness is seen in Fig. 14 and 15. Also, the effect of the Kirchoff number and Biot number is seen in Fig. 16 and 17. Here the tissue is cooled from freezing temperature to the lethal temperature. In particular, the lethal temperature for tissue destruc39

ACCEPTED MANUSCRIPT

tion usually begins around −400 C (Zhang, 2009) for diseased tissue like a

RI PT

tumor.

8. Conclusion

In this study, a two-dimensional mathematical bio-heat transfer model

of lung tumor tissue during the freezing process in cryosurgery has been de-

SC

veloped and then used the Modified Legendre wavelets Galerkin method to obtain the results. Firstly, our problem is transformed into non-dimensional

M AN U

form and then applying finite difference method in our problem to convert it into initial boundary value problem of ordinary matrix differential equation in stage 1 and moving boundary value problem in stage 2 and 3 by imposing on it at a constant temperature(I kind), a constant heat flux(II kind) and a constant heat transfer coefficient(III kind) in each stage. Af-

TE D

ter this, we obtained the system of Sylvester equations by using Legendre wavelet Galerkin method which are solved by generalized inverse technique (Bartels-Stewart algorithm). In each stage, we have applied this technique to carry out the temperature distribution. From our calculations, we see that

EP

the effect of relaxation time is negligible in the mushy and frozen region and present only in unfrozen region. We observe from the figures of moving layer

AC C

thickness and temperature distribution that there is slight variation not too much with the effect of different boundary condition. In Stage1, the tissue is cooled up to the liquidus temperature Tl (−10 C) where unfrozen region is

formed and then in Stage 2, the tissue is cooled up to the freezing temperature Ts (−80 C) where the mushy region is formed and in Stage3, the tissue is cooled up to the lethal temperature where the frozen region is formed. The effect of Stefan number on moving layer thickness is seen in stage 2 and 3. Also, the effect of Kirchoff number on moving layer thickness has been 40

ACCEPTED MANUSCRIPT

seen in Stage 3 of boundary condition II kind and the effect of Biot number on moving layer thickness has been seen in stage3 of boundary condition III

RI PT

kind. Moving layer thickness increases as the Kirchoff number increases and moving layer thickness decreases as the Biot number increases. Although a

two-dimensional mathematical bio-heat transfer model is considered in our investigation, it can further be extended for the three-dimensional mathe-

SC

matical bio-heat transfer model to attain more physically realistic results.

Lastly, this model is beneficial for the experimental analysis of the freezing

Acknowledgments

M AN U

process of cryosurgical treatment.

The research of the first author is supported by University Grant Commission New Delhi, India, grant no. 19/06/2016(i)EU-V. The authors express their sincere thanks to DST-Centre for Interdisciplinary Mathematical Sci-

TE D

ences, Banaras Hindu University Varanasi, India, for alloting the required facilities. The authors are grateful to the Reviewer for their valuable comments.

EP

References

Ahmadikia, H., Moradi, A., 2012. Non-Fourier phase change heat

AC C

transfer in biological tissues during solidification. Heat. Mass Transf. 48 (9), 1559–1568. Allington, H.V., 1950. Liquid nitrogen in the treatment of skin diseases. Calif. Med. 72 (3), 153–155. Antaki, P.J., 2005. New interpretation of Non-Fourier heat conduction in processed meat. ASME J. Heat Transf. 127, 189–193. Bartels R. and Stewart G., Solution of the Matrix Equation AX + XB = C: Algorithm 432, Comm. ACM, 15 ,pp.820–826(1972). 41

ACCEPTED MANUSCRIPT

Bischof, J.C., Bastacky, J., Rubinsky, B., 1992. An analytical study of cryosurgery in the lung. J. Biomech. Eng. 114 (4), 467–472.

RI PT

Budman, H., Shitzer, A., Dayan, J., 1995. Analysis of the inverse problem of freezing and thawing of a binary solution durin cryosurgical processes. Journal of Biomechanical Engineering 177,193–202.

Cattaneo, C., 1958. A form of heat conduction equation which elimi-

SC

nates the paradox of instantaneous propagation. Comptes Rendus 247, 431–433.

Chua, K.J., Chou, S.K., Ho, J.C., 2007. An analytical study on the

40 (1), 100–116.

M AN U

thermal effects of cryosurgery on selective cell destruction. J. Biomech.

Deng Z.S. , Liu J., Numerical simulation of selective freezing of target biological tissues following injection of solutions with specific thermal properties, Cryobiology, 50(2),183-192(2005).

TE D

Frayn, K.N., 1996. Metabolic Regulation:

A Human Perspective.

Portland Press, London, pp. 84–88. Gupta SC, The Classical Stefan Problems: Basic Concepts,Modelling and Analysis, Elsevier, 2003.

EP

Islami F, Goding Sauer A, Miller KD, Siegel RL, Fedewa SA, Jacobs EJ, McCullough ML, Patel AV, Ma J,Soerjomataram I, Flanders WD,

AC C

Brawley OW, Gapstur SM, Jemal A (January 2018). ”Proportion and number of cancer cases and deaths attributable to potentially modifiable risk factors in the United States”. Ca. 68 (1):

31–54.

doi:10.3322/caac.21440. PMID 29160902. Kotova T. G., Kochenov V. I., Tsybusov S. N. , Madai D. Y., Gurin A. V. , Calculation of Effective Freezing Time in Lung Cancer Cryosurgery based on Godunov Simulation, CTM, 8(1), 2016.

42

ACCEPTED MANUSCRIPT

Kumar, A., Kumar, S., Katiyar, V.K., Telles, S., 2017. Phase change heat transfer during cryosurgery of lung cancer using hyperbolic heat

RI PT

conduction model. Comput. Biol.Med. 84, 20–29. Math. Mech. 3 (3) (2007) 1–17.

Kumar, Ajay, Kumar, Sushil, Katiyar, V.K., Telles, Shirley, 2017.

Dual phase lag bio-heat transfer during cryosurgery of lung cancer:

SC

comparison of three heat transfer models. J. Therm. Biol. 69.

Kumar D. ,Singh S. , Rai K.N., Analysis of classical fourier, SPL and DPL heat transfer model in biological tissues in presence of metabolic

M AN U

and external heat source. Heat Mass Transf. 52 (6), 1089–1107(2016). Kumar D , kumar P., Rai K., N.,, A study on DPL model of heat transfer in bi-layer tissues during MFH treatment, Computers in Biology and Medicine 75 (2016) 160–172.

Kumar, Mukesh, Upadhyay, S., Rai, K.N.,2018. A study of cryosurgery

TE D

of lung cancer using Modified Legendre wavelet Galerkin method. J. of Thermal Biology. 78,356-366.

Kumar, P., kumar, D., Rai, K.N., 2015. A numerical study on dual phase lag model of bio heat transfer during hyperthermia treatment. J.

EP

Therm. Biol. 49, 98–105. Kumar, S., Katiyar V.K. , Numerical study on phase change heat

AC C

transfer during combined hyperthermia and cryosurgical treatment of lung cancer, Int. J. Appl., 2016. Lee, C.Y., Bastacky, J., 1995. Comparative mathematical analyses of freezing in lung and solid tissue. Cryobiology 32, 299–305. Liu, K.C., Chen, H., 2009. Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment. Int. J. Heat. Mass Transf. 52, 1185–1192.

43

ACCEPTED MANUSCRIPT

Majchrzak, E., 2010. Numerical solution of dual phase lag model of bioheat transfer using the general boundary element method. Comput.

RI PT

Model. Eng. Sci. 69 (1), 43–60. Paola Pasquali, Cryosurgery: A Practical Manual (2015).

Rabin, Y., Shitzer, A., 1995. Exact solution to the one-dimensional inverse-Stefan problem in non ideal biological tissue. Journal of Heat

SC

Transfer 117, 425–431.

Razzaghi M. and Yousefi S., Legendre wavelets operational matrix

495–502, (2001).

M AN U

of integration, International Journal of System Science,vol. 32, pp.

Tarwidi, D., 2015. Godunov method for computerized lung cancer cryosurgery planning with efficient freezing time. In:

Proceedings

of the 3rd International Conference on IEEE the Information and Communication Technology (ICoICT), 2015, pp. 494–499.

TE D

Tolar J, Neglia JP (June 2003). ”Transplacental and other routes of cancer transmission between individuals”. Journal of Pediatric Hematology/Oncology. 25 (6): 430–436. doi:10.1097/00043426-200306000-00002. PMID 12794519.

EP

Tzou, D.Y., 1995. A unified field approach for heat conduction from macro to microscales. J. Heat Trans. (117), 8–16.

AC C

Yadav S , Kumar D , Rai K.N., Finite Element Legendre Wavelet Galerkin Approach to inward solidification in simple body under most generalised boundary condition, Z. Naturforsch. 69a, 501-510(2014). Yadav S. , Upadhyay S. , Rai K.N. , A Mathematical Model for solidification of Binary Eutectic System including relaxation time, Computational Thermal Sciences,8,583-589(2016). Zhang, Y., 2009. Generalized dual-phase lag bio-heat equations based

44

ACCEPTED MANUSCRIPT

on non equilibrium heat transfer in living biological tissues. Int. J. Heat. Mass Transf. 52, 4829–4834.

RI PT

Zhou, J., Zhang, Y., Chen, J.K., 2009a. An axisymmetric dual-phase-lag bioheat model for laser heating of living tissues. Int. J. Therm. Sci. 48,

EP

TE D

M AN U

SC

1477–1485.

AC C

Figure 2: Temperature distribution T(x,t) for one-dimensional at time t=900 seconds calculated using physical properties(Table 1); T0 = 370 C, Tc = −1960 C; length = 0.5m.

45

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

Figure 3: Interface position for one-dimensional case considering the time calculated using physical properties(Table 1); T0 = 370 C, Tc = −1960 C; length = 0.5m.

46

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

47

Figure 4: Temperature distribution of unfrozen region in Stage1 at 200 sec (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

48

Figure 5: Temperature distribution of unfrozen region in Stage2 at 500 sec (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

49

Figure 6: Temperature distribution of mushy region in Stage2 at 500 sec (a) I kind (b) II kind (c) III kind.

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

Figure 7: Moving layer thickness λ2 in stage2.

50

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

Figure 8: Moving layer thickness λ1 in stage3.

51

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

Figure 9: Moving layer thickness λ2 in stage3.

52

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

53

Figure 10: Temperature distribution of unfrozen region in Stage3 at 900 sec (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

54

Figure 11: Temperature distribution of mushy region in Stage3 at 900 sec (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

55

Figure 12: Temperature distribution of frozen region in Stage3 at 900 sec (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

56

Figure 13: Moving layer thickness λ2 in stage2 with different Stefan number (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

57

Figure 14: Moving layer thickness λ1 in stage3 with different Stefan number (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

58

Figure 15: Moving layer thickness λ2 in stage3 with different Stefan number (a) I kind (b) II kind (c) III kind.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 16: Moving layer thickness in stage3 with different Kirchoff number (a).λ1 (F o3 ) (b). λ2 (F o3 )

59

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 17: Moving layer thickness in stage3 with different Biot number (a).λ1 (F o3 ) (b). λ2 (F o3 )

60

ACCEPTED MANUSCRIPT

Highlight

RI PT

SC

•

M AN U

•

TE D

•

EP

•

In this study, we have developed two-dimensional mathematical model describing bio-heat transfer during cryosurgery of lung cancer. We see the effect of imposition of different type of boundary condition on temperature distribution and moving layer thickness. We find the temperature distribution in all the three stages and moving layer thickness in stage 2 & 3. Our model will be very useful for the experimental analysis of freezing process of cryosurgical treatment in lung cancer. We will use this novel method in various type of moving boundary problem.

AC C

•

ACCEPTED MANUSCRIPT

Author’s Affiliation 1st and Corresponding author: Name: Mukesh Kumar

Email-id: [email protected] 2nd author:

SC

Name: Dr. Subrahamanyam Upadhyay

RI PT

Address: DST-CIMS, Institute of Science, Banaras Hindu University, Varanasi

Address: Department of Mathematics, Eternal University, Solan (H.P.)

3rd author: Name: Prof. K. N. Rai

M AN U

Email-id: [email protected]

Address: Department of Mathematical Sciences, IIT(BHU), Varanasi

AC C

EP

TE D

Email-id: [email protected]