Sensitivity analysis of burn integrals with respect to thickness of epidermis

SENSITIVITY ANALYSIS OF BURN INTEGRALS WITH RESPECT TO THICKNESS OF EPIDERMIS

Marek Jasiński

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

Abstract. In the paper the numerical analysis of thermal process proceeding in the domain of one-dimensional skin tissue subjected to an external heat source is presented. The degree of the skin burn can be predicted on the basis of Henriques integrals. Main subject of paper is the sensitivity analysis of these integrals with respect to the thicknesses of epidermis and dermis. On the stage of numerical realization the boundary element method has been used.

1. Governing equation

The skin is treated as a multilayer domain, in which one can distinguish the following sub-domains: epidermis Ω1 of thickness L1 − L0 [m], dermis Ω2 of thick- ness L² − L¹ and sub-cutaneous region Ω³ of thickness L³ − L² - Figure 1.

Fig. 1. Skin tissue domain

The transient bioheat transfer in the domain of skin is described by the follow- ing system of Pennes equations [1]:

2

: ^e 2^e ( )

e e e e B e me

T T

x c k T T Q

t x

∂ ∂

∈ Ω = λ + − +

∂ ∂ (1)

where e identifies the epidermis, dermis and sub-cutaneous region, λe [W/mK]

is the thermal conductivity and c_e [J/m³K] is specific heat per unit of volume, ke = GecB is the product of blood perfusion rate and volumetric specific heat of blood, T_B is the blood temperature and Q_me is the metabolic heat source. It should be pointed out that for the epidermis sub-domain (e = 1) G¹ = 0 and Qm1 = 0.

On the contact surfaces between sub-domains considered the continuity conditions are given, namely

1 1 , 1

1

: , 1, 2

e e

e e

e e

e e

T T

x x x e

T T

+ + +

+

∂ ∂

−λ = λ

∈ Γ  ∂ ∂ =

 =



(2)

Additionally

0 0

0

0

: ,

( ),

q t t

x q

T T^∞ t t

 ≤

∈ Γ = α − > (3)

where q = λ1∂T1/∂x, q0 is the given boundary heat flux, t0 is the exposure time, α is the heat transfer coefficient, T ^∞ is the ambient temperature. For conventionally assumed boundary limiting the system the no-flux condition

3: 3 0

x∈ Γ q = (4)

can be accepted. For t = 0 the initial temperature distribution is known, namely

1 1 2 2 3 3

0 : _p( ), _p( ), _p( )

t= T =T x T =T x T =T x (5)

A quadratic initial temperature distribution between 32.5°C at the surface and 37°C at the base of the sub-cutaneous region was introduced [1].

Thermal damage of skin begins when the temperature at the basal layer (the interface between epidermis and dermis) rises above 44°C (317 K). Henriques [2] found that the degree of skin damage could be predicted on the basis of the integrals

0

( ) exp d

b b b ( )

b

I P T E t

RT t

τ  ∆ 

= − 

 

∫

⁽⁶⁾

and

0

( ) exp d

d d d ( )

d

I P T E t

RT t

τ  ∆ 

= − 

 

∫

⁽⁷⁾

where ∆E/R [K] is the ratio of activation energy to universal gas constant, Pb, Pd

[1/s] are the pre-exponential factors, while T_b, T_d [K] are temperatures of basal

layer (the surface between epidermis and dermis) and dermal base (the surface between dermis and sub-cutaneous region).

First degree burn are said to occur when the value of the burn integral (6) is from the interval 0.53 < I_b ≤ 1, while the second degree burn when I_b > 1 [1, 2]. Third degree burn are said to occur when Id > 1. So, in order to determine the values of integrals I_b, I_d the heating and next the cooling curves for the basal layer and der- mal base must be known.

2. Shape sensitivity analysis

In the paper [3] the sensitivity analysis of temperature field in domain of skin tissue with respect to thermophysical parameters of skin has been presented. Here, the shape sensitivity analysis of temperature distribution and burn integrals is dis- cussed. Similar problem has been presented in [4], but only the sensitivities of burn integrals with respect to shape design parameter L0 have been calculated.

Here, the modification of parameter L1 is also discussed.

Using the concept of material derivative we can write [5, 6]

D D

e e e

s

s s

T T T

v

b b x

∂ ∂

= +

∂ ∂ (8)

where vs = vs(x, bs) is the velocity associated with design parameter b¹ = L⁰ or b2 = L1.

Forthematerialderivativefollowing formulascan bederived (c.f.equation(8))[4]:

D D

D D

e e e s

s s

T T T v

b x x b x x

 

∂ ∂ ∂ ∂

 =  −

 ∂  ∂ ∂ ∂

    (9)

D D

D D

e e

s s

T T

b t t b

 

∂ ∂

 =  

 ∂  ∂

    (10)

2 2 2 2

2 2 2 2

D

D 2

D D

e e e es e es

s s

T T T v T v

b x x b x x x x

 

∂  ∂ ∂ ∂ ∂ ∂

=  − −

 

∂ ∂ ∂ ∂ ∂ ∂

    (11)

If the direct approach of sensitivity method is applied [4-6] then the equations (1) are differentiated with respect to parameters bs, s = 1, 2.

Introducing the functions U_es = DT_e/Db_s and using the formulas (9), (10), (11) one has

2 2 2

2 2 2

: ^{es es} 2 ^{e s e s}

e e e e es

U U T v T v

x c k U

t x x x x x

 

∂ ∂ ∂ ∂ ∂ ∂

∈ Ω = λ  − − −

∂  ∂ ∂ ∂ ∂ ∂  (12)

or (c.f. equation (1))

2 2

2 2 ( ) 2

es es e s e s

e e e es e e B e me e

U U T v T v

c k U c k T T Q

t x t x x x

∂ ∂  ∂ ∂ ∂ ∂

= λ − −  − − −  − λ

∂ ∂  ∂  ∂ ∂ ∂ (13)

In similar way the boundary - initial conditions are differentiated with respect to shape parameters b_s. So, for surface of the skin one has

0

0

0 1

1

0

D 0,

D D D

: D D D

, D

e s

s s

s

q t t

T b x q

b b x T

t t b

 = ≤

∂ 

 

∈ Γ = λ  ∂  α=  >



(14)

this mean

1

1 0

1

0 1 1

1

1 1 0

, :

,

s s

s

s s

v T

t t

U x x

x W

v T

x U t t

x x

∂ ∂

λ ≤

∂  ∂ ∂

∈ Γ = λ ∂ = α − λ ∂ ∂ >

 ∂ ∂



(15)

or

0 0

0 1

1 1 0

, :

,

s

s

s s

q v t t

x W x

U q v t t

x

 ∂

 ∂ ≤

∈ Γ = 

α − ∂ >

 ∂



(16)

For boundary limiting the system:

3

3: 3 3 U ^s 0

x W

x

∂ 

∈ Γ = −λ  ∂ = (17)

Differentiating the continuity conditions (2) one obtains (c.f. formula (9))

1, 1

1,

: 1, 2

s s

es e e s e

e

es e s

v v

W q W q

x L x x e

U U

+ +

+

∂ ∂

 − = −

=  ∂ ∂ =

 =



(18)

where: q_e = −λ_e∂T_e/∂x, W_es = −λ_e∂U_es/∂x for x = L_e, e = 1, 2 and q_e = λ_e∂T_e/∂x, Wes = λe∂Ues/∂x for x = Le-1, e = 2, 3.

Consistent with the formula (8) the initial condition takes a form

ep

es s

T

U v

x

∂

=

∂ (19)

In the case of sensitivity analysis with respect to shape parameter b¹ = L⁰ we assume the following form of velocity

1

0 1

1 1

1 1

1 3

( , ) ,

0,

L x

L x L

L b

v x b

L x L

 − ≤ ≤

 −

=  ≤ ≤

(20)

while for the problem concerning the sensitivity analysis with respect to shape parameter b² = L¹

0

0 1

2 0

2

2 2 1 2

2 2

2 3

,

( , ) ,

0, x L

L x L

b L

L x

v x b L x L

L b

L x L

 − ≤ ≤

 −

 −

= ≤ ≤

 −

 ≤ ≤



(21)

Taking into account the forms (6), (7) of functionals Ib, Id, the sensitivity of these integrals with respect to the parameters b_s is calculated using the formulas

2 0

D exp d

D

r

r rs

s r r

I E E

P U t

b RT RT

τ ∆  ∆ 

= − 

 

∫

⁽²²⁾

where r = p or r = s and Tb = T¹(L¹, t) = T²(L¹, t), Td = T²(L², t) = T³(L², t), Ubs =

= U1s(L1, t) = U2s(L1, t), Uds = U2s(L2, t) = U3s(L2, t) (c.f. equations (2)).

The change of burn integrals connected with the change of parameters b_s results from the Taylor formula limited to the first-order sensitivity, this means

( ) ( ) D

D

r

r s s r s s

s

I

I b b I b b

± ∆ = ± b ∆ (23)

3. Boundary element method

The primary and also the additional problems resulting from the sensitivity analysis have been solved using the 1^st scheme of the BEM for 1D transient heat diffusion [7]. The boundary integral equations (for successive layers of skin - e = 1, 2, 3) corresponding to the primary problem and the transition t ^f−1 → t ^f are of the form [3, 4, 7]

1

1

1 1

1

1

1 1

1

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

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

1 ( , ) ( , , , ) d d

e f

f

e

e f

e

f

e e

e

e

t x L

f f

e e e

e t x L

x L L

t

f f f f

e e e e

et x L L

L

f f

e B e e me e

e L

T t T x t t q x t t

c

q x t t T x t t T x t t T x t x

c

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

c

−

−

−

−

−

−

=

∗

=

=

∗ ∗ − −

=

− ∗

 

ξ + ξ  =

 

 

 

= ξ  + ξ +

 

 

 

+  − +  ξ

∫

∫ ∫

∫

1 f

f

t

t

x

−

∫

(24)

where T_e^* are the fundamental solutions given by formulas

1 ( )2

( , , , ) exp

4 ( )

2 ( )

f

e f f

e e

x

T x t t

a t t

a t t

∗ ξ = π − − − ξ−  (25)

where ξ is the point in which the concentrated heat source is applied and ae = λe/ce. The heat fluxes resulting from the fundamental solutions are equal to

2 3/ 2

( , , , ) ( ) ( )

( , , , ) exp

4 ( )

4 ( )

f

f e e

e e f f

e e

T x t t x x

q x t t

x a t t a t t

∗

∗ ∂ ξ λ − ξ  − ξ 

ξ = −λ ∂ = π −  − − 

(26)

For ξ → L_e−1⁺ and ξ → Le− for each domain considered one obtains the system of equations

1 1

1 1

11 12 11 12

21 22 21 22

( ) ( )

( , ) ( , )

( ) ( )

( , ) ( , )

f f

e e e e

e e e e

e e e e

f f

e e e e

e e e e

e e e e

p L z L

q L t T L t

g g h h

p L z L

q L t T L t

g g h h

− −

− −

   

       

= + +

   

       

   

      (27)

and then the final form of resolving system results from the continuity conditions (2) ad conditions given for x = L0, and x = L3 (equations (3), (4)):

1 1 1

1 0

11 12 12

1 1 1

1

21 22 22

2 2 2 2

1

11 11 12 12

2 2 2 2

2

21 21 22 22

3 3 3

2

11 11 12

3 3 3

3 3

21 21 22

( , )

0 0 0

( , )

0 0 0

( , )

0 0

( , )

0 0

( , )

0 0 0

( , )

0 0 0

f

f b

f b

f d

f d

f

T L t

h h g

T L t

h h g

q L t

h g h g

T L t

h g h g

q L t

h g h

T L t

h g h

 

− − 

 

− −  

 

 

 − − 

 

 

− −  

 

 

 − − 



 

− − 

 

   

1

11 0 1 0 1 0

1

21 0 1 1 1 1

2 1 2 1

2 2 2 2

3 2 3 2

3 3 3 3

( ) ( )

( ) ( )

( ) ( )

( ) ( )

( ) ( )

( ) ( )

g q p L z L

g q p L z L

p L z L

p L z L

p L z L

p L z L

− + + 

− + + 

 

 + 

=  

 + 

 + 

  

  + 

(28)

In similar way one can solve the additional sensitivity problems.

4. Examples of computations

In numerical computations the following values of parameters have been assu- med [1]: λ¹ = 0.235 W/mK, λ² = 0.445 W/mK, λ³ = 0.185 W/mK, c¹ = 4.3068×

×10⁶ J/m³K, c2 = 3.96⋅10⁶ J/m³K, c3 = 2.674⋅10⁶ J/m³K, c_B = 3.9962⋅10⁶ J/m³K, TB = 37°C, G¹ = 0, Ge = 0.00125 (m³blood/ s)/m³tissue for e = 2, 3, Qm1 = 0, Q_me = 245 W/m³ for e = 2, 3 [1]. Pre-exponential factors: P_b = 1.43⋅10⁷² 1/s for Tb ≥ 317 K and Pb = 0 for Tb < 317 K, while Pd = 2.86⋅10⁶⁹ 1/s for Td ≥ 317 K and Pd = 0 for Td < 317 K. The ratio of activation energy to universal gas constant:

∆E/R = 55000 K. The thicknesses of sub-domains: epidermis L¹ − L⁰ = 0.0001 m, dermis L2 − L1 = 0.002 m and sub-cutaneous region L3 − L2 = 0.01 m. Heat transfer coefficient α = 8 W/m²K, and the ambient temperature T^∞ = 20°C. Time step:

∆t = 0.05 s. The sensitivity of burns integrals has been calculated under the assumption that ∆L⁰ = 10⁻⁵ m and ∆L¹ = 10⁻⁵ m (c.f. equation (23)).

In the first example of computations the heat flux q0 = 6500 W/m² on the skin sur- face has been assumed, the exposure time: t0 = 18 s. The successive skin layers have been divided into 10, 40 and 120 internal cells.

In Figure 2 the temperature distribution in the skin domain is shown. In Figures 3 and 4 the course of burn integral Ib and its courses found on the basis of sensitivity analysis with respect to parameters L0 and L1 are shown (c.f. equations (6), (23)).

The times to the first and second degree burns predicted for the basic value of epi- dermis thickness are equal 16.25 s and 17.8 s, respectively. It is visible, that the thickness of epidermis has essential influence on the times of burns appearance.

Fig. 2. Temperature distribution (q₀ = 6500 W/m², t₀ = 18 s)

Fig. 3. Course of burn integral I_b - change of L₀

Fig. 4. Course of burn integral I_b - change of L₁

In the second example of computations the heat flux q0 = 80000 W/m² on the skin surface has been assumed, the exposure time: t0 = 5 s. The successive skin layers have been divided into 5, 20 and 60 internal cells.

In Figure 5 the distribution of temperature in the skin domain is shown. In Figures 6 and 7 the course of burn integral Id and its courses found on the basis of sensi- tivity analysis with respect to the parameters L0 and L1 are shown (c.f. equations (7), (23)). The time to the third degree burn predicted for the basic values of epi- dermis and dermis thickness is equal 14.75 s. As previously, it is visible, that the thickness of epidermis has essential influence on the time of burn appearance.

Fig. 5. Temperature distribution (q₀ = 80000 W/m², t₀ = 5 s)

Fig. 6. Course of burn integral I_d - change of L₀

Fig. 7. Course of burn integral I_b - change of L₁

References

[1] Torvi D.A., Dale J.D., A Finite element model of skin subjected to a flash fire, Journal of Me- chanical Engineering 1994, 116, 250-255.

[2] Henriques F.C., Studies of thermal injures V. The predictability and the significance of thermally inducted rate processes leading to irreversible epidermal injury, Archives of Pathology 1947, 43, 489-502.

[3] Majchrzak E., Jasiński M., Sensitivity study of burn predictions to variations in thermophysical properties of skin, Advances in Boundary Element Techniques II, Hoggar, Geneva 2001, 273- -280.

[4] Majchrzak E., Jasiński M., Kałuża G., Sensitivity analysis of burn integrals with respect to the geometrical parameters of skin, Computer Methods in Mechanics CMM 2003, Gliwice/Wisła 2003.

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

[6] Dems K., Rousselet B., Sensitivity analysis for transient heat conduction in a solid body - Part II:

Interface modification, Structural Optimization 1999, 17, 46-54.

[7] Majchrzak E., Boundary element method in heat transfer, Wydawnictwo Politechniki Często- chowskiej, Częstochowa 2001.