Temperature field in burned and healthy tissue - sensitivity analysis with respectto the thermal parameters

TEMPERATURE FIELD IN BURNED AND HEALTHY TISSUE - SENSITIVITY ANALYSIS WITH RESPECT

TO THE THERMAL PARAMETERS

Katarzyna Freus ¹ , Sebastian Freus ² , Ewa Majchrzak ³

1 Institute of Mathematics,Czestochowa University of Technology, Częstochowa, Poland

2 Institute of Computer and Information Science, Czestochowa University of Technology Częstochowa, Poland

3 Institute of Computational Mechanics and Engineering Silesian University of Technology, Gliwice, Poland

1 katarzyna.freus@im.pcz.pl, ² sebastian.freus@icis.pcz.pl, ³ ewa.majchrzak@polsl.pl

Abstract. A non-homogeneous system being the composition of burn wound and healthy tissue is considered. The heat exchange between sub-domains and environment is described by the system of partial differential equations (the Pennes equations) supplemented by the assumed boundary conditions. Additional problems associated with sensitivity analysis with respect to thermal parameters occurring in the mathematical model are formulated.

Both the basic problem and additional ones concerning the sensitivity with respect to selected parameters are solved using the boundary element method. In the final part of the paper the results of computations are shown.

Keywords: bioheat transfer, Pennes equation, boundary element method, burn wounds

Introduction

In the paper the heterogeneous domain being the composition of burned and healthy layers of skin tissue is considered. The temperature distribution in the domains is described by the system of two Pennes equations [1-4] with the differ- ent thermophysical parameters. In the healthy layer the metabolic and perfusion heat sources are taken into account, while the burned layer is dead and the blood perfusion and metabolism do not occur in this region [5]. The system of equations is supplemented by appropriate boundary conditions.

Thermophysical parameters occurring in the presented mathematical model differ significantly because they are an individual feature of a person. Thus, the aim of the research presented is to estimate the temperature changes due to changes in these parameters. So, the governing equations are differentiated with respect to the parameters considered [6-8]. In this way, additional problems are formulated, whose number corresponds to the number of the parameters analyzed. To solve the basic problem and additional ones the boundary element method is used [4, 9-11].

In the final part the results of computations are shown.

1. Governing equations

The domain of healthy tissue and burn wound, as shown in Figure 1, is consid- ered. The steady temperature field in burned tissue is described by the Laplace equation

2

1 : λ 1 1 ( ) 0

x ∈ Ω ∇ T x = (1)

where λ 1 is the thermal conductivity of burned tissue, T 1 (x) is the temperature and x = {x 1 , x 2 } are the spatial coordinates.

The temperature field in healthy tissue is described by the Pennes equation [1-3]

[ ]

2

2 : λ 2 2 ( ) _{B B B} 2 ( ) _met 0

x ∈ Ω ∇ T x + G c T − T x + Q = (2)

where T 2 (x) is the tissue temperature, λ 2 is the tissue thermal conductivity, G B is the blood perfusion rate, c B is the specific heat of blood, T B is the arterial blood temperature, Q met is the metabolic heat source.

Fig. 1. Domain considered

On the surface between sub-domains the continuity of heat flux and temperature field is assumed

1 2

1 2

1 2

( ) ( )

λ λ

:

( ) ( )

c

T x T x

n n

x

T x T x

∂ ∂

 − =

 ∂ ∂

∈ Γ 

 =



(3)

where ∂T e /∂n, e = 1, 2 denotes the normal derivative and n = [cosα 1 , cosα 2 ].

On the internal surface Γ in (cf. Fig. 1) the Dirichlet condition is taken into account

: 2 ( )

in b

x ∈ Γ T x = T (4)

Heat exchange between the domain considered and the environment is determined by the well-known Robin condition

[ ]

1 1 1

: λ ( ) α ( )

ex a

x T x T x T

n

∈ Γ − ∂ = −

∂ (5)

where α is the heat transfer coefficient, T a is the ambient temperature.

For the other boundaries the no-flux conditions are assumed

1 2 3 4

: λ _e T ^e ( ) x 0 x

n

∈ Γ ∪ Γ ∪ Γ ∪ Γ − ∂ =

∂ (6)

2. Sensitivity analysis

In this chapter the sensitivity analysis of the process discussed with respect to the thermophysical parameters appearing in mathematical model (1)-(6) is presented.

Let p 1 = λ 1, p 2 = λ 2 , p 3 = G B and p 4 = Q met . Equations (1), (2) are differentiated with respect to the parameter p s , s = 1, 2, 3, 4. So

2 2

1 1

1 1

λ ( ) ( )

( ) λ 0

s s

x T x

p T x p

 

∂ ∂

∇ + ∇   =

∂  ∂  (7)

and

[ ]

2 2

2 2 2

2 2 2

λ ( ) ( )

( ) λ ^B _{B B} ( ) _{B B} ^met 0

s s s s s

G Q

T x T x

T x c T T x G c

p p p p p

  ∂ ∂

∂ ∂ ∂

∇ + ∇   + − − + =

∂  ∂  ∂ ∂ ∂ (8)

Information resulting from equations (1), (2) is of the form

2

1 ( ) 0 T x

∇ = (9)

and

[ ]

2

2 2

2

( ) 1 ( )

λ ^{B B B met}

T x  G c T T x Q 

∇ = −  − +  (10)

Introducing (9) into equation (7) and (10) into equation (8), one has

2

1 1

λ ∇ U _s ( ) x = 0 (11)

and

2

2 2 2 2

λ ∇ U _s ( ) x − G c U _{B B s} ( ) x + R _s = 0 (12)

where

1 2

1 2

( ) ( )

( ) , ( )

s s

s s

T x T x

U x U x

p p

∂ ∂

= =

∂ ∂ (13)

are the sensitivity functions, and

[ ] [ ]

2

2 2 2

2

λ 1

( ) ( )

λ

B met

s B B B met B B

s s s

G Q

R G c T T x Q c T T x

p p p

∂ ∂

∂  

= − ∂  − +  + ∂ − + ∂ (14)

Next, the boundary condition (3) is differentiated with respect to p s

1 1 1 2 2 2

1 2

1 2

λ ( ) ( ) λ ( ) ( )

λ λ

:

( ) ( )

s s s s

c

s s

T x T x T x T x

p n n p p n n p

x

T x T x

p p

 ∂ ∂ ∂  ∂  ∂ ∂ ∂  ∂ 

− − = +

 ∂ ∂ ∂  ∂  ∂ ∂ ∂  ∂ 

    

∈ Γ 

∂ ∂

 =

 ∂ ∂



(15)

or

1 2

1 1 2 2

1 2

1 2

( ) ( )

λ ( ) λ ( )

λ +λ

:

( ) ( )

s s

s s

c

s s

U x U x

T x T x

p n n p n n

x

U x U x

∂ ∂

∂ ∂ ∂ ∂

 − − =

 ∂ ∂ ∂ ∂ ∂ ∂

∈ Γ 

 =



(16)

Boundary conditions (4), (5) and (6) are also differentiated, namely

: 2 ( ) / 0

in s b s

x ∈ Γ U x = ∂ T ∂ p = (17)

1

1 1

1 1

( )

λ ( )

: λ ^s α ( )

ex s

s

U x

x T x U x

p n n

∂

∂ ∂

∈ Γ − − =

∂ ∂ ∂ (18)

1 2 3 4

( ) ( )

: λ ^{e e} λ _e ^{e s} 0

s

T x U x

x p n n

∂ ∂

∈ Γ ∪ Γ ∪ Γ ∪ Γ − ∂ − =

∂ ∂ ∂ (19)

Equations (16), (18), (19) can be written in the form

1 1 1 1 1 1 2 2

1 2

1 / λ λ / ( ) ( ) 1 / λ λ / ( ) ( )

: ( ) ( )

s s s s

c

s s

p q x W x p q x W x

x U x U x

∂ ∂ + = − ∂ ∂ −

∈ Γ  

 = (20)

1

1 1 1

1

λ

: ( ) α ( ) 1 ( )

ex s s λ

s

x W x U x q x

p

∈ Γ = − ∂

∂ (21)

1 2 3 4

1 λ

: ( ) ( ) 0

λ

e

es e

e s

x W x q x

p

∂

∈ Γ ∪ Γ ∪ Γ ∪ Γ = − =

∂ (22)

In this way, the equations (11), (12) supplemented by boundary conditions (17), (20)-(22) create additional problems associated with sensitivity analysis with respect to the successive parameters p s .

3. Boundary element method

The basic problem and additional ones connected with the sensitivity analysis have been solved using the boundary element method [9-11]. The boundary inte- gral equation corresponding to the equation (1) is the following:

I I

1 1 1 1 1

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

B T q x T ^∗ x T x q ^∗ x

Γ Γ

= ∫ Γ + ∫ Γ ⁽²³⁾

where ξ is the observation point, the coefficient B(ξ) is dependent on the location of source point ξ, T ^* (ξ, x) is the fundamental solution, q 1 ^* ( ξ, x ) = − λ 1 ∂ T 1 ^* ( ξ, x ) / ∂ n is the heat flux resulting from fundamental solution, ^{q x} ₁ ( ) ^{= − λ} ₁ ^{∂ T x} ₁ ( ) ^/ ∂ is the ⁿ heat flux. The fundamental solution of the problem discussed is of the form

( )

* 1

1

1 1

ξ, ln

2πλ

T x

= r (24)

where r is the distance between the points ξ and x.

For healthy tissue (equation (2)) the boundary integral equation is as follows:

( ) ( ) ( ) ( ) ( ) ( ) ( )

II II

* * *

2 2 2 2 2 2

ξ ξ ξ, d ξ, d ξ, d

B T q x T x T x q x Q T x

Γ Γ Ω

+ ∫ Γ = ∫ Γ + ∫ Ω ⁽²⁵⁾

where Q = G _B c _B T B + Q met , ^q ₂ ^* ( ^{ξ, x} ) ^{= − λ} ₂ ^{∂ T} ₂ ^* ( ^{ξ, x} ) ^{/ ∂ n , q} ₂ ( ) ^{x = − λ} ₂ ^{∂ T} ₂ ( ) ^{x / ∂ n .}

The fundamental solution in the case discussed has a form

( )

*

2 0

2 2

ξ, 1 K

2πλ λ

G c B B

T x r

 

 

=  

 

(26)

where K ₀ (·) is the modified Bessel function of second kind and zero order [9, 10].

To solve the equations (23) and (25) the boundary Γ is divided into N elements Γ j = 1,2,…,N and the interior Ω ₂ is divided into L internal cells, as shown in Figure 2.

Next, the integrals in the equations (23), (25) can be replaced by the sums of inte- grals over these elements. So

1 1

1 1 1 1 1

1 1

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

j j

N N

j j

j j

B T q x T ^∗ x T x q x ^∗

= Γ = Γ

= ∑ ∫ Γ + ∑ ∫ Γ ⁽²⁷⁾

and

( ) ( ) ( ) ( )

( ) ( ) ( )

1

1

*

2 2 2

1

* *

2 2 2

1 1

ξ ξ ξ, d

ξ, d ξ, d

j

j l

N

j j N

N L

j l

j N l

B T q x T x

T x q x Q T x

= + Γ

= + Γ = Ω

+ Γ =

Γ + Ω

∑ ∫

∑ ∫ ∑ ∫

(28)

where N ₁ is the number of elements on the boundary limiting domain Ω ₁ .

After the mathematical manipulations one obtains the following systems of algebraic equations corresponding to the burned tissue

1 1 = 1 1

G q H T (29)

and healthy sub-domain

2 2 = 2 2 +

G q H T P (30)

The way of calculation of matrix G ₁ , H ₁ , G ₂ , H ₂ , P elements is described in detail in [10].

Fig. 2. Discretization of boundaries and interior Ω ₂

For the needs of further considerations concerning the temperature field computa- tions the following denotations are introduced (cf. Figs. 1 and 2)

−

1 2 1 2

1 , 1 , 1 , 1 , 1 , 1

ex ex

T T T q q q are the vectors of functions T and q at the boundary Γ ₁ ∪Γ 2 ∪Γ ^ex of domain Ω ₁ ,

− ₁ , ₂ , ₁ , ₂

c c c c

T T q q are the vectors of functions T and q on the contact surface Γ c between sub-domains Ω ₁ and Ω ₂ ,

− ₂ ³ , ₂ ⁴ , ₂ , ³ ₂ , ⁴ ₂ , ₂

in in

T T T q q q are the vectors of functions T and q at the boundary Γ ₃ ∪Γ 4 ∪Γ in of domain Ω ₂ .

The condition (3) can be written in the form

1 2 , 1 2

c = − c = c = c =

q q q T T T (31)

Using the above symbols, one obtains the following system of equations [11]

1 1

1 2 1

1 2

1 1 1 1 1 1 1

3 4

2 2 2 2 2 2

3 2

2 4 2

α α

ex

ex ex ex

c c a

in in

c c b

in

T T

 

 

 

 

 

 

 − − − −   

=

 

   

− − − −   +

   

   

 

 

 

 

 

 

T T

T

H G H H H G 0 0 0 T G

0 0 0 H G H G H q G P

T q T

(32) The system of equations (32) allows one to find the “missing” boundary values.

In the case of an additional problems solution (associated with the sensitivity anal- ysis) the procedure is similar. The systems of equations corresponding to equations (11), (12) and resulting from the boundary element method application have the following form for the burned tissue

1 1 s = 1 1 s

G W H U (33)

and healthy sub-domain

2 2 s = 2 2 s + 2 s

G W H U R (34)

or

1 1

1 1

1 2 1 1 2 1

1 1 1 1 2 1 1 1 1 2

1 1

1 1

s s

ex ex

ex s ex s

c c

s s

c s c s

   

   

   

    =    

   

   

   

   

W U

W U

G G G G H H H H

W U

W U

(35)

and

2 2

3 3

2 2

3 4 3 4

2 2 2 2 2 2 2 2 2

2 2

4 4

2 2

c s c s

s s

in in

c in c in s

s s

s s

   

   

   

    =     +

   

   

   

   

W U

W U

G G G G H H H H R

W U

W U

(36)

Taking into account the dependence (31) the boundary condition (20) can be written as

1 2

2 1

1 2

1 2

λ λ

1 1

: λ λ

c s c s

c s s

c s c s s

x p p

  ∂ ∂ 

= − − +

  

∈ Γ   ∂ ∂ 

 = =



W W q

U U U

(37)

Introducing (21), (37) into (35), (36) one obtains

1 1

1

1

1 1 1

1 2 1 2 1

1 1 1 1 1 1 1 1 1 2

2 1

1 1

λ 1 λ

s

s

ex ex ex

ex s ex s

s

c c

s s

s c s

p

 

 

 

 

 α − ∂   

 

  ∂ =    

     

 

 

 

   

 

W U

U q U

G G G G H H H H

W U W U

(38)

and

1 2

1

1 2

3 4 3

2 2 2 2 2

2 4 2

3

3 4 2

2 2 2 2 2

2 4 2

λ λ

1 1

λ λ

c s

s s

in

c s

in s

s

s in s

c in s

s

s

p p

  ∂ ∂  

− − +

   

∂ ∂

 

 

 

  =

   

 

 

 

 

 

 

 

    +

 

 

 

 

W q

G G G G W

W W U

H H H H U R

U U

(39)

After the introduction of remaining boundary conditions (17), (22) the systems of equations (38) and (39) take a form

1 1

1

1 2 2 1

1 1 1 1 1 1 1 1 1

1

1

λ 1 λ

s ex s

ex ex ex ex

c c s

s s

c s

p

 

 

 

  ∂

 − α − − −  =

    ∂

 

 

 

 

U U

H G H H H G U G q

U W

(40)

and

1

3 4 3 1 2

2 2 2 2 2 2 2 2

1 2

2 4 2

λ λ

1 1

λ λ

s c s in

c c s c s

s s

in s s

p p

 

 

 

 ∂ ∂ 

 

 − − − −  =  +  +

       ∂ ∂ 

 

 

 

U W

H G H G H U G q R

W W

(41)

Equations (40), (41) coupling gives

1 1

1 2 1

1 2

1 1 1 1 1 1

3 4

2 2 2 2 2 1

3 2 2 4 2

1

1 1

1

1 2

2 2

1 2

α

λ 1 λ

λ λ

1 1

λ λ

s ex

s

ex ex

c c s

in c s

c c

s in s

s ex ex

s

c s

s s

p

p p

 

 

 

 

 

 

 − − − − 

  =

 

− − − −  

 

 

 

 

 

 

 

 

 ∂ 

 ∂ 

 

  ∂ ∂  

+ +

   

∂ ∂

   

 

U U U

H G H H H G 0 0 0 U

0 0 0 H G H G H W

U W

U G q

G q R

(42)

It should be pointed out that the main matrix of the system of equations (42) asso- ciated with the sensitivity functions is the same as in the case of the basic problem solution (cf. equation (32)).

4. Results of computations

The domain of dimensions 0.04 m × 0.02 m has been considered. The following input data have been assumed: thermal conductivity of burned tissue λ ₁ = 0.1 W/(mK), thermal conductivity of healthy tissue λ ₂ = 0.2 W/(mK) [5], blood perfusion rate W _B = 0.5 kg/(m ³ s), specific heat of blood c _B = 4200 J/(kgK), arterial blood tempera- ture T _B = 37°C, metabolic heat source Q _met = 200 W/m ³ , heat transfer coefficient α = 10 W/(m ² K), ambient temperature T _a = 20°C.

Figure 3 illustrates the temperature distribution in the domain considered, while

Figures 4-7 show the distributions of sensitivity functions with respect to the suc-

cessive parameters.

Fig. 3. Temperature distribution

Fig. 4. Distribution of sensitivity function with respect to the parameter p ₁ = λ 1

Fig. 5. Distribution of sensitivity function with respect to the parameter p ₂ = λ 2

Fig. 6. Distribution of sensitivity function with respect to the parameter p ₃ = G _B

Fig. 7. Distribution of sensitivity function with respect to the parameter p ₄ = Q met

Conclusions

The application of sensitivity analysis methods gives essential information con- cerning the influence of thermophysical parameters perturbations on the changes of temperature field on the homogeneous or heterogeneous tissue domain. In this way one can estimate which parameter significantly affects the final result and is insig- nificant in regards to which one. For example, the results presented above show the visible influence of thermal conductivity of burned tissue, while on the other hand the value of metabolic heat source capacity has a little importance for the tempera- ture distribution. Such knowledge can be useful at the stage of bioheat problems modeling.

References

[1] Pennes H.H., Analysis of tissue and arterial blood temperatures in the resting human forearm, Journal of Applied Physiology 1948, 1, 93-122.

[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., Paruch M., Identification of electromagnetic field parameters assuring the cancer destruction during hyperthermia treatment, Inverse Problems in Science & Engineering 2011, 19, 1, 45-58.

[4] Majchrzak E., Application of different variants of the BEM in numerical modeling of bioheat transfer processes, MCB: Molecular & Cellular Biomechanics 2013, 10, 3, 201-232.

[5] Romero Mendez R., Jimenez-Lozano J.N., Sen M., Gonzalez F.J., Analytical solution of a Pennes equation for burn-depth determination from infrared thermographs, Mathematical Medicine and Biology 2010, 27, 21-38.

[6] Kleiber M., Parameter Sensitivity, J.Wiley & Sons Ltd., Chichester 1997.

[7] Majchrzak E., Mochnacki B., Sensitivity analysis and inverse problems in bio-heat transfer mod- elling, Computer Assisted Mechanics and Engineering Sciences 2006, 13, 85-108.

[8] Jasiński M., Sensitivity analysis of transient bioheat transfer with perfusion rate dependent

on tissue injury, Computer Assisted Mechanics and Engineering Science 2009, 16, 267-277.

[9] Brebbia C.A., Dominguez J., Boundary Elements an Introductory Course, CMP, McGraw-Hill Book Company, London 1992.

[10] Majchrzak E., Boundary Element Method in Heat Transfer, Publ. of the Techn. Univ. of Czest., Czestochowa 2001 (in Polish).

[11] Freus K., Freus S., Majchrzak E., Determination of the temperature field in burned and healthy

tissue using the boundary element method - part I, Journal of Applied Mathematics and Compu-

tational Mechanics 2013, 12(3), 39-46.