Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds

JBUR-4189; No. of Pages 10 burns xxx (2013) xxx–xxx

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/burns

Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds Shong-Leih Lee *, Yung-Hsiang Lu Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 30013, Taiwan

article info

abstract

Article history:

Heat transfer in a unit three-dimensional skin tissue with an embedded vascular system of

Accepted 15 October 2013

actual histology structure is computed in the present work. The tissue temperature and the

Keywords:

temperature of the tissue over the cross-section of the unit skin area is evaluated. The

Vascular model

resulting one-dimensional function is regarded as the temperature of healthy tissue (or

blood temperatures in artery and vein vessels are solved with a multi-grid system. The mean

Bioheat equation

injured skin but the blood perfusion is still normally working) for large area of skin in view of

Skin burn wound

the symmetric and periodic structure of the paired artery–vein vessels in nature. A three-

Active dynamic thermography

dimensional bioheat equation then is formulated by the superposition of the skin burn wound effect and the healthy skin temperature with and without thermal radiation exposure. When this bioheat equation is employed to simulate ADT process on burn wounds, the decaying factor of the skin surface temperature is found to be a sharply decreasing function of time in the self-cooling stage after a thermal radiation heating. Nevertheless, the boundary of non-healing (needing surgery) and healing regions in a large burn wound can be estimated by tracking the peak of the gradient of decaying factor within 30 s after the thermal radiation is turned off. Experimental studies on the full ADT procedure are needed to justify the assumptions in the present computation. # 2013 Elsevier Ltd and ISBI. All rights reserved.

1.

Introduction

The first attempt to quantitatively describe heat transfer in human tissue with blood flow effect was presented by Pennes [1]. He added a source term vb rb ðc p Þb ð1  kÞðTa  TÞ to the heat conduction equation for tissue temperature by considering heat transfer from blood to tissue, where Ta and T denote arterial blood and tissue temperatures, while rb, (cp)b and vb are density, specific heat and perfusion of blood, respectively. The equilibrium factor k¼

Tv  T Ta  T

(1)

is a prescribed constant in the range 0  k  1 throughout the tissue. Pennes [1] assigned k = 0 by assuming thermal equilibrium between venous blood and tissue (Tv = T). This is the wellknown Pennes equation. The Pennes equation has gained widespread acceptance ever since it was published in 1948, although its validity has been seriously questioned in many applications [2,3]. One of the major problems is the countercurrent heat exchange taking place between artery and vein in paired artery–vein vessels. The net heat lost to the tissue from the vessel pairs was found to be the predominant mode of bioheat transfer for vessels of 50–200 mm in diameter [4–6]. This behavior has been confirmed by the numerical investigation for a branching

* Corresponding author. Tel.: +886 3 5728230; fax: +886 3 5753081. E-mail address: [email protected] (S.-L. Lee). 0305-4179/$36.00 # 2013 Elsevier Ltd and ISBI. All rights reserved. http://dx.doi.org/10.1016/j.burns.2013.10.013 Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

2

burns xxx (2013) xxx–xxx

countercurrent network [7,8]. To take this effect into account, some nonvascular models were proposed to improve the Pennes equation by altering the equilibrium factor such as the efficiency function [3] and the correction coefficient [6]. The artery and vein vessels in skin are typically paired. Their diameters are less than 200 mm in general [6]. Hence, the countercurrent heat exchange in the vessel pairs might be significant. Henriques and Moritz [9] were the pioneers in the study of skin burn. They found that a first degree burn occurs when the skin is maintained at a temperature above 44 8C. Since then, a few mathematical models of heat transfer in skin have been proposed for some thermal hazards including contact burn, scald burn, flash fire, and thermal radiation exposure [10–16]. The grades of skin burn wound are classified according to the severity of the injured skin. The burn wounds of grades I and IIa would heal spontaneously within 3 weeks of the burn, whereas both grades IIb and III need surgery. In clinical practice even an inexperienced doctor has no difficulty in distinguishing grades I and III burn wounds. However, differentiation between grade IIa and IIb is still problematic. Accurate prognosis is only 50–70% in clinical evaluation based on visual inspection [17]. To assess the burn wound grade more accurately, Renkielska et al. [18] developed a noninvasive diagnostic method termed as active dynamic thermography (ADT). They conducted an experiment on young domestic pigs (each weighted approximately 20 kg) in view of the high degree of functional and structural similarity of pig skin to human skin [19]. The wounds were inflicted by an aluminum rod which was applied to the skin at controlled temperature and time. One day after the burn the ADT experiment was performed with a short optical excitation (halogen lamps of 1000 W). The optical heating resulted in a surface temperature rise of about 2.5 8C followed by a self-cooling stage. The skin surface temperature was assumed to decrease exponentially with a constant time constant during the self-cooling stage. The time constant was calculated for each of the burn wounds, while the wound depth was determined by the method of histological analysis [20]. In this experiment, Renkielska et al. [18] found that burn wound having a time constant longer than 10.125 s would heal after 3 weeks of burn spontaneously, otherwise it would be unhealed. It is noted that the thermo-physical properties of the skin tissue do not significantly change 24 h after the burn is inflicted [13,15,21]. Hence, the time constant variation from one skin burn wound to another that Renkielska et al. observed in their experiment [18] can be attributed to the disabled blood perfusion (or destruction of vascular system). Bejan and coworkers [22,23] seems to be the only investigator in the literature that computes heat transfer in 3D triplelayered skin embedded with artery and vein vasculature. However, their mathematical model does not properly reflect the structure of the actual vascular system in skin. In the present work, the heat transfer in a unit 3D skin tissue with an embedded vascular system of actual histology structure is computed. Based on the 3D numerical result, a simple bioheat equation is developed for healthy skin (or injured skin but the blood perfusion is still normally working). The effect of disabled blood perfusion is then taken into account for burned skin by a superposition technique. The bioheat equation then

is employed to simulate the ADT process on burn wounds as in the experiment by Renkielska et al. [18]. A new parameter is proposed to estimate the boundary of non-healing and healing regions in burn wounds of large area from the different responses of the healthy and burned skins.

2.

Three-dimensional vascular model

2.1.

Governing equation

Skin consists of three layers, namely epidermis, dermis, and hypodermis. Fig. 1 shows a histological diagram [24] and a schematic vascular system [6–8] of human skin. The blood vessels (including artery and vein) beneath the muscle are known as primary vessels (see Fig. 1(b)). The blood circulates between the primary vessels and the cutaneous vessels by separate riser vessels. The blood enters the secondary artery in the bottom of the dermis with a temperature Ta0 approximately the same as that in the primary artery. Next, the blood rises to the top of the dermis by the terminal artery, and then flows to the terminal vein through a capillary bed. Subsequently, the blood descends to the secondary vein by the terminal vein. The blood temperature in the secondary vein Tv0 is slightly smaller than Ta0. There is no blood vessel in epidermis. The counterparts of the secondary and terminal vessels illustrated in Fig. 1(b) are observable from the histological diagram [24] of human skin shown in Fig. 1(a). The secondary vessels are the major supplier of blood perfusion for the skin. Destruction of the secondary vessels would lead to the necrosis of the whole skin. Surgery is needed under this situation. It is noted that the secondary vessels are located in the bottom of the dermis layer (see Fig. 1(a)). Thus the differentiation of grades IIa and IIb for a burn wound should depend on the survival of the secondary vessels rather than merely consider the burn depth of the injured skin. The terminal vessels are roughly 20–40 mm in diameter. They form a pair of countercurrent heat exchanger. The spacing between the terminal vessel pairs is the typical length of the capillary bed, 500–1000 mm. The typical diameter of the secondary vessels is 50–100 mm [6,7]. In the present study, diameter and spacing of the terminal vessels are assumed to be 30 and 750 mm, respectively, while the distance between the terminal artery and terminal vein is 30 mm (60 mm center-tocenter). The diameter of the secondary vessels is 75 mm. The average thicknesses of the epidermis, the dermis, and the hypodermis are, respectively, 75, 1500, and 10,000 mm [15]. Fig. 2 shows a simple 3D model for a unit skin area that contains just one single pair of terminal artery and vein. The dimensionless coordinates (x, y, z) is normalized with the thickness of the dermis L = 1500 mm such that the dimensionless thicknesses of the hypodermis and the epidermis are b = 6.667 and d = 0.05, respectively. Due to symmetry, the computational domain b  x  b, 0  y  b, 0  z  1 + d is employed, where b = 0.25. The radius of the secondary vessel is c = 0.025. All of the thermophysical properties in the tissue are assumed constant. After imposing the assumptions and introducing the dimensionless transformation,

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10 burns xxx (2013) xxx–xxx

3

Fig. 1 – (a) The actual histological diagram (reprinted with permission from Pearson Custom Publishing) and (b) a schematic vascular system of human skin.

t¼

ad t ; L2

u¼

T  T1 ; DT

qm ¼

Q m L2 kd DT

(2)

the energy equation inside the tissue of the dermis and the epidermis is expressible as g

      @u @ @u @ @u @ @u ¼ k k k þ þ þ qm @t @x @x @y @y @z @z

(3)

where ad and kd are, respectively, the thermal diffusivity and thermal conductivity of the dermis. Qm is the heat generation by metabolism. The dermis and the epidermis are represented by 0  z  1 and 1  z  1 + d, respectively. The thermal conductivity k is a step function across the two skin layers. Its value is unity in the dermis, and jumps to another constant ked/kd in the epidermis. Similarly, the specific heat g jumps

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

4

burns xxx (2013) xxx–xxx

vessel structure is reflected in the present study by simplifying the secondary artery and vein as a flat surface of width 2c (c  y  c; i.e. diameter of the secondary vessels is 2c) and temperature ub0 ¼ ðua0 þ uv0 Þ=2 on the bottom of the dermis (z = 0), which is the mean temperature of the secondary artery and vein. The boundary condition on the surface z = 0 thus assumes the form @uðx; y; 0Þ ¼ 0 for y > c @z uðx; y; 0Þ ¼ ub0

for 0  y  c

(5)

On the skin surface (z = 1 + d), the convective boundary condition is @uðx; y; 1 þ dÞ þ Biuðx; y; 1 þ dÞ þ ev ¼ 0 @z Bi ¼

Fig. 2 – A 3D vascular model for one single pair of terminal artery and vein.

from unity in the dermis to ðrc p Þed =ðrc p Þd in the epidermis. The reference temperature difference is assigned as DT = Ta0  T1 with T1 being the ambient temperature. Both blood temperatures inside the terminal artery and vein are assumed one-dimensional. Based on the total heat transfer across the wall of the blood vessels, one has  I  dua @u ¼  (4a) ds l @n dz a

l

duv ¼ dz

l¼

 I  @u  ds @n v

m ˙ b ðc p Þb kd L

(4b)

(4c)

where ua and uv denote the blood temperatures inside the terminal artery and vein, respectively. The notation –@u/@n represents the normal heat flux on the tissue side of the vessel surface. The perfusion parameter l is proportional to m ˙ b which is the mass flow rate of blood inside the blood vessels.

2.2.

Boundary conditions

The hypodermis possesses a large thickness and a low thermal conductivity. It provides a good insulation for the muscle under the skin such that the temperature inside the dermis is essentially maintained by the secondary vessels. Hence, the temperature inside the hypodermis is not computed in the present work for simplicity. As observable from the histological diagram in Fig. 1(a), the secondary artery and vein always work in pair in the bottom of the dermis layer. This particular

h1 L ; ked

ev ¼

LEv ked DT

(6a) (6b)

where h1 is the ambient heat transfer coefficient. Ev ¼ m ˙ s h fg is the heat transfer due to sweat evaporation with m ˙ s being the rate of sweat evaporation per square meter and hfg the latent heat of sweat. The boundary conditions in the y-coordinate assume the form @uðx; 0; zÞ ¼ 0; @y

@uðx; b; zÞ ¼0 @y

(7)

due to symmetry, while the periodic boundary condition is imposed in the x-coordinate uðb; y; zÞ ¼ uðb; y; zÞ  uðx; y; zÞ ¼

uðx þ 2b; y; zÞ if x <  b uðx  2b; y; zÞ if x  b

(8)

There are a few 3D vascular models for bioheat equation that take a unit tissue cuboid/cylinder/cone containing just one pair of artery and vein. They include at least unit tissue cuboid in skin, cylinder in muscle and extremity, and cone in kidney. The boundary condition on the lateral surface of the unit tissue was assumed adiabatic in skin [22,23], extremity [7], and kidney [25], while the Dirichlet condition was assumed in muscle [26]. Nevertheless, the two different types of boundary conditions give rise to essentially the same numerical results. The reason will be discussed later. In the present study, both terminal vessels and secondary vessels are formulated in the model such that the geometry of the vascular system inside the unit tissue cuboid is no longer symmetric. Hence, the periodic boundary condition (8) is used instead of the conventional adiabatic boundary condition. The temperatures on the surfaces of the terminal artery and vein are, respectively, uðx; y; zÞ ¼ ua ðzÞ;

uðx; y; zÞ ¼ uv ðzÞ

(9)

The ‘‘initial’’ condition for the artery blood temperature (4a) is ua ð0Þ ¼ ua0

(10)

Owing to the very tiny vessels, the blood temperature inside the capillary bed is essentially the same as the

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

5

burns xxx (2013) xxx–xxx

surrounding tissue at z = 1. Hence, one assigns the ‘‘initial’’ condition uv ð1Þ ¼ u¯ ð1Þ u¯ ðzÞ ¼

1 ad

(11)

ZZ

uðx; y; zÞ dxdy

(12)

ad

for the blood temperature in the terminal vein (4b), where ad is the cross-section area of the dermis.

2.3.

Numerical method and the results

In the present study, a uniform fine grid system is employed in the vicinity of the terminal vessels (0.05  x  0.05 and 0  y  0.05), xi ¼ 0:05 þ ði  1ÞDx;

y j ¼ ð j  1ÞDy

(13)

where Dx = Dy = 0.001, i = 1, 2, 3, . . ., 101, and j = 1, 2, 3, . . ., 51 after a series of grid tests. To conserve computational effort, the other coarse grid system with a grid mesh of 5Dx  5Dy xˆ i ¼ 0:25 þ 5ði  1ÞDx;

yˆ j ¼ 5ð j  1ÞDy

(14)

where i = 1, 2, 3, . . ., 101, and j = 1, 2, 3, . . ., 51, is used for the remainder of the physical domain (0.25  x  0.25 and 0  y  0.25). The temperature of the outer boundary of the fine grid system at x = 0.05 and y = 0.05 is provided by the solution of the coarse grid system. Similarly, the temperature of the inner boundary of the coarse grid system at x = 0.045 and y = 0.045 is determined by the solution from the fine grid system. Both gird systems have the same gird size in z-coordinate, zk ¼ ðk  1ÞDz

(15)

where Dz = 0.01105, k = 1, 2, 3, . . ., 96. The interface of the dermis and the epidermis (z = 1) is located at the middle of the two successive grid points z91 and z92. Through this particular grid arrangement, the step function k is efficiently handled with the integration scheme [27] across the interface to gain a good numerical stability. The thermophysical properties of skin and blood from the literature [11,15] are employed in the present computation, i.e. kd ¼ 0:37 W m1 K1 ; rd ¼ 1200 kg m3 ; ðc p Þd ¼ 3200 J kg1 K1 ;

Fig. 3 – (a) Isotherms on the plane y = 0. (b) Isotherms on the plane z = 0.3.

ked ¼ 0:21 W m1 K1 ; red ¼ 1200 kg m3 ; ðc p Þed ¼ 3580 J kg1 K1 ;

planes y = 0 and z = 0.3 are shown in Fig. 3(a) and (b), respectively. In Fig. 3(a) the increments of the isotherms employed in dermis and epidermis are 0.002 and 0.001, respectively. From Fig. 3(a), one sees that the blood temperature in the terminal artery decreases from ua(0) = 1 to ua(1) = 0.959 while that in the terminal vein increases from

khd ¼ 0:16 W m1 K1 ; rb ¼ 1100 kg m3 ðc p Þb ¼ 3300 J kg1 K1

The typical sweat evaporation rate under normal condition is Ev ¼ 10 W m2 , while the metabolism heat is negligible in skin (Qm = 0) [15]. Under normal condition the blood perfusion rate is vb = 0.024 m3 s1 m3 [15]. It is equivalent to m ˙b ¼ 2:1  108 kg s1 or l = 0.125 in the present configuration. The other two dimensionless parameters in Eq. (6b), the sweat evaporation and the Biot number, would be ev = 0.006803 and Bi = 0.05, respectively, if the reference temperature difference is DT = 10.5 8C and the ambient heat transfer coefficient is h1 = 7 W m2 K1 [15]. Based on these parameters, the steady-state solution of the problem is solved with the conventional central difference scheme and the SIS solver [28]. The resulting isotherms on the

(16)

uv ð1Þ ¼ 0:9583 to uv ð0Þ ¼ 0:994. This reveals a strong countercurrent heat exchange in the paired artery–vein vessels through the gap between them as expected (see Fig. 3(b)). From Fig. 3(b), one sees also that the tissue temperature is quite uniform (approximately 0.9800) in the region of jxj > 0:1 and jyj > 0:1. This finding indicates that use of the periodic boundary condition (8) agrees with the adiabatic boundary conditions commonly adopted by the previous investigators [7,22,23,25] on the lateral surface of the unit tissue. Fig. 4 shows the mean temperature u¯ ðzÞ over the cross-section of the dermis (12) at two representative perfusion parameters l = 0.1 and l = 0.2. The perfusion parameter after exercise could be twice as large as in rest. Surprisingly, the mean temperature of

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

6

burns xxx (2013) xxx–xxx

u¯ ðzÞ ¼ cðzÞ

(17)

where

  q cðzÞ ¼ cð1Þ þ ked ð1  zÞ for 0  z  1 ked   q cðzÞ ¼ cð1Þ þ ð1  zÞ for 1  z  1 þ d ked  cð1Þ ¼

dþ

1 Bi



q ked

 

ev þ f1 Bi

Bi f 2 þ ev q ¼ ked 1 þ Bið0:8ked þ dÞ f 1 ¼ 0:01812  0:03995 Bi þ 0:08439 Bi2  0:08799 Bi3 þ 0:03359 Bi4 f 2 ¼ 0:9797  0:2651 Bi þ 0:1411 Bi2  0:03627Bi3 (18)

Fig. 4 – The mean tissue temperature u¯ ðzÞ and blood temperatures in artery ua(z) and in vein uv(z) at l = 0.1 and l = 0.2.

the dermis u¯ ðzÞ is essentially independent of the perfusion rate. Fig. 5 shows the effect of the Biot number on the mean temperature of the tissue u¯ ðzÞ under normal condition (l = 0.125). The Biot number is seen to have a significant influence on the tissue temperature. Moreover, the mean temperature of the tissue u¯ ðzÞ is always a linear function of z in both dermis and epidermis, while the epidermis has a smaller thermal conductivity (ked/kd = 0.5676) and thus a larger temperature gradient. This implies that the net heat transfer normal to the lateral surface of the unit skin tissue is zero. For convenience, the mean temperature of the tissue u¯ ðzÞ is correlated by

The correlation is valid for 0 < Bi < 1 and ev ¼ 0:006803. The particular function c(z) is the solution of the one-dimensional problem   d dc k ¼0 dz dz dcð0Þ ¼ q; dz

dcð1 þ dÞ þ Bi cð1 þ dÞ þ ðev  Bi f 1 Þ ¼ 0 dz

(19)

It is interesting to note that the present 3D vascular model reveals an effect on the one-dimensional heat conduction problem (19) through the additional term Bi f1. The value of Bi f1 increases from 0.001492 to 0.008560 as the Biot number increases from 0.1 to unity. It has a magnitude of the same order as the sweat evaporation rate. The blood perfusion in dermis seems to prevent heat loss to the environment when the Biot number becomes large.

3. Burned skin with thermal radiation exposure As mentioned earlier, the differentiation of non-healing and healing regions should depend on the survival of the secondary vessels due to the fact that the secondary vessels are the major supplier of blood perfusion for the skin. Destruction of the secondary vessels would lead to the necrosis of the whole skin. In this section, the secondary vessels are assumed to have been destructed in a round area r  rb inside a large injured skin. The heat transfer in the annulus region (rb + 1)  r  r0 will be treated as that in healthy skin where the secondary vessels are still working (healing area). Between them there is a buffer zone rb < r < (rb + 1) where the blood perfusion is partially disabled. For convenience, the round area r  rb is referred to as ‘‘surgery region’’ because surgery is needed there (non-healing region). The bioheat equation for the injured skin tissue covering the surgery region can be written as     @u¯ k @u¯ @ @u¯ @ @u¯ þ k k þ (20) g ¼ @t r @r @r @r @z @z

Fig. 5 – Influence of Biot number on the mean tissue temperature u¯ ðzÞ at l = 0.125.

where (r, z) is a dimensionless cylindrical coordinate system normalized with reference length L. The definition of the mean

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

7

burns xxx (2013) xxx–xxx

temperature of the skin tissue u¯ ðr; z; tÞ is the same as that in the previous section. The formulation is performed in the region of r  r0 and 0  z  1 + d that might cover thousands of terminal vessel pairs. The boundary condition for the region (rb + 1)  r  r0 is 

@u¯ ðr; 0; tÞ ¼q @z

(21)

by noting that the heat flux q is maintained by the secondary vessels (see Eq. (18)). In contrast, the burned secondary vessels provide no heat energy inside the surgery region (r  rb). The heat flux at z = 0 would come from the core temperature uc (at z =  b) beneath the hypodermis. Hence, the boundary condition is assumed 

@u¯ ðr; 0; tÞ ¼ Vðuc  u¯ ðr; 0; tÞÞ; @z

V¼

khd bkd

(22)

Boundary condition (22) can be regarded as an insulation surface because the pseudo-Biot number is very small (V  0.065) in general. In the buffer zone, the boundary condition is assumed to be a linear variation between the two boundary conditions (21) and (22). On the skin surface (z = 1 + d), the convective boundary condition is employed du¯ þ Biu¯ þ ðev  Bi f 1 Þ  z ¼ 0 dz

(23a)

Fig. 6 – Isotherms of u¯ ðr; z; tÞ for the case of rb = 5, Bi = 0.05, ua = 1.02, and z = 0.

by substituting the superposition u¯ ðr; z; tÞ ¼ fðr; z; tÞ þ cðzÞ

LQ tr z¼ ked DT

(23b)

where Qtr is the net intensity of a thermal radiation absorbed by the skin surface. It comes from some heat source such as halogen lamps. It is noted that Eq. (23a) has been corrected by adding the term Bi f1 to match Eq. (19) when the mean temperature u¯ is used instead of the local temperature u for the tissue. Finally, the boundary conditions at r = 0 and r = r0 are assigned as @u¯ ð0; z; tÞ ¼ 0; @r

@u¯ ðr0 ; z; tÞ ¼0 @r

(24)

The axially-symmetrical transient problem (20)–(24) should reduce to the one-dimensional steady problem (19) in the absence of surgery region (and buffer zone) and thermal radiation exposure (z = 0). Therefore, one gets the governing equation     @f k @f @ @f @ @f ¼ þ k k þ (25) g @t r @r @r @r @z @z and the associated boundary conditions @fð0; z; tÞ @fðr0 ; z; tÞ ¼ 0; ¼0 @r @r @fðr; 1 þ d; tÞ þ Bi fðr; 1 þ d; tÞ  zðtÞ ¼ 0 @z @fðr; 0; tÞ  Vfðr; 0; zÞ ¼ q  Vðuc  cð0ÞÞ for 0  r  rb @z @fðr; 0; tÞ  sVfðr; 0; zÞ ¼ sq  sVðuc  cð0ÞÞ @z @fðr; 0; tÞ ¼ 0 for ðrb þ 1Þ  r  r0 @z

for rb < r < ðrb þ 1Þ (26)

(27)

into Eqs. (20)–(24), where s = rb  r + 1. The unsteady term on the left-hand-side of Eq. (25) arises from the time-dependent thermal radiation z(t). It appears that the solution of Eqs. (25) and (26) is f(r, z, t) = 0 for healthy skin without thermal radiation exposure (rb = 0 and z = 0). Fig. 6 shows the resulting isotherms of the solution u¯ ðr; z; tÞ for the case of rb = 5, Bi = 0.05, uc = 1.02, and z = 0, while the computation domain is truncated at r0 = 50. The surgery region has a low temperature as compared to the healthy skin due to the burned secondary vessels. The minimum temperature occurring at the center point of the surgery area is u¯ ð0; 1 þ dÞ ¼ 0:76 on the skin surface. It corresponds to 34 8C if the reference temperature difference and the ambient temperature are DT = 10.5 8C and T1 = 26 8C, respectively. Next, let the solution revealed in Fig. 6 be the initial condition. As suggested by Renkielska et al. [18] a thermal radiation of intensity Qtr = 441 W m2 is applied on the skin surface (covering the burn wound) for 30.42 s to raise the skin surface temperature by about 2.5 8C, then turn off the halogen lamps and leave the skin cooling-off gradually. The thermal radiation in dimensionless form is  0:3 1:303  t  0 (28) zðtÞ ¼ 0 t>0 based on the reference time L2/ad = 23.4 s. The time step employed in the computation is Dt = 0.05 (or 1.17 s). The resulting skin surface temperature is shown in Fig. 7. The surgery region (r  rb) is seen to have a larger temperature decaying factor a¼

@u=@t u

(29)

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

8

burns xxx (2013) xxx–xxx

Fig. 7 – Variation of skin surface temperature after thermal radiation exposure.

as compared to the healthy skin (see Fig. 8). This finding is consistent with the experiment [18]. However, the decaying factor seems to be a sharply decreasing function of time. It decreases from 0.476 to 0.035 within 8.19s at r = 0 (see the curves with t = 1.17 s and t = 9.36 s in Fig. 8). Hence, it is not relevant to grade the burn wounds with the value of a (or the time constant employed in [18]). In the present study, the gradient of the decaying factor @a/@r is proposed for the grade assessment of skin burn wounds. Fig. 9(a) and (b) shows the distribution of @a/@r in the periods of 0  t  30.42 s and 30.42 s  t  152.10 s for the case of Bi = 0.05 and rb = 5. From Fig. 9 – Variation of S@a/@r with time after thermal radiation exposure for Bi = 0.05 and rb = 5.

Fig. 8 – Decaying factor of skin surface temperature after thermal radiation exposure.

Fig. 9(a) the gradient of the decaying factor is seen to have a sharp peak at r = rp near the boundary of the surgery region (rp  rb). As time elapses the peak moves left gradually, while the second peak emerges at t = 143 s (see Fig. 9(b)). For convenience, the ratio rp/rb as a function of time is presented in Fig. 10 for various sizes of surgery region. Fig. 10 indicates that rp is stably closed to rb at the beginning of the cooling stage. The surgery region (r  rb) can be estimated by tracking the first peak of the decaying factor gradient of the skin surface temperature (r = rp) within 30 s after the halogen lamps are turned off. Finally, it is mentioned that most thermo-physical properties of the skin tissue depend on the hydration rate, e.g. thermal conductivity, density, specific heat, absorptivity and emissivity, etc. [29,30] that would give rise to significant influence on the tissue temperature [31]. Great water content in full thickness burns of skin might occur within a few minutes post injury [32]. Moreover, the alteration of blood flow

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10 burns xxx (2013) xxx–xxx

9

(b) The blood perfusion in dermis prevents heat loss to the environment when the Biot number becomes large. (c) The decaying factor of the skin surface temperature is a sharply decreasing function of time in the self-cooling stage after a thermal radiation heating. Thus, it is not relevant to grade the burn wounds with the time constant. (d) The boundary of non-healing and healing regions in a large skin burn wound area can be estimated by tracking the peak of the decaying factor gradient of the skin surface temperature within 30 s after a thermal radiation heating.

Conflict of interest statement The authors declare that there are no conflicts of interest.

Acknowledgment Fig. 10 – The ratio rp/rb as a function of time for various sizes of burn wound.

The authors wish to express their appreciation to the National Science Council of Taiwan for the financial support of this work through the contract NSC 100-2221-E-007-085.

could draw edema formation in the adipose tissue if deep layer skin is damaged [33]. Fortunately, the water content in the skin might have recovered at the time (24 h after burn) when the ADT process is applied for evaluating the burn depth. In this connection, the thermo-regulation feedback effect such as vasodilation and vasoconstriction [34] might have returned to the normal condition too. Hence, all of the thermo-physical properties of the skin tissue are assumed constant in the present study for simplicity. Due to the lack of reliable model, the absorption and scattering coefficients of the skin surface are assumed uniform (z is assumed constant in Eq. (23)) over the entire test area of the ADT process when the thermal radiation exposure is applied. Experimental studies on the full ADT procedure are needed to justify the assumptions.

references

4.

Conclusion

In the present study, a one-dimensional function is obtained for the tissue temperature of healthy skin based on a threedimensional solution of a vascular model. Next, a threedimensional bioheat equation is formulated by superposition of the skin burn wound effect and the normal skin temperature. Finally, the proposed bioheat equation is employed to simulate the method of active dynamic thermography for grade assessment of skin burn wounds. Based on the numerical results, the following conclusions are drawn. (a) The effect of blood perfusion in dermis based on the present 3D vascular model is expressible as an additional term in the convective boundary condition if the tissue temperature is formulated with a conventional onedimensional heat conduction problem. This particular term has a magnitude of the same order as the sweat evaporation rate.

[1] Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol 1948;1:93–122. [2] Wissler EH. Pennes’ 1948 paper revisited. J Appl Physiol 1998;85:35–41. [3] Brinck H, Werner J. Efficiency function: improvement of classical bioheat approach. J Appl Physiol 1994;77:1617–22. [4] Weinbaum S, Jiji LM. A new simplified bioheat equation for the effect of blood flow on local average tissue temperature. J Biomech Eng Trans ASME 1985;107:131–9. [5] Baish JW, Ayyaswamy PS, Foster KR. Heat transport mechanism in vascular tissue: a model comparison. J Biomech Eng Trans ASME 1986;108:324–31. [6] Weinbaum S, Xu LX, Zhu L, Ekpene A. A new fundamental bioheat equation for muscle tissue: Part I—Blood perfusion term. J Biomech Eng Trans ASME 1997;119:278–88. [7] Brinck H, Werner J. Estimation of the thermal effect of blood flow in a branching countercurrent network using a three-dimensional vascular model. J Biomech Eng Trans ASME 1994;116:324–30. [8] Ji Y, Liu J. Vasculature based model for characterizing the oxygen transport in skin tissue – analogy to the WeinbaumJiji bioheat equation. Heat Mass Transfer 2004;40:627–37. [9] Henriques FC, Moritz AR. Studies of thermal injury: the conduction of heat to and through skin and the temperature therein. A theoretical and experimental investigation. Am J Pathol 1947;23:531–49. [10] Diller KR, Hayes LJ. A finite element model of burn injury in blood-perfused skin. J Biomech Eng Trans ASME 1983;105:300–7. [11] Torvi DA, Dale JD. A finite element model of skin subjected to a flash fire. J Biomech Eng Trans ASME 1994;116:250–5. [12] Ng EYK, Chua LT. Comparison of one- and twodimensional programmes for predicting the state of skin burns. Burns 2002;28:27–34. [13] Autrique L, Lormel C. Numerical design of experiment for sensitivity analysis—application to skin burn injury prediction. IEEE Trans Biomed Eng 2008;55:1279–90.

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013

JBUR-4189; No. of Pages 10

10

burns xxx (2013) xxx–xxx

[14] Hirata A, Fujiwara O. Modeling time variation of blood temperature in a bioheat equation and its application to temperature analysis due to RF exposure. Phys Med Biol 2009;54:N189–96. [15] Ng EYK, Tan HM, Ooi EH. Boundary element method with bioheat equation for skin burn injury. Burns 2009;35:987– 97. [16] Xu F, Wang PF, Lin M, Lu TJ, Ng EYK. Quantifying the underlying mechanism of skin thermal damage: a review. J Mech Med Biol 2010;10:373–400. [17] Heimbach D, Engrave L, Grube B, Marvin J. Burn depth: a review. World J Surg 1992;16:10–5. [18] Renkielska A, Nowakowski A, Kaczmarek M, Ruminski J. Burn depths evaluation based on active dynamic IR thermal imaging—a preliminary study. Burns 2006;32:867– 75. [19] Meyer W, Schwarz R, Neurand K. The skin of domestic mammals as a model for human skin with special reference to domestic pig. Curr Probl Dermatol 1978;7:39– 52. [20] Singer AJ, Berruti L, Thode HC, McClain SA. Standardized burn model using a multiparametric histologic analysis of burn depth. Acad Emerg Med 2000;7:1–6. [21] Museux N, Perez L, Autrique L, Agay D. Skin burns after laser exposure: histological analysis and predictive simulation. Burns 2012;38:658–67. [22] Wang H, Dai W, Bejan A. Optimal temperature distribution in a 3D triple-layered skin structure embedded with artery and vein vasculature and induced by electromagnetic radiation. Int J Heat Mass Transfer 2007;50:1843–54. [23] Dai W, Wang H, Jordan PM, Mickens RE, Bejan A. A mathematical model for skin burn injury induced by radiation heating. Int J Heat Mass Transfer 2008;51:5497–510.

[24] Marieb EN. Human anatomy and physiology laboratory manual, 7th ed., San Francisco, CA: Pearson Custom Publishing; 2006: 68. [25] Chen C, Xu LX. A vascular model for heat transfer in an isolated pig kidney during water bath heating. J Heat Transfer Trans ASME 2003;125:936–43. [26] Zhu L, Xu LX, He Q, Weinbaum S. A new fundamental bioheat equation for muscle tissue: Part II—Temperature of SAV vessels. J Biomech Eng Trans ASME 2002;124:121–32. [27] Lee SL, Tzong RY. An enthalpy formulation for phase change problems with a large thermal diffusivity jump across the interface. Int J Heat Mass Transfer 1991;34:1491– 502. [28] Lee SL. A strongly-implicit solver for two-dimensional elliptic differential equations. Numer Heat Transfer B 1989;16:161–78. [29] Takata AN. Development of criterion for skin burns. Aerosp Med 1974;45:634–7. [30] Zhang JZ, Shen YG, Zhang XX. A dynamic photo-thermal model of carbon dioxide laser tissue ablation. Laser Med Sci 2009;24:329–38. [31] Deng ZS, Liu J. Mathematical modeling of temperature mapping over skin surface and its implementation in thermal disease diagnostics. Comput Biol Med 2004;34:495– 521. [32] Leape LL. Early burn wound changes. J Pediatric Surg 1968;3:292–9. [33] Sakurai H, Nozaki M, Traber LD, Hawkins HK, Traber DL. Microvascular changes in large flame burn wound in sheep. Burns 2002;28:3–9. [34] Xu F, Lu TJ, Seffen KA, Ng EYK. Mathematical modeling of skin bioheat transfer. Appl Mech Rev Trans ASME 2009;62:050801.

Please cite this article in press as: Lee S-L, Lu Y-H. Modeling of bioheat equation for skin and a preliminary study on a noninvasive diagnostic method for skin burn wounds. Burns (2013), http://dx.doi.org/10.1016/j.burns.2013.10.013