• Nie Znaleziono Wyników

Heating of tissue by means of the electric field - numerical model basing on the BEM

N/A
N/A
Protected

Academic year: 2022

Share "Heating of tissue by means of the electric field - numerical model basing on the BEM"

Copied!
12
0
0

Pełen tekst

(1)

HEATING OF TISSUE BY MEANS OF THE ELECTRIC FIELD - NUMERICAL MODEL BASING ON THE BEM

Ewa Majchrzak1, 2, Joanna Drozdek2, Marek Paruch1

1 Department for Strength of Materials and Computational Mechanics Silesian University of Technology, Poland

2 Institute of Mathematics, Czestochowa University of Technology, Poland ewa.majchrzak@polsl.pl, jdrozdek@imi.pcz.pl, marek.paruch@polsl.pl

Abstract. The domain of tissue is subjected to the action of electrodes located on the skin surface. External electric field causes the heat generation in tissue domain. The distribution of electric potential in domain considered is described by the Laplace equation, while the tem- perature field is described by the Pennes equation. These problems are coupled by source function being the additional component in Pennes equation and resulting from the electric field action. The coupled problem is solved using the boundary element method. In the final part of the paper the examples of computations are shown.

1. Governing equations

In Figure 1 a typical radio frequency (RF) hyperthermia system is shown [1]. The mathematical model of the process analyzed consists of two parts [1, 2]. The electric part concerns the Laplace equation to obtain the electric field distribution. The ther- mal part is connected with the bioheat transfer equation to obtain the temperature distribution. In the bioheat transfer equation the additional source term associated with the heat generation caused by electric field distribution appears.

Fig. 1. Action of electric field on biological tissue

(2)

The potential inside the tissue is described by the Laplace equation

(

x y,

)

∈Ω ∇: ε

(

x y,

) (

φ x y,

)

=0 (1) where ε(x, y) [C2 /(Nm2 )] is the dielectric permittivity of tissue.

On the external surface of tissue being in a contact with the electrodes the following condition is accepted:

( ) ( )

( )

12

( )

, : φ ,

, : φ ,

x y x y U

x y x y U

∈Γ =

∈Γ = − (2)

where U [V] is the electric potential of the electrode relative to the ground.

On the remaining external boundary of the tissue the ideal electric isolation is as- sumed:

( )

φ ,

ε x y 0

n

− ∂ =

∂ (3)

The electric field inside the tissue is described by equation

( ) ( )

( )

( )

φ ,

, φ ,

φ , x y x y x y x

x y y

∂ 

 

 ∂ 

= −∇ = −

∂ 

 

 ∂ 

E (4)

Heat generation Q [W/m3 ] due to the electromagnetic dissipated power in tissue depends on the conductivity σ [S/m] and the electric field E [2]

(

,

)

σ

(

,

)

σ φ

(

,

)

2 φ

(

,

)

2

2 2

x y x y x y

Q x y

x y

∂  ∂  

 

= =   + 

∂ ∂

    

 

E (5)

The temperature field in the domain considered is described by the Pennes equation [1, 2]

(

x y,

)

∈Ω: λ2T x y

(

,

)

+G cB BTBT x y

(

,

)

+Qmet+Q x ye

(

,

)

=0 (6) where T denotes the temperature, λ [W/(mK)] is the thermal conductivity, GB [1/s] is the perfusion rate, cB [J/(m3K)] is the volumetric specific heat of blood, TB is the supplying arterial blood temperature which is treated as a constant, Qmet is the metabolic heat source.

On the upper and lower surfaces of tissue domain (skin surface) the Dirichlet condition is assumed T(x, y)=Tb, where Tb is known temperature, on the remaining internal boundaries of tissue the adiabatic condition can be taken into account:

−λ∂T(x, y)/∂n= 0.

(3)

2. Boundary element method - electric field

The boundary integral equation corresponding to the equation (1) is following [3]

( ) ( ) ( ) ( )

( ) ( )

*

*

ξ, η φ ξ, η ψ , φ ξ, η, , d φ , ψ ξ, η, , d

B x y x y

x y x y

Γ

Γ

+ Γ =

Γ

(7)

where (ξ, η) is the observation point, the coefficient B(ξ, η) is dependent on the location of source point (ξ, η), ψ

( )

x y, = − ∂ε φ

( )

x y, /n.

Fundamental solution of the problem discussed has the following form

( )

* 1 1

φ ξ, η, , ln

x y 2πε

= r (8)

where r is the distance between points

( )

ξ, η and

(

x y . Differentiating the func-,

)

tion φ ξ, η, ,*

(

x y with respect to the outward normal

)

n=

[

cos α, cosβ

]

the func- tion ψ ξ, η, ,*

(

x y is obtained

)

( )

*

( )

*

2

φ ξ, η, ,

ψ ξ, η, , ε

x y d

x y n r

= − ∂ =

∂ (9)

where

(

ξ cos α

) (

η cosβ

)

d = −x + y− (10)

The boundary of the domain is divided into N boundary elements. For constant boundary elements it is assumed that

( ) ( ) ( )

( ) ( )

φ , φ , φ

, :

ψ , ψ , ψ

j j j

j

j j j

x y x y

x y

x y x y

 = =

∈Γ 

= =

 (11)

and then one obtains the following approximation of Equation (7)

1 1

ψ φ , 1, 2,...,

N N

ij j ij j

j j

G H i N

= =

= =

∑ ∑

(12)

where

( )

φ ξ , η , ,* d

j

ij i i j

G x y

Γ

=

Γ (13)

and (i ≠ j)

(4)

( )

ψ ξ , η , ,* d

j

ij i i j

H x y

Γ

=

Γ (14)

while Hi i = –1/2.

The system of Equations (12) can be written in the matrix form

=

(15)

This system allows to determine the 'missing' boundary values of functions φ , ψj j. Next, the values of function ϕ at the internal points

(

ξ , ηi i

)

can be determined using the formula

1 1

φ φ ψ , 1, 2,...,

N N

i ij j ij j

j j

H G i N N N L

= =

=

= = + + + (16)

It should be pointed out that in order to determine the electric field inside tissue (Eq- uation (5)) the partial derivatives ∂φe

(

x y,

)

/x, ∂ψe

(

x y,

)

/y must be known.

One of the possibilities is application of equation (7) for internal nodes

( )

ξ, η

(

B

( )

ξ, η =1

)

and then

( ) ( )

*

( ) ( ) (

*

)

φ ξ, η ψ ξ, η, , φ ξ, η, ,

φ , d ψ , d

ξ ξ η

x y x y

x y x y

Γ Γ

∂ ∂ ∂

= Γ − Γ

(17)

and

( ) ( )

*

( ) ( ) (

*

)

φ ξ, η ψ ξ, η, , φ ξ, η, ,

φ , d ψ , d

η η ξ

x y x y

x y x y

Γ Γ

∂ ∂ ∂

= Γ − Γ

(18)

where

* *

2 2

φ ξ φ η

ξ 2πε , η 2πε

x y

r r

∂ = − ∂ = −

∂ ∂ (19)

and

( ) ( )

* *

4 2 4 2

2 ξ 2 η

ψ 1 cos α ψ 1 cosβ

ξ 2π , η 2π

x d y d

r r r r

 −   − 

∂ =  −  ∂ =  − 

∂   ∂   (20)

Applying previously presented discretization of the boundary of domain, numeri- cal calculations of partial derivatives are not difficult to obtain. These derivatives are determined at the internal nodes.

,

(5)

3. Dual reciprocity boundary element method - temperature field The Pennes equation (6) can be written in the form

(

x y,

)

∈Ω: λ2T x y

(

,

)

gT x y

(

,

)

+Q x y

(

,

)

=0 (21)

where

( ) ( )

, , ,

B B B met e

g=G c Q x y =gT +Q +Q x y (22) The standard boundary element method algorithm leads to the following integral equation [2, 3]

( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( )

*

* *

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

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

B T q x y T x y

T x y q x y Q x y gT x y T x y

Γ

Γ

+ Γ =

Γ +  −  Ω

∫ ∫

(23)

where

( )

* 1 1

ξ, η, , ln

T x y 2πλ

= r (24)

and

( )

*

( )

*

2

ξ, η, ,

ξ, η, , λ

n 4π

T x y d

q x y

r

= − ∂ =

∂ (25)

while q x y

(

,

)

= − ∂λ T x y

(

,

)

/n.

It should be pointed out that the function T*

(

ξ, η, ,x y fulfills the equation

)

( ) ( )

2 *

λ∇T ξ, η, ,x y = −δ ξ, η, ,x y (26) where δ ξ, η, ,

(

x y is the Dirac function.

)

The solution of Pennes equation (21) can be written as a sum

(

,

) (

ˆ ,

) (

,

)

T x y =T x y +U x y (27)

where the first function is the solution of Laplace's equation

( )

2ˆ

λ∇T x y, =0 (28)

and U(x, y) is a particular solution of equation

(6)

( ) ( ) ( )

λ∇2U x y, −gU x y, +Q x y, =0 (29) From Equations (27), (28), (29) results that

(

,

) (

,

)

T x y =U x y (30)

In the dual reciprocity method the following approximation is proposed [4]

( ) ( ) ( )

1

, , ,

N L k k k

Q x y gU x y a f x y

+

=

− = −

(31)

where ak are unknown coefficients and fk(x, y) are approximating functions fulfilling the equations (c.f. Equation (29))

( ) ( )

λ∇2Uk x y, = fk x y, (32)

Putting (32) into (31) one obtains

( ) ( )

2

( )

1

, , λ ,

N L

k k

k

Q x y gU x y a U x y

+

=

− = −

(33)

In Equations (31), (33) N+L corresponds to the total number of nodes, where N is the number of boundary nodes and L is the number of internal nodes.

We consider the last integral in equation (23)

(

,

) (

,

) (

* ξ, η, ,

)

d

D Q x y gT x y T x y

=

 −  Ω (33)

Taking into account the dependences (30), (33) one obtains

( ) ( )

2 *

1

λ , ξ, η, , d

N L

k k

k

D a U x y T x y

+

= Ω

 

= −

∑ ∫

 ∇  Ω (34)

Using the second Green formula [3] one has

( ) ( )

( ) ( )

( ) ( )

2 *

1

*

1

*

1

λ ξ, η, , , d

λ ξ, η, , , d

ξ, η, ,

λ , d

N L

k k

k N L

k k

k N L

k k

k

D a T x y U x y

U x y

a T x y

n

T x y

a U x y

n

+

=

+

= Γ

+

= Γ

 

= −  ∇  Ω −

∂ Γ +

∂ Γ

∑ ∫

∑ ∫

∑ ∫

(35)

Because (c.f. formula (26))

(7)

( ) ( )

( ) ( ) ( ) ( )

2 *

λ ξ, η, , , d

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

k

k k

T x y U x y

x y U x y B U

 ∇  Ω =

 

− Ω = −

∫∫

(36)

so

( ) ( )

( ) ( ) ( ) ( )

1

* *

1

ξ, η ξ, η

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

N L

k k

k N L

k k k

k

D a B U

a T x y W x y U x y q x y

+

= +

= Γ

= +

 −  Γ

 

∑ ∫

(37)

where

(

,

)

λ k

(

,

)

k

U x y W x y

n

= − ∂

∂ (38)

Taking into account the formula (38) the Equation (23) can be written in the form

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

*

*

1

* *

1

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

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

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

N L

k k

k N L

k k k

k

B T T x y q x y

q x y T x y a B U

a T x y W x y q x y U x y

Γ

+

Γ = +

= Γ Γ

+ Γ =

Γ + +

 

Γ − Γ

 

 

∫ ∑

∑ ∫ ∫

(39)

For constant boundary elements the following approximation of equation (40) can be taken into account (i = 1, 2, ...,N, N+1,..., N+L)

1 1 1 1 1

ˆ ˆ

N N N L N N

i i ij j ij j k i ik ij jk ij jk

j j k j j

B T P q R T a B U P W R U

+

= = = = =

 

+ = +  + − 

 

∑ ∑ ∑ ∑ ∑

(40)

where

( )

* ξ , d

j

ij i j

P T x

Γ

=

Γ (41)

and

( )

ˆ * ξ , d

j

ij i j

R q x

Γ

=

Γ (42)

while Bi=B

(

ξ , ηi i

)

. We define [5, 6]

(8)

2 3

4 9

jk jk

jk

r r

U = + (43)

where

( ) (

2

)

2

2

jk k j k j

r = xx + yy (44)

Using the formula (39) one obtains

1 1

λ cos α cosβ λ

2 3

jk

j

jk j j jk jk

jk

j

U x

W d r

U y

 

 

∂  

 

 

= −   ∂ =  + 

 

 

(45)

where

( )

cos α

( )

cosβ

jk k j j k j j

d = xx + yy (46)

Because

2 2

2

2 2

1 2

sk sk 1

sk sk

k k

U U

U r

x x

∂ ∂

∇ = + = +

∂ ∂ (47)

so on the basis of equation (32) one has

(

,

) (

λ 1

)

sk k s s sk

f = f x y = +r (48)

Taking into account the dependencies (30), (31) one obtains

1 N L

s s k sk

k

gT Q a f

+

=

− =

(49)

where Ts =T x y

(

s, s

)

and Qs =Q x y

(

s, s

)

. The system of Equations (50) can be writ- ten in the matrix form

11 12 1,

1 1 1

,1 ,2 ,

,1 , ,

...

... ... ... ...

... ...

...

... ... ... ...

... ...

...

N L

N N N N L

N N N

N L N L L N L N L

N L N L N L

f f f

gT Q a

f f f

gT Q a

f f f

gT Q a

+

+

+ + + +

+ + +

−  

   

 

   

 

   

 

 − =  

 

   

 

   

 

 −   

    

(50)

or

( )

gT Q− = → =fa a f1 gT Q − (51)

(9)

The following matrices of dimensions N+L×N+L can be defined

11 12 1,

,1 ,2 ,

1,1 1,2 1,

,1 ,2 ,

... 0 ... 0

... ... ... ... ... ... ...

... 0 ... 0

... 0 ... 0

... ... ... ... ... ... ...

... 0 ... 0

N

N N N N

N N N N

N L N L N L N

P P P

P P P

P P P

P P P

+ + +

+ + +

 

 

 

 

= 

 

 

 

 

 

P (52)

11 12 1,

,1 ,2 ,

1,1 1,2 1,

,1 ,2 ,

... 0 ... 0

... ... ... ... ... ... ...

... 0 ... 0

... 1 ... 0

... ... ... ... ... ... ...

... 0 ... 1

N

N N N N

N N N N

N L N L N L N

R R R

R R R

R R R

R R R

+ + +

+ + +

 

 

 

 

= 

 − 

 

 

 

 

R (53)

where

ˆ , ˆ 0.5,

ij ij

ij

R i j

R

R i j

 ≠

=

− =

 (54)

and

11 12 1, 1, 1 1,

,1 ,2 , , 1 ,

1,1 1,2 1, 1, 1 1,

,1 ,2 , , 1 ,

... ...

... ... ... ... ... ... ...

... ...

... ...

... ... ... ... ... ... ...

... ...

N N N L

N N N N N N N N L

N N N N N N N N L

N L N L N L N N L N N L N L

U U U U U

U U U U U

U U U U U

U U U U U

+ +

+ +

+ + + + + + +

+ + + + + + +

 

 

 

 

= 

 

 

 

U



(55)

11 12 1, 1, 1 1,

,1 ,2 , , 1 ,

... ...

... ... ... ... ... ... ...

... ...

0 0 ... 0 0 ... 0

... ... ... ... ... ... ...

0 0 ... 0 0 ... 0

N N N L

N N N N N N N N L

W W W W W

W W W W W

+ +

+ +

 

 

 

 

= 

 

 

 

 

 

W (56)

So, the system of Equations (41) can be written in the matrix form

(10)

( )

= + −

Pq RT PW RU a (57)

or (c.f. formula (52))

( ) (

1 g

)

= + − −

Pq RT PW RU f T Q (58)

where

1

...

0 ...

0

N

q

q

 

 

 

 

= 

 

 

 

 

 

q (59)

This system of equations allows to find, among others, the temperatures at the boundary and internal nodes.

4. Results of computations

The rectangular domain of dimensions 0.08m×0.04m has been considered. The heating area is described as {0.032 ≤ x ≤ 0.048, y = 0 m}, {0.032 ≤ x ≤ 0.048, y = 0.04m} and the voltage applied on these surfaces is 10V and −10V, respec- tively. For biological tissue the following parameters have been assumed: thermal conductivity λ = 0.5 W/(mK), perfusion rate GB = 0.0005 1/s, metabolic heat source Qmet = 420 W/m3, blood temperature TB = 37°C, volumetric specific heat of blood cB = 4.2 MJ/(m3K) [2]. On the skin surface the temperature T = 32.5°C has been accepted. At first, the temperature distribution in the tissue without electric field influence under the assumption that on the external surface 60 constant boundary elements have been distinguished (Fig. 2).

Fig. 2. Temperature distribution without electric field

(11)

The distibution of electric field is shown in Figure 3, while Figure 4 illustrates the temperature field in the tissue subjected to electric field. The source function Qe under the assumption that σ = 0.4 S/m (Equation (5)) is shown in Figure 5.

Fig. 3. Electric field distribution

Fig. 4. Temperature distribution

Fig. 5. Source function due to the electric field

(12)

Summing up, the boundary element method has been applied to solve the cou- pled problem connected with the biological tissue heating. The simplified 2D mathematical model based on the Pennes equation supplemented with an equation determining the electric field due to the external electrodes action has been con- sidered.

Acknowledgement

This work was supported by Grant No N N501 3667 34 from the Polish Ministry of Education and Science.

References

[1] Lv Y.G., Deng Z.S., Liu J., 3D numerical study on the induced heating effects of embedded micro/nanoparticles on human body subject to external medical electromagnetic field, IEEE Transactions on Nanobioscience 2005, 4, 4, 284-294.

[2] Majchrzak E., Dziatkiewicz G., Paruch M., The modelling of heating a tissue subjected to exter- nal electromagnetic field, Acta of Bioengineering and Biomechanics 2008, 10, 2, 29-37.

[3] Majchrzak E., Boundary element method in heat transfer, Publ. of the Czestochowa University of Technology, Czestochowa 2001.

[4] Drozdek J., Majchrzak E., Numerical solution of bioheat transfer equation by means of the dual reciprocity method, Scientific Research of the Institute of Mathematics and Computer Science of Czestochowa University of Technology 2007, 1(6), 47-56.

Cytaty

Powiązane dokumenty

Dla odróżnienia, która z form narratywnych powinna być oddana w języku polskim odpowiednikiem leksykalnym, a w której zachodzi neutralizacja modalności

The estimated heat exchange area in each effect equals the heat exchange area of a single effect evaporator with the same evaporation capacity as the whole multiple

When the extraction voltage is applied, the surface chemi- cal potential is no longer a function of the local slope alone but also of the energy density of the local electric field:

Do pierw szej k ateg o rii zaliczają się te, k tó re są łatw o rozszyfrow alne, zrozum iałe dla każdego czytelnika in teresu jąceg o się lite ra tu rą... Nagi,

Każdy z tych sposobów pisania, jakkolw iek nie podporząd­ kow uje się całkowicie żadnej z obficie wykorzysty­ wanych gatunkowych konw encji, niem niej na oso­

In the post-Riverdance era, Venable (2001) states there has been a global growth of Irish dance that continues to invent tradition. 286), “One of the most ironic

(Received 30 January 2017; revised manuscript received 23 May 2017; published online 30 June 2017) The quantum mechanical theory of spectral parameters and dynamic conductivity

Prowincjał polskiej prowincji pasjonistów dr Łukasz Andrzejewski CP, przed- stawił referat Wkład Ojca Profesora Damiana Henryka Wojtyski CP do histo-... 543