Shape sensitivity analysis of temperature distribution in heating tissue

(1)

in domain of biological tissue is presented (2D problem). The temperature distribution results from the action of external heat source emanating the boundary heat flux.

The sensitivity analysis of this process with respect to the geometrical parameters of biological tissue is considered. The explicit differentiation method basing on the material derivative approach is used. On the stage of numerical computations the boundary element method is applied. In the final part of the paper the results obtained are shown.

1. Formulation of the problem

The temperature field in domain of biological tissue (Fig. 1) subjected to the action of boundary heat flux q₀ is described by the following equation [1] and con- ditions:

( ) ( ) ( )

( )

( ) ( )

( ) ( ) ( )

2 2

1

2

0 0

, ,

, : λ λ , 0

, : ,

, : , λ , 0

, : , λ ,

B m

b

T x y T x y

x y k T T x y Q

x y

x y T x y T

x y q x y T x y

x y q x y T x y q

∂ ∂

∈ Ω ∂ + ∂ +  −  + =

∈ Γ =

∈ Γ = − ⋅ ∇ =

n n

(1)

where λ [W/(mK)] is the thermal conductivity, k =cBGB (cB [J/(m³K)] is the volu- metric specific heat of blood, GB [m³blood/s/m³tissue] is the perfusion rate), TB is the blood temperature, Qm[W/m³] is the metabolic heat source, T is the temperature, x, y are the geometrical co-ordinates. The domain Ω is the square of dimen- sions L×L, while the boundary heat flux q₀ changes according to the formula

( )

² ²

0 m 2

L y

q y q

L

= − (2)

where q₀(0) = q_m is the maximum value of this heat flux.

(2)

Fig. 1. Domain considered

We assume that the shape parameter b corresponds to the half of square diago- nal, this means b=L√2/2. The aim of investigations is to estimate the changes of temperature due to the expansion (or contraction) of the domain considered.

2. Explicit differentiation method

Using the concept of material derivative we can write [2, 3]

D

D ^x ^y

T T T T

v v

b b x y

∂ ∂ ∂

= + +

∂ ∂ ∂ (3)

where b is the shape parameter, vx=vx(x,y,b) and vy=vy(x,y,b) are the velocities associated with design parameter b.

If the direct approach of sensitivity method is applied [2, 3] then the equation (1) are differentiated with respect to shape parameter b, and then:

( )

( ) ( )

2 2

1

2

0 0

D D D

, : λ λ 0

D D D

D

, : D 0

D D

, : λ 0

D D

, : D λ =

D D D

b

T T T

x y k

b x b y b

T x y T

b b

q T

x y b b

T q

x y q

b b b

∂  ∂ 

∈ Ω  +  − =

∂ ∂

   

∈ Γ = =

∈ Γ = − ⋅ ∇ =

∈ Γ = − ⋅ ∇ n n

(4)

(3)

2 2 2

2 2

D D

2 ^y

x x

b y y b y y

v v v

T T T

x y x y y y y

 

∂ ∂ ∂

− −

∂ ∂ ∂ ∂ ∂ ∂ ∂

(6)

while

( )

^T

D D T

D D

T T

T v v v T

b b

⋅ ∇ = ⋅ ∇ + ⋅ ∇ ⋅ ⋅ ∇ ⋅ − ⋅ ∇ + ∇ ⋅ ∇

n n n n n n (7)

where

x y

v v

x x

v v v

y y

∂ ∂ 

 ∂ ∂ 

 

∇ = ∂ ∂ 

 ∂ ∂ 

 

(8)

In the case of expansion or contraction of the domain considered it is assumed that the velocities vx = x/b, vy = y/b. So

2 2 2

D D 2

D D

T T T

b x x b b x

∂  = ∂  − ∂

 ∂  ∂   ∂

  (9)

and

2 2 2

D D 2

D D

T T T

b y y b b y

∂  = ∂  − ∂

 ∂  ∂   ∂

  (10)

while

( )

D D 1

D D

T T

b b b T

⋅ ∇ = ⋅ ∇  − ⋅ ∇

 

n n n (11)

(4)

Introducing (9), (10), (11) into (4) and taking into account that from equation (1) results that

( )

2 2

λ λ

B m

T T

k T T Q

x y

∂ ∂

+ = − −

∂ ∂ (12)

one obtains:

( )

2 2

1

2

0 0

D D D 2 2

, : λ λ 0

D D D

, : D 0

D

, : λ D 0

D

, : λ D

D

T T T

x y k T Q

x b y b b b b

x y T

b x y T

b q x y T

b b

∂   ∂  

∈ Ω ∂  + ∂   − − + =

∈ Γ =

 

∈ Γ − ⋅ ∇  =

 

 

∈ Γ − ⋅ ∇  =

 

n n

(13)

where Q = kTB + Qm. Finally:

( )

2 2

1

2

0 0

2 2

, : λ λ 0

, : 0

, :

U U

x y k U k T Q

x y b b

x y U

x y W

x y W q

b

∂ ∂

∈ Ω + − − + =

∂ ∂

∈ Γ =

(14)

where U = DT/Db is the sensitivity function and

λ D λ

D

W T U

b

 

= − ⋅ ∇  = − ⋅ ∇

 

n n (15)

3. Boundary element method

In order to solve the basic problem (1) and additional problem (14) connected with the shape sensitivity analysis the boundary element method has been applied.

So, we consider the following equation

( ) ( )

( )

2 2

, ,

λ F x y λ F x y , 0

S x y

x y

∂ ∂

+ + =

∂ ∂ (16)

(5)

(ξ, η) is the observation point, F (ξ, η, x, y) is the fundamental solution

(

^{ξ, η, ,}

)

¹ ^ln¹

2 π λ

F x y

r

∗ = (18)

where r is the distance between the points (ξ, η) and (x, y)

(

^ξ

)

²

(

^η

)

²

r = x − + y − (19)

Function J ^∗(ξ, η, x, y) is defined as follows

(

^{ξ, η, ,}

)

^λ

(

^{ξ, η, ,}

)

J^∗ x y = − n⋅ ∇F^∗ x y (20)

and it can be calculated analytically

(

^{ξ, η , ,}

)

₂

2 π

J x y d

r

∗ = (21)

where

(

^{ξ co s α +}

) (

^{η c o s β}

)

d = x− y − (22)

In equation (17): J(x, y) = −λ n ⋅ ∇F(x, y).

In order to solve the equation (17) the boundary is divided into N boundary elements and the interior is divided into L internal cells. For constant boundary elements we assume that

( ) ( ) ( )

( ) ( )

, ,

, :

, ,

j j j

j

j j j

F x y F x y F

x y

J x y J x y J

 = =

∈ Γ 

= =

 (23)

and for constant internal cells

(

x y^,

)

∈ Ωl ^: S x y

(

^,

)

= S x y

(

l^, l

)

= Sl (24)

(6)

In this case the following approximation of equation (17) is obtained

( )

( ) ( )

1

1 1

1 ξ , η , , d

2

ξ , η , , d ξ , η , , d

j

j l

N

i j i i j

j

N L

j i i j l i i l

j l

F J F x y

F J x y S F x y

∗

= Γ

∗ ∗

= Γ = Ω

+ Γ =

Γ + Ω

∑ ∫

∑ ∫ ∑ ∫∫

(25)

where (ξ_i, η_i) denotes the boundary node and i = 1, 2,…, N.

The system of equations (25) can be written in the form

1 1 1

N N L

i j j i j j i l l

j j l

G J H F P S

= = =

= +

∑ ∑ ∑

⁽²⁶⁾

where

(

^{ξ , η , ,}

)

^d

j

i j i i j

G F^∗ x y

Γ

=

∫

Γ ⁽²⁷⁾

and

(

^{ξ , η , ,}

)

^d ^,

1, 2

j

i i j

i j

J x y i j

H

i j

∗

Γ

 Γ ≠

= 

 − =



∫

(28)

while

(

^{ξ , η , ,}

)

^d

l

i l i i l

P F^∗ x y

Ω

=

∫∫

Ω ⁽²⁹⁾

After solving the system of equations (26) the boundary values of functions Fj and J_j are known. The values F_i at the internal points (ξ_i, η_i)∈Ω are calculated using the formula

1 1 1

N N L

i i j j i j j i l l

j j l

F H F G J P S

= = =

=

∑

−

∑

+

∑

⁽³⁰⁾

It should be pointed out that in the case when the source function S is a function of T or U the basic algorithm should be supplemented by adequate iterative procedure [4].

(7)

illustrates the distribution of sensitivity function U.

Fig. 2. Temperature distribution

Fig. 3. Distribution of function U

(8)

Function T can be expanded into Taylor’s series taking into account two com- ponents

(

^{, ,}

) (

^{, ,}

) (

^{, ,}

)

T x y b+ ∆b =T x y b +U x y b ∆ b (31) In Figure 4 the temperature distribution obtained from equation (31) under the assumption that ∆b = 0.1b is shown.

Fig. 4. Temperature distribution for b +∆b

Summing up, the shape sensitivity analysis allows, among others, to estimate the changes of temperature in the case when the geometry of the domain considered is changed.

The paper has been sponsored by KBN (Grant No 3 T11F 018 26).

References

[1] Budman H., Shitzer A., Dayan J., J. of Biomechanical Engineering 1995, 117, 193.

[2] Kleiber M., Parameter sensitivity, J. Wiley & Sons Ltd., Chichester 1997.

[3] Dems K., Rousselet B., Structural Optimization 1999, 17, 36-45.

[4] Kałuża G., Doctoral thesis, Silesian University of Technology, Gliwice 2005.

[5] Majchrzak E., BEM in heat transfer, Publ. of the Techn. Univ. of Czest., Czestochowa 2001.