Galerkin's finite element-laplace transform technique to transient heat migration in human skin and subcutaneous tissues

Bulletin of Mathematical Biology Vol.48, No. 2, pp. 125-136,1986.

Printedin GreatBritain.

0092-8240/86S3.00+ 0.00 PergamonPressLtd. © 1986Societyfor MathematicalBiology

GALERKIN'S FINITE ELEMENT-LAPLACE TRANSFORM TECHNIQUE TO TRANSIENT HEAT M I G R A T I O N IN H U M A N SKIN A N D SUBCUTANEOUS TISSUES •

S. PAL and D. S. PAL Department of Mathematics and Statistics, Punjab Agricultural University, L u d h i a n a - - 141004, India

Galerkin's finite element-Laplace transform technique (GAFELTTE) has been used to study transient temperature distribution i n human skin and subcutaneous tissues. This study incorporates heat conduction, heat carried by perfusion of blood in the capillary beds and metabolic heat generation in the tissues. Different values of various quantities have been considered in all three parts, namely epidermis, dermis and subcutaneous tissues, depending on physiological considerations. The GAFELTTE provides interface temperatures for a wide range of the values of skin surface temperatures. These values have been used to obtain temperature profiles in the region considered. Steadystate temperature distribution has been deduced from the solution obtained by GAFELTTE and has been compared with the results obtained by using different methods.

1. Introduction. The complexity in the structure of skin and subcutaneous tissues (SST) prevents the use of well-known mathematical and ordinary numerical methods in obtaining a useful solution for the study of heat migration in SST. However, the Galerkin's finite element-Laplace transform technique (GAFELTTE) proves to be very useful and provides a solution which incorporates almost all the physical considerations. The core temperature of a human body remains constant in the range 13-57°C of environmental temperature. The SST plays an important role in temperature control of the human body. The skin consists of two layers: epidermis and dermis. The epidermis is composed of living stratum Malaphighii which rests upon the dermis and a dead, horny superficial stratum corneum. The outer dead layer of the epidermis, the stratum corneum, is composed of flattened cells. The dermis, which is the deeper and thicker of the two major layers of the skin, is composed of matted masses of connective tissues and elastic fibres through which pass numerous blood capillaries, lymphatics and nerves. The subcutaneous tissues contain fat cells, connective tissue, larger blood vessels, lymphatics and nerves. In principle, we can say that the physical quantities like the rate of blood mass flow (BMF), rate of metabolic heat generation (MHG) and thermal conductivity of the tissue are different in different layers of SST. 125

126

S. PAL AND D. S. PAL

In this paper we study transient temperature distribution in SST taking into account heat conduction, heat carried by perfusion of blood in the capillary beds and metabolic heat generation in the tissues. We assume that the outer surface of the skin is exposed to a dry and cool environment with negligible perspiration. The physical quantities are assumed to have different values in all the three regions depending on the physiology of SST. Our chief interest in the present study is to obtain interface temperature in the SST for the transient case for various values o f the skin surface temperatures. If T denotes the temperature at time t and position x measured perpendicularly into the tissue from the skin surface, then the rate of change of temperature is governed by

~--~-(K~T~+ m c b ( T b - - T) + S = p c t ~

~x \ ax /

(1)

where K = thermal conductivity of the tissues; m = rate of BMF; cb = specific heat of the blood; Tb = temperature of the blood when it enters the tissue element; S = rate of MHG; p = density of tissue and ct = specific heat of tissue. Chao et al. (1973) and Chao and Yang (1975) have used equation (1) to study temperature variation in SST by considering all parameters as constant with different values o f thermal conductivity for the two layers o f SST. An unsteady-state simplified problem was solved by Stekette and Vander Hock (1979) taking SST as a single region. Saxena and Arya (1981) have discussed steady-state heat distribution by taking three regions of SST. Here, we investigate the transient temperature distribution in SST by using GAFELTTE. The temperature distribution obtained by Saxena and Arya has been deduced as a particular case.

2. Galerkin's Finite Element-Laplace Transform Technique. We take the thickness of epidermis, dermis and subcutaneous regions as a, b -- a and c -- b, respectively. In view o f the properties o f SST, it seems quite reasonable to assume that m and S are negligible in the epidermis and variable in the dermis and subcutaneous tissues. We consider the following linear variations for BMF and MHG. ruder

= moll(x);

msubcu = mofz(X);

Sder

= Sofl(X)

(2)

S~ubcu = Sofz(x)

where fl(x) = (x -- a)/(b -- a) and f2(x) = [ 1 + (x -- b)/(c -- b ) ] . m0 and So are constants. As we know that Tb represents the temperature of the blood entering the region, we can reasonably assume it to be equal to the b o d y

FINITE ELEMENT-LAPLACE TRANSFORM TECHNIQUE

127

core temperature (37 °C). The initial and boundary conditions are considered as T(x, O) = Ta -- (Ta Tb )(X/C) (3) T(O, t) = Ta; T(c, t) = Tb -

-

where Ta is the outer skin surface temperature (Saxena, 1978). Let the trial solution of equation (1) be 4

T(x, t) = ~,Ti(t)4)i(x)

(4)

i=1

where ~i(x) (i = 1, 2, 3, 4) are the basis functions associated with the nodes xi. The nodes xl ,x2, x3 and x4 are associated with x = 0, x = a, x = b and x = c, respectively. The basis function satisfies the essential boundary conditions of the original partial differential equation. Tl(t), T2(t), T3(t) and T4(t) are the temperatures associated with the nodes x l, x2, x3 and x4, respectively. Defining 02T aT L ( T ) = K ~ x 2 + mCb(Tb -- T) + S -- OCt O---t"

(5)

The approximating function T from equation (4) when substituted in L ( T ) will result in a residual. By applying Galerkin's finite element technique, this residual is made small with the help o f orthogonality condition with respect to the basis functions (Connor and Brebbia, 1977),

( L ( T ) , ~b~(x)) = 0.

(6)

Four conditions o f orthogonality can be satisfied as four basis functions have been selected. These are

(t)(~/(x),~bi(x

= O,

where

(7)

i = 1, 2, 3, 4.

Solution o f the above equations results in the values of T2(t) and T3(t), since Tl(t) and T4(t) are known (boundary conditions). As the region of interest is divided into elements defined by nodes, the values of these undetermined coefficients are the values of the dependent variables at the nodes. Substituting equation (5) into (7), we get d ~/ 4 Tj(t)--~ + M Tb - - ~ , Tj(t)~j(x) \

j=l

j=l

where M = me b and a = pct or

4~dT ' + S -- a~A--~, q51(x), c&(x j = l u~

= 0

128

S. PAL AND D. S. PAL 4

l=a

\

dx

-~(t)

~j(x),~(x + {S, ~dx))

4 dr1

- ~ - ~ - <~j(x), ~dx)) = o.

(8)

Integrating the first term of equation (8) by parts, we get

ZTj(t)K

j=l

x=o

i=i

\

d 'j, dx

~/

4

-- ~ Tl(t)(M¢i(x), ¢i(x) ) + (S, ~dx) ) j=l

--~-; a , (otOl(x), ¢,(x)) = O.

(9)

]=1 u t

Since @ is always constant for x = 0 and x = c, therefore d@/dx will always be zero at these nodes. Now writing the above equation in the matrix form, we get [[A] + [B]] [T(t)] + [ C ] [ d T / d t ] = [D] + [E]

(10)

where [A ] 4x4 =

(a~i) = ~K dx'

(11)

[ g ] 4X4 ~- (bii) ~- (M¢i(x), ~)i(x) )

(12)

[ C ] 4 x 4 = ( c i ] ) -~- (og~](x), el(X))

(13)

FTI(t)~

= [ T2(t) ]

(14)

LT,(t)J [ D ] 4 x 1 = (MTb, ¢i(x)) [ g ] 4 × 1 = (S, ¢i(x)).

After finding the values of the above matrices, we get

(15) (16)

FINITE ELEMENT-LAPLACETRANSFORM TECHNIQUE

11 A12 0 21 A22 A23 A32 A33 0 A4a

Ii

TI T2

A3,

ql +

G1

~

C22

/
C2a

c~: c~3 c~4

A,,A

o

=Tb

~

C12 0

2

+

129

[ d~ldt I

c,~ c~[ La~lat_j

2

(17)

LD41] LE~,_I where A l l = K1/a;

A12 =A21 =--K1/a;

A22 = [ 12(/1K1 + aK2) + mol~a]/1211a; A23 =A32 = (mol~ -- 12K2)/1211;

A33 = [12(/2k2 + laK3) + mo(3/~ + 512)]/121~12; A a4 = A4a = (m o l~ -- 4K3)/412;

Cll = C21 = C12 = ala/3;

C23 =

C32 =/10L2/6;

Ca4 i C43 ~_. /2 a3/6; = mo(ll + 512)/6;

A 44 = (12K3 + 7m ol~)/12/2 ;

C22 = (a0tl +/lOt2)/3;

C33 = (/1012 + 120t3)/3;

C44 = 12aa/3; D21 = moll/6; /)41 = 5mo12/6;

E31 = So(/1 + 5/2)/6; E41=5So12/6;

g21 = So/a/6;

ll=b--a,

12=c--b.

K1, K2 and/£3 are the thermal conductivities of the epidermis, dermis and subcutaneous tissues, respectively. After applying the boundary conditions T(0, t) = T1 = Ta and T(c, t) = T4 = 7~, we get

r1oo o]i ] r,oo trol A21 A22 A2a

[~

0

A320 A3a 0 Aa4 1

T2 "t- C21 C22 C23

T3 T4

t00

C32 0 Ca3 0

~

dT2/dt

~T3/dtJ

130 S.PALANDD. S.PAL

(18)

LITb DalTb+ E3,J Further rearranging equation (18) to symmetric form, we get

;1

L23 IT2(t)l __ EL4 Ls] 2 L3

kT3(t)J

s

L6

[dT2/dt] = [L~] kdTa/dtJ Ls

(19)

where L1 =A22;

L2 =/123;

L3 =A3a;

L4 ~-----C22; Ls = --C23;

L6 = --C33;

L7 =TbD21 + E21 --A21Ta ; Ls = TbDal + E31 --Aa4Tb. Equation (19) gives the following set of equations: dT2 dt

L1T2 + L2T3 -- L 4 - - - -

Ls

dT3 - L7 dt

dT2

dT3

dt

dt

L2T2 + LaTa --Ls - - - - L 6

(20)

-- Ls.

(21)

Applying Laplace transform to equations (20) and (21), we get L1T2 + L2T3 --L4 [PT2 -- 7'2(0)1 --Ls [pTa -- T3(0)] -

L2 T2 +

L7

(22)

P

L3 T3 --Ls [P]~2-- T2(0)] -- L6 [pTa -- Ta(0)]

= L__~.s (23) , P where a bar denotes the Laplace transform of the unknowns. Using initial condition [equation (3)] in equations (22) and (23) and after solving these simultaneous equations, we obtain T2 (p) =

Fs(p + F2/2F1) [(p + F2/2F1)2--w: ]

FsF9

Fs(2F4 + F2) 1 2F1 [(p + F2/2Fa)2--w 2 ]

1

F~ [(p + F2/2Fx) 2 -- w 2 ] _t_F7

F6 F1 [(p + F2/2F1) 2 -- w 2 ]

1

F1 p[(p + F2/2F1) 2 -- w 2 ]

(24)

FINITE ELEMENT-LAPLACE TRANSFORM TECHNIQUE

Ta(p) = _ & F l o Fa

131

1

[(p + F2/2Fa) 2 ~ w 21

F1Fg(p + F2/2F~)+ ( F g F n - ½F2F9-- F~2) F1 [(p + F2/2Fx) 2 -- w 2 ]

+FIa 1 Fx p[(p +Fz/2F1)=--w 2]

(25)

where Fx =L4L6--L2s; F4 = L a L 4 - -

F2=2L2Ls--LaL4--LaL6;

LzLs ; Fs = LzL6 -- LaLs ;

F6 =

Fa=LxL3--L~; L6L7 -- LsLs ;

F7 = LaL7 -- LzL8 ; 1

/78 = - [ ( c --a)Ta + 37a]; e

1

F9 = -[(e -- b)Ta + 37b1 ; Flo = L1Ls--L2L4;

Fax = L2L s - L x L 6 ;

C

Fa2=L4Ls--LsL7;

Faa=LxL8--L2L7;

w=(x/F~--4FxFa)/2F1.

Taking inverse Laplace Transformation of equations (24) and (25) we get,

T2(t) = e x p \ F7

T3(t) = exp

2Fa / Fscosh wt + F1---~[FsFg--½Fs(F2 + 2F4) -- F6lsinh wt [. exp(wt)

)-~-~-1t

exp(--wt) l/+

Fgcoshwt

+ 2F~w I(w - - & / 2 & )

[ F22

2~

+ &w(FgFn--½F~Fg--FsF~o--&2)

+ (we+ N / 2 & lI+ F ~ / \ 4 F ~

w

sinhwt

(27)

where Fx > 0" F 2 > 0 . Deduction. When t is very large (t -~ ~), equations (26) and (27) give us the temperatures for the steady-state case as follows: = F , f F 2 -- ~-1 7'2 Fa \4F2 w (28)

132

S. PAL AND D. S. PAL

(29)

T3=FTT2,

In the steady-state case T is independent of time (t) and equation (1) reduces to

3. Direct Formulation for the Steady-state Case. d (K~x)+M(Tb dx

T)+S

O.

(30)

We assume that the trial solution of equation (30) be 4

T(x) = ~ ] ~ ¢ i ( x ) .

(31)

i=l

Applying Galerkin's finite element technique and proceeding in a similar way to the previous case, equation (30) reduces to [[A] + [B]] [T] = [D] + [E]. (32) Use of the boundary conditions reduces equation (32) to ~ 41 ITbl T2 = ITbD2a+E2 I IA51 A22 0 A23 0 1 11.

"" ' o"" o =,

(33)

"

Simplification of equation (33)gives A32 A33J /'3

H2

(34)

where

Ha = TbD21+ E2a --A21Ta;

1t2 = TbD31 +E3a--A34Tb.

Equation (34) gives the followingvalues of T2 and/'3" 1 T2 = -: (A3j/~ - A 2 J / 2 ) ix

1

T2 = -7(A22H2-- A32H1) where A ---=A 2 2 A 3 3 - - A 2 3 A 3 2 .

(35)

(36)

FINITE ELEMENT-LAPLACE TRANSFORM TECHNIQUE

133

4. Numerical Results and Discussion. To solve the system of equations numerically, we make use of following physical constants (Chao et al., 1973; Guyton, 1973). /£1 = 0.5 X 10-3cal/cmsec°C K2 = 0.75 X 10 -3 cal/cm sec °C /£3 = 1.00 X 10-3cal/cmsec°C M = 0.33 X 10 -3 cal/cm 3 sec °C So = 0.88 X 10-Scal/cm3sec P = 1.05 gm/cm 3 • ct = 0.83 cal/g. We can assign any value to the constants a, b and c, depending upon the sample under study. In the present case, we take a = 0.10, b = 0.45 and c = 0.65 cm. The unsteady and steady temperatures computed are shown in Table I. The curves are plotted: (i) temperature vs skin depth and (ii) temperature vs time. The pattern o f variation of temperature (with respect to skin depth) is almost uniform for different skin surface temperatures. We have taken the graph for temperature vs skin depth at Ta = 21 °C for our discussion. There is a steep rise in temperature in epidermis where there is no BMF and MHG (Fig. 1 ). The rate o f rise of temperature is slightly decayed in the T o

=

21oc

37 --

--

2 0 sec

, o , e c

....

6 0 sec

ii

!

8 0 sec 51

Epidermis

el

0

~ / " J/.'~'7"

Dermis

I

0.1

I

0,2

I

L

0.3 0.4 Skin depth (cm)

~/"

c,..^,,+~,~ +~ . . . . Subcutaneous tissue

I

0.5

Figure 1. Temperature vs skin depth at various times.

I

0.6

134

S. PAL A N D D. S. PAL

~,,~.,~oo ,="~,"-~,"'~Cq("q¢"-,I

o r~ 0.)

/

°,..~ g.w

tt~ ©

FINITE ELEMENT-LAPLACE TRANSFORM TECHNIQUE

135

upper parts o f the dermis and subcutaneous tissue owing to different metabolic activities in all the three layers. As the time increases, the temperature curves are getting closer and-closer which indicates the possibility of steady state around 100 sec. T2 and T3 temperatures increase gradually with respect to time (Fig. 2). At a particular time, there is a significant increase in T2 temperature for different values of skin surface temperatures whereas /'3 increases slowly for the same conditions. This indicates that 7"3 temperature, in addition to physiological activities, is also dominated by the body core temperature. 37 _

_

To=IS

----

Ta=15

-----

Ta=I7

....

Ta=I9

- -

To=el

33

29 m

E

25

21

17

~3

I

0

I0

L 20

I

I

I

I

30

40

50

60

L "tO

I

b

80

90

I00

Time

Figure 2. Temperature vs time. The steady-state temperatures computed from equations (28), (2"9)and (35) and (36) are the same. These steady-state temperatures have been compared with those of Saxena and Arya. It is found that T2 and T3 temperatures differ to the extent of 0.08 °C and 0.03 °C, respectively. This difference is due to the consideration o f more suitable assumptions for MHG and BMF and other biophysical data. The useful feature o f the present study has been the prediction of interface temperatures at various times, which are not accessible to a thermograph. The study provides the insight into the way the temperature distribution is affected by MHG and BMF. This investigation may be useful to physiologists and biomedical engineers in formulating the mathematical models for the heat flow problems to detect various disorders in the human body. The number of nodes considered in the present model gives temperatures at the interfaces of epidermis-dermis and dermis-subcutaneous tissues.

136

S. PAL AND D. S. PAL

F o r calculating t h e t e m p e r a t u r e o t h e r t h a n at t h e interfaces, the n u m b e r o f n o d e s has t o be increased, a n d this will a d d s o m e f u r t h e r n u m e r i c a l calculations. G a l e r k i n ' s finite e l e m e n t t e c h n i q u e is well suited to the t r a n s i e n t t e m p e r a t u r e d i s t r i b u t i o n in SST w h e r e t h e region o f i n t e r e s t is divided i n t o more than one part owing to the heterogeneous behaviour of biophysical data. This t e c h n i q u e has t h e a d v a n t a g e o v e r o r d i n a r y n u m e r i c a l t e c h n i q u e s as it can be applied to the s t r u c t u r e o f skin having curved b o u n d a r i e s .

LITERATURE Chao, K. N. and W. J. Yang. 1975. "Response of Skin and Tissue Temperature in Sauna and Steam Baths." Biomechanics Symposium, Vol. 69. New York: American Society of Mechanical Engineers. --, J. G. Eisley and W. J. Yang. 1973. "Heat and Water Migration in Regional Skins and Subcutaneous Tissues." Biomechanics Symposium, Vol. 69. New York: American Society of Mechanical Engineers. Connor, J. J. and C. A. Brebbia. 1977. Finite Element Technique for Fluid Flow, p. 3. London: Butterworth. Guyton, A. C. 1973. Text Book of Medical Physiology, pp. 379, 904. Philadelphia, PA: W. B. Saunders. Saxena, V. P. 1978. "Application of Similarity Transformation to Unsteady State Heat Migration Problems in Human Skin and Subcutaneous Tissues." Proc. Int. Heat Transfer Conference, Vol. 3, p. 65. and D. Arya. 1981. "Steady State Heat Distribution in Epidermis, Dermis and Subdermal Tissues." J. theor. Biol. 89,423. Stekette, J. and M. J. Vander Hock. 1979. "Thermal Recovery of Skin After Cooling." Phys. Med. Biol. 24, 583. R E C E I V E D 7-1 1-83 R E V I S E D 7-30-85