JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

44, 2, pp. 381-392, Warsaw 2006

SENSITIVITY ANALYSIS OF BIOLOGICAL TISSUE

FREEZING PROCESS WITH RESPECT TO THE RADIUS

OF SPHERICAL CRYOPROBE

Ewa Majchrzak

Grażyna Kałuża

Department for Strength of Materials and Computational Mechanics, Silesian University of

Technology, Gliwice

e-mail: ewa.majchrzak@polsl.pl; grazyna.kaluza@polsl.pl

Application of the shape sensitivity analysis to the case of problems of the
biological tissue freezingprocess is discussed.The freezingprocessdescribed
by a strongly non-linear bioheat transfer equation in which an additional
term controlling the evolution of latent heat appears. Using the approach
called ’the one domain method’, one finally obtains a partial differential
equation containing the substitute thermal capacity of tissue.Theboundary
and initial conditions determine the thermal interaction between the tissue
and cryoprobe tip.
In the paper, we consider a spherical internal cryoprobe. In order to es-
timate the influence of cryoprobe geometry on the course of the process,
the shape sensitivity analysis is applied. In particular, a direct approach is
used (explicit differentiation method). The results of numerical modeling
(the boundary element method is applied) allow one to formulate essential
practical conclusions concerning the course of cryosurgery treatments.

Key words: freezing process, sensitivity analysis, boundary elementmethod

1. Governing equations

From the mathematical point of view, the freezing process belongs to a
group ofmoving boundary problems because the shape and dimensions of the
frozen region are time-dependent. In the case of internal spherical cryoprobe
application (Fig.1), the following equationwritten in the spherical co-ordinate
system should be taken into account

R1 <r<R2 : C(T)
∂T(r,t)

∂t
=
1

r2
∂

∂r

[

λ(T)r2
∂T(r,t)

∂r

]

(1.1)



382 E.Majchrzak, G.Kałuża

where C(T) [J/(m3K)] is the substitute thermal capacity per unit volume,
λ(T) [W/(mK)] is the thermal conductivity, T , r, tdenote temperature, spatial
co-ordinate and time. The courses of C(T) and also λ(T) are presented in
Fig.2 (Comini and Giudice, 1976; Majchrzak and Dziewoński, 2000).

Fig. 1. Spherical cryoprobe

Fig. 2. (a) Substitute thermal capacity; (b) thermal conductivity

Equation (1.1) is supplemented by the following boundary conditions














r=R1 : T(r,t) =Tc

r=R2 :
∂T(r,t)

∂r
=0

(1.2)

and initial condition
t=0 : T(r,t)=T0 (1.3)



Sensitivity analysis of biological tissue... 383

where R1 is the cryoprobe radius, R2 is the conventionally assumed exter-
nal radius of the domain considered, Tc is the temperature of the cryoprobe
surface, T0 is the initial temperature of the tissue.
In order to use the boundary elementmethod, linearization of the problem

discussedmustbe introduced.Here, the artificial heat sourcemethod is applied
(Majchrzak and Mochnacki, 1996). This method requires transformation of
the governing equations and boundary-initial conditions by the introduction
of Kirchhoff’s function, namely

V (T)=

T
∫

Tr

λ(µ) dµ (1.4)

where Tr is an arbitrary assumed reference level.
Using this function, we can transform governing equations (1.1), (1.2),

(1.3) to the form

Φ[T(V )]
∂V (r,t)

∂t
=
1

r2
∂

∂r

[

r2
∂V (r,t)

∂r

]

(1.5)

where (cf. Fig.3)

Φ[T(V )] =
C[T(V )]

λ[T(V )]
(1.6)

Boundary and initial conditions (1.2), (1.3) should be also transformed and
then























r=R1 : V (r,t)=Vc

r=R2 :
∂V (r,t)

∂r
=0

t=0 : V (r,t)=V0

(1.7)

where Vc =V (Tc), V0 =V (T0).

2. Sensitivity analysis

In order to estimate the influence of the cryoprobe radius on the course of
the freezing process, a sensitivity model is constructed. Using the concept of
material derivative (Dems, 1987; Kleiber, 1997) one can write

DV

Db
=
∂V

∂b
+
∂V

∂r
v (2.1)



384 E.Majchrzak, G.Kałuża

Fig. 3. Course of the function Φ[T(V )]

where v= v(r,b) is the velocity associated with the design parameter b=R1.
Equation (1.5) can be written in the form

Φ[T(V )]
∂V

∂t
=
∂2V

∂r2
+
2

r

∂V

∂r
(2.2)

Using the direct approach of sensitivity analysis (Dems, 1987; Kałuża, 2005;
Kleiber, 1997), equation (2.2) is differentiated with respect to the shape pa-
rameter b, namely

DΦ[T(V )]

Db

∂V

∂t
+Φ[T(V )]

D

Db

(∂V

∂t

)

=
D

Db

(∂2V

∂r2

)

+2
D

Db

(1

r

∂V

∂r

)

(2.3)

Because (c.f. formula (2.1))

D

Db

(∂V

∂r

)

=
∂

∂b

(∂V

∂r

)

+
∂

∂r

(∂V

∂r

)

v=
∂

∂r

(∂V

∂b

)

+
∂2V

∂r2
v (2.4)

and

∂

∂r

(DV

Db

)

=
∂

∂r

(∂V

∂b
+
∂V

∂r
v
)

=
∂

∂r

(∂V

∂b

)

+
∂2V

∂r2
v+
∂V

∂r

∂v

∂r
(2.5)

therefore
D

Db

(∂V

∂r

)

=
∂

∂r

(DV

Db

)

−
∂V

∂r

∂v

∂r
(2.6)

Next, the material derivative of the component ∂2V/∂r2 is calculated. From
equation (2.1), one obtains

D

Db

(∂2V

∂r2

)

=
D

Db

[ ∂

∂r

(∂V

∂r

)]

=
∂

∂r

[ D

Db

(∂V

∂r

)]

−
∂2V

∂r2
∂v

∂r
(2.7)



Sensitivity analysis of biological tissue... 385

this means
D

Db

(∂2V

∂r2

)

=
∂2

∂r2

(DV

Db

)

−2
∂2V

∂r2
∂v

∂r
−
∂V

∂r

∂2v

∂r2
(2.8)

In a similar way, one finds

D

Db

(1

r

∂V

∂r

)

=
1

r

[ ∂

∂b

(∂V

∂r

)

+
∂

∂r

(∂V

∂r

)

v
]

−
1

r2
∂V

∂r
v (2.9)

Taking into account formula (2.1) and dependence (2.6), we have

D

Db

(1

r

∂V

∂r

)

=
1

r

D

Db

(∂V

∂r

)

−
1

r2
∂V

∂r
v=
1

r

∂

∂r

(DV

Db

)

−
1

r

∂V

∂r

∂v

∂r
−
v

r2
∂V

∂r
(2.10)

Thematerial derivative of the component ∂V/∂t is also calculated

D

Db

(∂V

∂t

)

=
∂

∂b

(∂V

∂t

)

+
∂

∂r

(∂V

∂t

)

v=
∂

∂t

(∂V

∂b
+
∂V

∂r
v
)

=
∂

∂t

(DV

Db

)

(2.11)

From equation (2.2), it results that

∂2V

∂r2
=Φ[T(V )]

∂V

∂t
−
2

r

∂V

∂r
(2.12)

Using formulas (2.8), (2.10) and (2.11), equation (2.3) takes form

Φ[T(V )]
∂U(r,t)

∂t
=
∂2U(r,t)

∂r2
+
2

r

∂U(r,t)

∂r
+

(2.13)

−
(

2Φ[T(V )]
∂v

∂r
+
DΦ[T(V )]

Db

)∂V (r,t)

∂t
+
(2

r

∂v

∂r
−
∂2v

∂r2
−
2v

r2

)∂V (r,t)

∂r

where

U(r,t) =
DV (r,t)

Db
(2.14)

is the sensitivity function, while

DΦ[T(V )]

Db
=

1

λ2[T(V )]

(dC

dT
−
dλ

dT
Φ[T(V )]

)

U(r,t) (2.15)

Boundary and initial conditions (1.7) are also differentiated with respect
to b







































r=R1 :
DV

Db
=
DVc
Db
=0

r=R2 :
Dq

Db
=−
D

Db

(∂V

∂r

)

=−
∂

∂r

(DV

Db

)

+
∂V

∂r

∂v

∂r
=0

t=0 :
DV

Db
=
DV0
Db
=0

(2.16)



386 E.Majchrzak, G.Kałuża

this means














r=R1, : U(r,t)= 0

r=R2 : W(r,t)= 0

t=0 : U(r,t)= 0

(2.17)

where

W(r,t)=−
∂U(r,t)

∂r
(2.18)

In equation (2.13), the velocity field v(r,b) associated with the design para-
meter b=R1 is defined as follows

v= v(r,b) =
R2−r

R2− b
(2.19)

It should be pointed out that the additional problem (c.f. equations (2.13),
(2.17)) connected with the sensitivity function U(r,t) is coupled with the
basic one (c.f. equations (1.5), (1.7)) by the functions V (r,t) and Φ[T(V )].

3. Method of solution

The basic problem and the additional one have been solved using a combi-
nedvariant ofBEM(Brebbia andDominiguez, 1992;Kałuża, 2005;Majchrzak,
2001) supplemented by the artificial heat source procedure (Majchrzak and
Mochnacki, 1996).

Let us consider the following equation

Φ[T(V )]
∂F(r,t)

∂t
=
1

r2
∂

∂r

[

r2
∂F(r,t)

∂r

]

+R(r,t) (3.1)

where for the basic problem F(r,t) =V (r,t), R(r,t) = 0 (cf. equation (1.5)),
while for the additional problem F(r,t) = U(r,t) and (cf. equations (2.13),
(2.15))

R(r,t)=
2Φ[T(V )]

R2− b

∂V (r,t)

∂t
+

1

λ2[T(V )]

(dλ

dT
Φ[T(V )]−

dC

dT

)

U(r,t)
∂V (r,t)

∂t
+

(3.2)

−2
[ 1

r(R2− b)
+
R2−r

r2(R2− b)

]∂V (r,t)

∂r



Sensitivity analysis of biological tissue... 387

Using the artificial heat sourcemethod ((Majchrzak andMochnacki, 1996) the
function Φ[T(V )] is expressed as a sum of two components, which entails a
constant part Φ0 and a certain increment ∆Φ[T(V )], see Fig.3

Φ[T(V )] =Φ0+∆Φ[T(V )] (3.3)

Equation (3.1) can be written in the form

Φ0
∂F(r,t)

∂t
=
1

r2
∂

∂r

[

r2
∂F(r,t)

∂r

]

+S(r,t) (3.4)

where

S(r,t)=R(r,t)−∆Φ[T(V )]
∂F(r,t)

∂t
(3.5)

is the source function (capacity of internal heat sources).
The essential feature of equation (3.4) consists in the fact that leaving out

the last term, we obtain a linear form of energy equation. The calculation of
the source function requires the introduction of a certain iterative procedure
(Majchrzak andMochnacki, 1996).
In order to solve equation (3.4), BEM using discretization in the time

domain is applied (Brebbia and Dominiguez, 1992; Kałuża, 2005; Majchrzak,
2001). At first, the time derivative appearing in equation (3.4) is substituted
by a differential quotient

t∈ [tf−1, tf] :
∂F(r,t)

∂t
=
F(r,tf)−F(r,tf−1)

∆t
(3.6)

and then equation (3.4) takes form

∂

∂r

[

r2
∂F(r,tf)

∂r

]

−
Φ0
∆t
r2F(r,tf)+

Φ0
∆t
r2F(r,tf−1)+r2S(r,tf)= 0 (3.7)

Next, the weighted residual criterion is applied

R2
∫

R1

{ ∂

∂r

[

r2
∂F(r,tf)

∂r

]

−
Φ0
∆t
r2F(r,tf)+

(3.8)

+
Φ0
∆t
r2F(r,tf−1)+r2S(r,tf)

}

F∗(ξ,r) dr=0

where ξ is the observation point and F∗(ξ,r) is the fundamental solution

F∗(ξ,r)=
1

2rξ

√

∆t

Φ0

[

exp
(

−|r−ξ|

√

Φ0
∆t

)

− exp
(

−|r+ ξ|

√

Φ0
∆t

)]

(3.9)



388 E.Majchrzak, G.Kałuża

After mathematical transformations of criterion (3.8), one obtains

F(ξ,tf)+
[

r2F∗(ξ,r)J(r,tf)
]r=R2

r=R1
=
[

r2J∗(ξ,r)F(r,tf)
]r=R2

r=R1
+

(3.10)

+
Φ0
∆t

R2
∫

R1

r2F(r,tf−1)F∗(ξ,r) dr+

R2
∫

R1

r2S(r,tf)F∗(ξ,r) dr

where

J∗(ξ,r)=−
∂F∗(ξ,r)

∂r
=
1

r
F∗(ξ,r)+

(3.11)

+
1

2rξ

[

sgn(r− ξ)exp
(

−|r− ξ|

√

Φ0
∆t

)

− sgn(r+ ξ)exp
(

−|r+ ξ|

√

Φ0
∆t

)]

and

J(r,tf)=−
∂F(r,tf)

∂r
(3.12)

Equation (3.10) can be written in the form

F(ξ,tf)+R22F
∗(ξ,R2)J(R2, t

f)−R21F
∗(ξ,R1)J(R1, t

f)=
(3.13)

=R22J
∗(ξ,R2)F(R2, t

f)−R21J
∗(ξ,R1)F(R1, t

f)+P(ξ)+Z(ξ)

where

P(ξ)=
Φ0
∆t

R2
∫

R1

r2F(r,tf−1)F∗(ξ,r) dr

(3.14)

Z(ξ)=

R2
∫

R1

r2S(r,tf)F∗(ξ,r) dr

For ξ→R+1 and ξ→R
−

2 , one obtains the following system of equations

[

−R21F
∗(R1,R1) R

2
2F
∗(R1,R2)

−R21F
∗(R2,R1) R

2
2F
∗(R2,R2)

][

J(R1, t
f)

J(R2, t
f)

]

=

(3.15)

=

[

−R21J
∗(R+1 ,R1)−1 R

2
2J
∗(R+1 ,R2)

−R21J
∗(R−2 ,R1) R

2
2J
∗(R−2 ,R2)−1

][

F(R1, t
f)

F(R2, t
f)

]

+

[

P(R1)
P(R2)

]

+

[

Z(R1)
Z(R2)

]



Sensitivity analysis of biological tissue... 389

fromwhich theunknownboundaryvalues J(R1, t
f),F(R2, t

f) aredetermined.
Values of the function F at the internal points ξ ∈ (R1,R2) are calculated
using the formula

F(ξ,tf)=R22J
∗(ξ,R2)F(R2, t

f)−R21J
∗(ξ,R1)F(R1, t

f)+
(3.16)

−R22F
∗(ξ,R2)J(R2, t

f)+R21F
∗(ξ,R1)J(R1, t

f)+P(ξ)+Z(ξ)

4. Results of computations

A cryoprobe of diameter R1 = 0.005m has been considered. The surface
temperature is assumed as −90◦C. The external radius of the tissue domain:
R2 = 0.025m. The initial temperature of tissue T0 = 37

◦C. The domain of
tissue has been divided into 600 linear internal cells, time step ∆t=0.5s and
constant Φ0 =4 ·10

7 (see equation (3.4) and Fig.3).
In Figure 4a, the temperature distribution in the tissue for times instants

5, 10, 15, ..., 60s is shown. Figure 4b illustrates cooling curves at the points
r=0.006, 0.007, 0.008, 0.009 and 0.01m. In Figure 5a, the distribution of the
sensitivity function U =DV/Db in the domain considered for instants 5, 10,
15, ..., 60s is presented. Figure 5b shows courses of the function U at selected
points in the tissue domain.

Fig. 4. (a) Temperature distribution; (b) cooling curves

The function V (r,t,R1) is expanded intoTaylor series taking into account
the first two components

V (r,t,R1±∆R1)=V (r,t,R1)±U(r,t,R1)∆R1 (4.1)



390 E.Majchrzak, G.Kałuża

Fig. 5. (a) Distribution of the function U =DV/Db; (b) courses of the function U
at selected points

Thus, using the solutions V (r,t,R1) and U(r,t,R1), it is possible to obtain
solutions corresponding to cryoprobe radii R1−∆R1 and R1+∆R1. Figures 6
present results obtained this way on the assumption that ∆R1 =0.1R1.

Fig. 6. (a) Temperature distribution for R1−∆R1; (b) temperature distribution for
R1+∆R1

Summing up, the solutions to the basic and additional problem allow
one (using the Taylor formula) to rebuilt the results obtained on the infi-
nite number of other solutions corresponding to the changed cryoprobe radii.
In this way, it is possible to analyze mutual connections between the radius
of the cryoprobe and the course of freezing process proceeding in the tissue
domain.



Sensitivity analysis of biological tissue... 391

This paper is a part of project No. 3T11F01826 sponsored by the State Com-

mittee for Scientific Research (KBN).

References

1. Brebbia C.A., Dominiguez J., 1992, Boundary Elements, an Introductory
Course, ComputationalMechanics Publications,McGraw-Hill BookCompany,
London

2. Comini G., Giudice L.D., 1976, Thermal aspects of cryosurgery, Journal of
Heat and Mass Transfer, 12, 543-549

3. Dems K., 1987, Sensitivity analysis in thermal problems – II: structure shape
variation, Journal of Thermal Stresses, 10, 1-16

4. Kałuża G., 2005, Application of the sensitivity analysis methods in bioheat
transfer, Doctoral thesis Silesian University of Technology, Gliwice, 2005

5. Kleiber M., 1997,Parameter Sensitivity in Nonlinear Mechanics, J.Wiley &
Sons, Chichester

6. Majchrzak E., 2001, Boundary Element Method in Heat Transfer, Publ. of
CzestochowaUniversity of Technology, Częstochowa (in Polish)

7. MajchrzakE., Dziewoński M., 2000,Numerical simulation of freezing pro-
cess using the BEM, Computer Assisted Mechanics and Engineering Sciences,
7, 667-676

8. Majchrzak E., Mochnacki B., 1996, The BEM application for numerical
solution of non-steady and non-linear thermal diffusion problems, Computer
Assisted Mechanics and Engineering Sciences, 3, 4, 327-346

Analiza wrażliwości procesu zamrażania tkanki biologicznej ze względu

na promień kriosondy sferycznej

Streszczenie

Praca dotyczy zastosowania analizy wrażliwości kształtu w zagadnieniachmode-
lowania procesu zamrażania tkanki biologicznej. Proces ten jest opisany silnie nieli-
niowym równaniem przepływu biociepła, w którym pojawia się dodatkowy składnik
związany z wydzielaniem się utajonego ciepła zamrażania. Formalne przekształcenia
wyjściowegorównaniaopisującegoprocesprowadzido równania różniczkowego,wktó-
rym występuje zastępcza pojemność cieplna tkanki (tzw. metoda jednego obszaru).



392 E.Majchrzak, G.Kałuża

Odpowiednio przyjęte warunki brzegowe determinują oddziaływania termiczne mię-
dzy tkanką a końcówką kriosondy.
W pracy rozważa się sferyczną kriosondęwewnętrzną.W celu oszacowaniawpły-

wu geometrii kriosondy na przebieg procesu zamrażaniawykorzystanopodejście bez-
pośrednie analizy wrażliwości kształtu.Wyniki obliczeń numerycznych otrzymane za
pomocą metody elementów brzegowych pozwoliły na sformułowanie wniosków przy-
datnych w praktyce kriochirurgicznej.

Manuscript received December 16, 2005; accepted for print January 23, 2006