Application of the DRBEM for numerical solution of Cattaneo-Vernotte bioheat transfer equation

(1)

APPLICATION OF THE DRBEM FOR NUMERICAL SOLUTION OF CATTANEO-VERNOTTE BIOHEAT TRANSFER EQUATION

Ewa Majchrzak ^{1, 2}, GraŜyna KałuŜa ¹, Jolanta Poteralska ¹

1 Silesian University of Technology, Gliwice, Poland

2 Czestochowa University of Technology, Czestochowa, Poland ewa.majchrzak@pol.sl.pl

Abstract. The 1D problem of heat transfer in the tissue subjected to action of external heat

source is considered. This phenomenon is described by Cattaneo-Vernotte equation supplemented by adequate boundary and initial conditions. To solve the problem formulated the dual reciprocity boundary element method (DRBEM) is applied. In the final part of the paper the examples of computations are shown.

1. Formulation of the problem

According to the newest opinions the heat conduction proceeding in the biological tissue domain should be described by the hyperbolic equation (Cattaneo- -Vernotte equation [1, 2]) in order to take into account its nonhomogeneous inner structure. So the following bio-heat transfer equation is considered

2 2

( , ) ( , ) ( , ) ( , )

( , )

T x t T x t T x t Q x t

c Q x t

t t x t

 τ ∂ + ∂  = λ ∂ + + τ ∂

 

∂ ∂ ∂ ∂

  (1)

where c, λ denote the volumetric specific heat and thermal conductivity of tissue, Q (x, t) is the capacity of internal heat sources due to metabolism and blood perfusion, τ is the relaxation time (for biological tissue it is a value from the scope 20-35 s), T is the tissue temperature, x, t denote the spatial co-ordinates and time.

The function Q (x, t) is equal to

[ ]

( , )

_B _B _B

( , )

_m

Q x t = G c T − T x t + Q (2)

where G

B

is the blood perfusion rate, c

B

is the volumetric specific heat of blood, T

B

is the artery temperature and Q

m

is the metabolic heat source.

The equation (1) is supplemented by the boundary conditions

1

2

0 : ( , ) ( ) : ( , )

b

x T x t T t x L T x t T

= =

= = (3)

(2)

and initial ones

0

( , )

0 : ( , ) , 0

t

T x t t T x t T

t

₌

= = ∂ =

∂ (4)

where T

b1

(t), T

b2

are known boundary temperatures and T

0

is known initial temperature of biological tissue.

It should be pointed out that for τ = 0 the equation (1) reduces to the well- known Pennes bioheat equation.

2. Boundary element method

Taking into account formula (2) the equation (1) can be written in the form

[ ]

2 2

( , ) ( , )

( , )

m B B B B B

T x t T x t

c t t

T x t T x t

Q G c T T x t G c

x t

 τ ∂ + ∂  =

 

∂ ∂

 

∂ ∂

λ + + − − τ

∂ ∂

(5)

or

[ ]

2 2

( , ) ( , ) ( , )

( )

( , ) 0

B B

B B B m

T x t T x t T x t

c G c c

x t t

G c T T x t Q

∂ ∂ ∂

λ − + τ − τ +

∂ ∂ ∂

− + =

(6)

For t

= t

^f

one has

2 2

( , )

( , ) 0

f

T x t

f

S x t x

λ ∂ − =

∂ (7)

where

2 2

( , ) ( , )

( , ) ( )

( , )

f f

f

B B

t t t t

f

B B B m

T x t T x t

S x t c G c c

t t

G c T T x t Q

= =

∂ ∂

= + τ + τ −

∂ ∂

 −  −

 

(8)

Application of the standard boundary element method leads to the following

equation [3, 4]

(3)

* *

* * *

0

( , ) ( , ) ( , ) ( , 0) (0, ) ( , ) ( , ) ( , 0) (0, ) ( , ) ( , ) d

f f f

L

f f f

T t T L q L t T q t

q L T L t q T t S x t T x x

ξ + ξ − ξ =

ξ − ξ − ∫ ξ ⁽⁹⁾

where ξ is the observation point, T

³^*

(ξ, x) is the fundamental solution, q(x , t

^f

) =

−λ ∂ T (x, t

^f

) / ∂ x is the heat flux, q

^*

(ξ, x) =

− λ ∂ T

^*

(ξ, x) / ∂ x is the heat flux resulting from the fundamental solution.

For the problem considered the fundamental solution has following form

*

1 ( , ) ( )

T ξ x = 2 L − − ξ x

λ (10)

Heat flux resulting from the fundamental solution can be calculated analytically

*

( , ) sgn ( )

( , )

2 T x x

q x

x

∂ ξ − ξ

ξ = −λ =

∂ (11)

It should be pointed out that the function T

^*

(ξ, x) fulfills the equation

2 * 2

( , )

T x

x x

∂ ξ

λ = −δ ξ

∂ (12)

where δ (ξ, x) is the Dirac function.

3. Dual reciprocity boundary element method

The solution of Cattaneo-Vernotte equation (6) written for time t

^f

(Equation (7)) can be expressed as a sum [5, 6]

( ,

^f

) ˆ ( ,

^f

) ( ,

^f

)

T x t = T x t + U x t (13)

where the first function is the solution of Laplace equation

2 2

ˆ( , ) 0 T x t

f

x

λ ∂ =

∂ (14)

and U (x, t

^f

) is a particular solution of Equation (7), this means

[ ]

2 2

( , ) ( , ) ( , )

( )

( , ) 0

B B

B B B m

U x t U x t U x t

c G c c

x t t

G c T U x t Q

∂ ∂ ∂

λ − + τ − τ +

∂ ∂ ∂

− + =

(15)

(4)

From Equations (6), (13), (14), (15) results that

2 2

( , ) ( , )

( ) ( , )

( , ) ( , )

( ) ( , )

f

B B B B

t t

B B B B

t t

U x t U x t

c G c c G c U x t

t t

T x t T x t

c G c c G c T x t

t t

=

 + τ ∂ + τ ∂ +  =

 

∂ ∂

 

 + τ ∂ + τ ∂ + 

 

∂ ∂

 

(16)

So, the source function (8) can be expressed as follows

2 2

( , ) ( , )

( , ) ( )

( , )

f f

f

B B

t t t t

f

B B B m

U x t U x t

S x t c G c c

t t

G c T U x t Q

= =

∂ ∂

= + τ + τ −

∂ ∂

 −  −

 

(17)

It is generally difficult to find the solution U (x, t

^f

). In the dual reciprocity method, at first, the following approximation of function S (x, t

^f

) (c.f. Equation (8)) is proposed [5]

1

( , ) ( ) ( )

K

f f

k k

k

S x t a t P x

=

= ∑ ⁽¹⁸⁾

where a

k

(t

^f

) are unknown coefficients and P

k

(x) are approximating functions ful- filling the equations

2 2

( , ) ( )

f k k

U x t

P x x

= λ ∂

∂ (19)

In Equation (18), K corresponds to the total number of nodes, where 2 is the number of boundary nodes and K – 2 is the number of internal nodes.

The last integral in equation (9) can be expressed as follows

2

* *

2

0 0 1

( , ) ( , ) d ( ) ( ) ( , ) d

L L K

f f k

k k

U x

D S x t T x x a t T x x

=

x

= − ξ = − λ ∂ ξ

∑ ∂

∫ ∫ ⁽²⁰⁾

Integrating twice by parts with respect to x one obtains

2 * 2

1 0

*

1 0

( , )

( ) ( )

( ) ( , )

( ) ( , ) ( )

K L f

k k

k K x L

f k

k k

k x

T x

D a t U x dx

x

U x T x

a t T x U x

x x

=

= =

 ∂ ξ 

= −   λ ∂   −

 λ ξ ∂ − λ ∂ ξ 

 ∂ ∂ 

 

∑ ∫

∑

(21)

(5)

Taking into account the property (12) one has

* *

0 1

( ) ( ) ( , ) ( ) ( ) ( , )

K

f x L

k k k k x

k

D a t U T x W x U x q x

⁼₌

=

 

= ∑  ξ + ξ − ξ  ⁽²²⁾

where

( )

^k

( )

k

U x

W x x

= −λ ∂

∂ (23)

Finally, the Equation (9) can be written in the form

* * *

* *

1

* * *

( , ) ( , ) ( , ) ( , 0) (0, ) ( , ) ( , )

( , 0) (0, ) ( ) ( ) ( , ) ( )

( , 0) (0) ( , ) ( ) ( , 0) (0)

f f f f

K

f f

k k k

k

k k k

T t T L q L t T q t q L T L t

q T t a t U T L W L

T W q L U L q U

=

ξ + ξ − ξ = ξ

− ξ +   ξ + ξ −

ξ − ξ + ξ  

∑

(24) We define [5]

2 3

( ) 4 9

k k

k

x x x x

U x − −

= + (25)

from which results that

( ) ( ) 1

2 3

k

k k

x x

W x x x  − 

= λ −  + 

  (26)

On the basis of (25) the functions (19) are calculated 1 2

( ) 2 3

k k

x x

P x  − 

= λ  + 

  (27)

The following approximations of time derivatives appearing in formula (8) can be applied

( , ) ( , ) ( ,

1

)

f

f f

t t

T x t T x t T x t

t t

−

=

∂ = −

∂ ∆ (28)

2 1 2

2 2

( , ) ( , ) 2 ( , ) ( , )

( )

f

f f f

t t

T x t T x t T x t T x t

t t

− −

=

∂ = − +

∂ ∆ (29)

So, we have (c.f. Equation (18))

(6)

1

1 2

2

1

( , ) ( , )

( )

( , ) 2 ( , ) ( , )

( )

( , ) ( ) ( )

f f

B B

f f f

K

f f

B B B B B m k k

k

T x t T x t c G c

t T x t T x t T x t

c t

G c T x t G c T Q a t P x

−

− −

=

+ τ − +

∆

− +

τ +

∆

− − = ∑

(30)

or

1 2

1

( , ) ( , ) ( , ) ( ) ( )

K

f f f f

k k

k

AT x t BT x t

⁻

C T x t

⁻

D a t P x

=

+ + + = ∑ ⁽³¹⁾

where

2 2

2

, 2 ,

( ) ( )

( ) ,

B B B B

B B

B B B m

c G c c c G c c

A G c B

t t t t

C c D G c T Q

t

 

+ τ τ + τ τ

= + + = −  + 

∆ ∆  ∆ ∆ 

= τ = − −

∆

(32)

Let x

1

= 0, x

2

= L, and x

i

are the internal nodes, i = 3, 4,..., K. For these nodes the following equations are obtained

1 2

1

( ) , 1, 2, ...,

K

f f f f

i i i k k i

k

AT BT

⁻

C T

⁻

D a P x i K

=

+ + + = ∑ = ⁽³³⁾

This system of equations can be written in the matrix form

1 2

1 1 1 1 1 2 1 1

1 2

2 2 2 1 2 2 2 2

1 2

3 3 3 1 3 2 3 3

1 2

( ) ( ) ... ( )

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

...

( ) ( ) ... ( )

f f f

K

f f f

K

f f f

K

f f f

K K K K

K K K

AT BT C T D P x P x P x

P x P x P x

AT BT C T D

− −

 + + +  

 

+ + +

 

 + + +  =

 

+ + + 

 

 

1

2

3

...

f

f K

a a a

a

 

  

   

     

(34)

or

1 2

f f f f

A T + B T

⁻

+ C T

⁻

+ = D P a (35) from which results that

1

(

1 2

)

f

=

⁻

A

f

+ B

f⁻

+ C

f⁻

+ D

a P T T T (36)

(7)

The Equation (24) for the nodes x

i

, i

= 1, 2,..., K can be expressed in the form

1 2 1 2

1

( , ) (0, ) ( , ) (0, ) ( , )

( ) ( ) (0) ( ) (0) ( )

f f f f f

i i i i i

K f

k k i i k i k i k i k

k

T x t G q t G q L t H T t H T L t

a t U x G W G W L H U H U L

=

+ + = + +

 + + + + 

 

∑ ⁽³⁷⁾

or

1 1 2 2 1 1 2 2

1

( ) ( ) ( ) ( ) ( )

f f f f f

i i i i i

K

k k i i k i k i k i k

k

T G q G q H T H T

a U x G W x G W x H U x H U x

=

+ + = + +

 + + − − 

 

∑ ⁽³⁸⁾

where

1 2 1 2

, , 1

2 2 2

i i

i i i i

L x x

G = − − G = H = H =

λ λ (39)

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

( )

f

=

f

+ −

f

G q H T G W H U a (40)

where

11 12

21 22

31 32

1 2

0 ... 0 0 ... 0 0 ... 0 ... ... ... ... ...

0 ... 0

K K

G G

 

 

=  

 

 

G (41)

11 12

21 22

31 32

1 2

1 0 ... 0

1 ... 0 ... ... ... ... ...

0 ... 1

K K

H H

−

 

 

 − 

 

= −

 

 − 

 

H (42)

and

1 1 2 1 3 1 1

1 2 2 2 3 2 2

1 3 2 3 3 3 3

1 2 3

( ) ( ) ( ) ... ( )

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

( ) ( ) ( ) ... ( )

K

K K K K K

U x U x U x U x

 

 

=  

 

 

U (43)

(8)

1 1 2 1 3 1 1

1 2 2 2 3 2 2

( ) ( ) ( ) ... ( )

0 0 0 ... 0

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

0 0 0 ... 0

K

W x W x W x W x

 

 

=  

 

 

W (44)

while

1 1 1

2 2 2

3

, 0 ,

3

... ... ...

0

f f f

f f f f f

f f

K K

T q a

T a

     

     

= = =

     

 

   

T q a (45)

Putting (36) into (40) one has

1 2

( )

f f f f f

A B

⁻

C

⁻

D

= + + + +

G q H T M T T T (46)

where

( )

⁻1

= −

M G W H U P (47)

The boundary conditions (3) should be introduced to the system of equations (46) and next this system can be solved. A start point of numerical simulation process results from the initial conditions (4), in particular T

i0

= T

i1

= T

0

for i = 3, 4,... K.

4. Results of computations

The layer of tissue (L = 1 cm) limited by the skin surface and conventionally assumed internal one is considered. The following input data have been introduced λ = 0.5 W/(mK), c = 4.2 MJ/(m

³

K), c

B

= 3.9962 MJ/(m

³

K), G

B

= 0.002 1/s, T

B

= 37°C, Q

m

= 420 W/m

³

, τ = 35 s, T

0

= 37°C. The boundary condition on the skin surface has been assumed in the form T

b1

(t) = 37 + 0.25t, while for internal surface x

= L: T

b 2

= 37°C.

The problem has been solved by means of the DRBEM under the assumption that in the interior of domain 99 internal points have been distinguished and ∆t = 0.5.

In Figure 1 the temperature distribution in the tissue for times 10 and 20 s for

τ = 35 s (Cattaneo-Vernotte model) and τ = 0 (Pennes model) is shown. Figure 2

illustrates the heating curves at two points selected from the domain considered.

(9)

Fig. 1. Comparison of Cattaneo-Vernotte and Pennes models - temperature distribution

Fig. 2. Comparison of Cattaneo-Vernotte and Pennes models - heating curves for x = 0.0001 m and x = 0.001 m

The numerical solution of Cattaneo-Vernotte equation in comparison with the

Pennes one leads to the visible different results. Introduction of relaxation time

causes that the process proceeds slower.

(10)

Acknowledgement

This work was sponsored by Grant No N N501 366 734.

References

[1] Kaminski W., Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, Journal of Heat Transfer 1990, 112, 555-560.

[2] Liu J., Xu L.X., Boundary information based diagnostics on the thermal states of biological bodies, Journal of Heat and Mass Transfer 2000, 43, 2827-2839.

[3] Brebbia C.A., Dominguez J., Boundary elements, an introductory course, Computational Me- chanics Publications, McGraw-Hill Book Company, London 1992.

[4] Majchrzak E., Boundary element method in heat transfer, Publ. of the Czestochowa University of Technology, Czestochowa 2001 (in Polish).

[5] Partridge P.W., Brebbia C.A., Wrobel L.C., The dual reciprocity boundary element method, CMP, London, New York 1992.

[6] 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.