## Abstract

In this paper, an axisymmetric model of the human skin is developed to simulate the steady-state temperature distribution during contact with a hot solid. Simulations are carried out using the boundary element method. This study seeks to investigate the feasibility of using the boundary element method in the studies of burn. A sensitivity analysis is carried out to examine the effects of various parameters on the temperature distribution inside the skin during burn. Furthermore, a statistical analysis based on the Taguchi method is performed to determine the combination of factors that produce the desired outcome (least increase in temperature). In order to validate the accuracy of the numerical scheme, results obtained using the boundary element method are compared with the solutions obtained using the more established finite-element method.

## 1. Introduction

In the USA, two million cases of burn injuries are reported annually. Of the two million cases reported, about 50 per cent are treated by medical professionals and 25 per cent are hospitalized (Garner & Magee 2005). Of the 25 per cent of patients who are hospitalized, approximately 20 000 of them have burn injuries that involve more than 25 per cent of their total body surface burnt. An estimated 10 000 fatalities as a result of burn injuries are reported while almost one million patients will have substantial or permanent disabilities resulting from their burn injury.

Various factors may contribute to burn injuries, such as heat, cold, electricity, chemicals, light, radiation and friction. The degree of thermal damage inflicted onto the human skin during burn may differ from a minor first-degree wound to a more severe case of death. A serious thermal injury generally affects the local and systemic responses of the human physiology. According to Henriques & Moritz (1947), when the temperature of skin is maintained at 44^{°}C, the rate of injury caused is higher than that of recovery, i.e. at a temperature above 44^{°}C, a first-degree burn injury occurs. The degree of thermal damage sustained by the skin is also dependent on both the heating temperature and the duration of heating. An increase in tissue temperature or a prolonged heating duration may result in second- or third-degree burns.

Mathematical modelling has proved to be a versatile approach in the studies of the temperature distribution sustained by human tissue during burn (Diller *et al.* 1983; Charny 1992). This may be due to the advancement in computational technology, which allows accurate numerical simulation of biological systems. The difficulty in carrying experimental investigations on living human tissues may also be a contributing factor. One of the more commonly used mathematical models for describing the temperature distribution inside biological tissues is the Pennes bioheat equation (Pennes 1948) or the bioheat equation in short. It is practically a simple model that distinguishes itself from the classical heat diffusion equation by including the thermal effects owing to tissue blood perfusion and metabolic heat generation.

In simplified terms, the bioheat equation may be mathematically expressed as
1.1
where *κ* and *T* are the thermal conductivity and temperature of the tissue, respectively, *ρ*_{bl}, *c*_{bl} and *ω*_{bl} are the density, specific heat and perfusion of blood, respectively, *T*_{bl} is the temperature of blood and *Q*_{m} is the metabolic heat generation. The second and the third term on the left-hand side of equation (1.1) are associated with the thermal effects owing to blood perfusion and metabolic heat generation, respectively. Although various improved and alternative models have emerged, the Pennes bioheat model is still widely in use owing to its accuracy and its simplicity in numerical coding. It has since become the starting point of scientific modelling work in bioheat transfer.

Mathematical models that adopt the bioheat equation to investigate the temperature distribution inside skin tissues during burns have been reported in the literature. Among the more notable works include the models developed by Diller *et al.* (1983) and Ng & Chua (2000, 2001, 2002). These models were developed using the finite-difference method and the finite-element method. These numerical methods are generally categorized as *domain methods*, i.e. numerical methods that require a fill domain discretization in order to obtain a numerical solution.

The boundary element method is a numerical technique that is relatively new when compared with the finite-difference method and the finite-element method. Unlike the domain methods, the boundary element method is known as a *boundary method*, i.e. the method requires only the boundary of the solution domain to be discretized. As a result, this numerical approach generally requires lesser computational demand (in terms of computational memory and execution time) than domain methods. Although the boundary element method has not been used in engineering analysis as frequently as the finite-difference and finite-element methods, its role as a leading alternative in numerical analysis has remained unchallenged (Cheng & Cheng 2005).

The boundary element method has been used to solve various problems pertaining to the heat transfer inside biological tissues (Ooi *et al.*2007, 2008). However, to the best of the authors’ knowledge, boundary element models of the human skin that are developed specifically for studying burns owing to contact with a hot solid have not been reported. It may be interesting to examine the feasibility of this method in simulating the temperature distribution inside the human skin during burns.

In the present paper, a model of the human skin is developed to predict the temperature distribution during contact with a hot solid. Simulations are performed numerically using the boundary element method. Ideally, by accurately predicting the temperature changes within the skin tissues during burns using proper computer-programmed codes, we may estimate the seriousness of different burns effectively. Nonetheless, in this preliminary study, we shall only consider the steady-state temperature changes inside the human skin.

## 2. Mathematical model

### (a) Model of the human skin

The human skin is modelled as a multi-layered cylindrical region using the axisymmetric (*O**r**z*) coordinate system. Three distinct regions are considered, namely the epidermis (top layer), dermis (middle layer) and subcutaneous fat (bottom layer), which we denote by *R*_{1}, *R*_{2} and *R*_{3}, respectively. These regions along with their dimensions are shown in figure 1. The skin model is enclosed by six boundaries, namely *C*_{1}, *C*_{2}, *C*_{3}, *C*_{4}, *C*_{5} and *C*_{6}. Note that in the axisymmetrical model, the axis of symmetry, i.e. the *z*-axis, does not form part of the boundary.

To simulate burns caused by contact with a hot solid, a heating disc, which is at a constant temperature, is placed on the surface of the epidermis as illustrated in figure 1.

### (b) Governing equation

As described in the previous section, the bioheat equation may be used to describe the temperature distribution inside the human skin model during burns. However, since blood perfusion exists only inside the layer of dermis (Ng & Chua 2000, 2001, 2002), the term associated with the thermal effects of blood perfusion may be omitted when describing the layers of epidermis and the subcutaneous fat. Likewise, it is assumed that the metabolic heat generated by the skin tissues has negligible thermal effects in the presence of a heat source. Consequently, the equations governing the steady-state temperature distribution inside the skin with reference to the model shown in figure 1 may be expressed as
2.1
and
2.2
where *κ*_{i} and *T*_{i} are the thermal conductivity and temperature of region *R*_{i}, respectively, *ρ*_{bl}, *c*_{bl} and *T*_{bl} are the density, specific heat and temperature of blood, respectively, and *ω*_{bl} is the blood perfusion rate inside region *R*_{i}.

### (c) Boundary conditions

Boundary conditions are specified on boundaries *C*_{1}–*C*_{6} based on the anatomical and physiological observations of actual human skin.

Boundary *C*_{1} represents the bottom surface of the human skin. The temperature here is assumed to be the same as the body core and is not affected by the heating disc placed on the surface of the human skin. Hence, we may write
2.3
where *T*_{b} is the body core temperature.

Boundaries *C*_{2}, *C*_{3} and *C*_{4} are assumed to be thermally insulated (Ng & Chua 2000, 2001, 2002), such that
2.4
where ∂*T*_{i}/∂*n* denotes the rate of change of *T*_{i} in the direction-outward unit vector normal to boundaries *C*_{2}, *C*_{3} and *C*_{4}.

Boundary *C*_{5} is exposed to the environment. For an ambient temperature that is lower than the skin surface temperature, heat loss from the skin to the environment takes place via convection and radiation. Cooling is aided by the evaporation of sweat from the surface of the skin. Assuming that the heat loss owing to radiation is negligible, the boundary condition here may be written as
2.5
where *h*_{amb} is the ambient convection coefficient, *T*_{amb} is ambient temperature and *E*_{vap} is the amount of heat loss owing to the evaporation of sweat.

The heating disc is placed on boundary *C*_{6}. Consequently, we may write
2.6
where *T*_{disc} is the temperature of the heating disc.

At the interfaces between two contiguous regions, the continuity condition applies, i.e.
2.7
where *I*_{ij} denotes the interface between two adjacent regions, *R*_{i} and *R*_{j}.

### (d) Material properties

Values of the parameters used to describe equations (2.1)–(2.6) may be obtained from the literature. These are summarized in table 1.

The value of skin evaporation rate presented in table 1 merits some discussion. According to Hedge (2008), heat loss owing to evaporation contributes to approximately 15 per cent of the total heat loss from the surface of the skin. Therefore, if we assume convection and evaporation to be the only two mechanisms responsible for the cooling of the human skin surface, we may simply deduce that 85 per cent of the heat loss is caused by convection, which we may calculate using
where *q*_{conv} is the heat flux owing to convection and *T*_{s} is the temperature at the surface of the skin exposed to the environment. Assuming that the average temperature at the surface of the skin is 34^{°}C (Haraldson & Kopp 1983), and using values of *h*_{amb} and *T*_{amb} given in table 1, *q*_{conv} is calculated to be 63 W m^{−2}. Using this value, the 15 per cent heat loss owing to evaporation is estimated to be approximately 10 W m^{−2}.

Values presented in table 1 are taken as the control, where they are varied accordingly in the sensitivity analysis carried out in §5.

## 3. Numerical implementation

To solve equations (2.1) and (2.2) subjected to equations (2.3)–(2.7) using the boundary element method, the boundary integral representations of equations (2.1) and (2.2) have to be derived. Following the steps outlined by Brebbia *et al.* (1984), we arrive at
3.1
and
3.2
where Γ_{i} denotes the curve boundary of region *R*_{i}, *d**s*(*r*,*z*) denotes the length of an infinitesimal part of curve Γ_{i} and *d**R*(*r*,*z*) denotes the area of an infinitesimal portion of region *R*_{i}, *λ*(*ξ*,*η*)=1 if (*ξ*,*η*) lies at the interior of *R*_{i}, *λ*(*ξ*,*η*)=1/2 if (*ξ*,*η*) lies on a smooth part of Γ_{i} and Φ(*r*,*z*;*ξ*,*η*) and ∂Φ(*r*,*z*;*ξ*,*η*)/∂*n* are the axisymmetric fundamental solution of the Laplace equation and its normal derivative, respectively, which is given by
3.3
where *n*_{r} and *n*_{z} are the components of the outward unit normal vector on Γ_{i} in the *r* and *z* direction, respectively, *K* and *E* denote the complete elliptic integral of the first and the second kind, respectively, as defined by Abramowtiz & Stegun (1972) and the variables *m*, *a* and *b* are expressed as
Details on the derivation of equations (3.1)–(3.3) may be found in textbooks on the boundary element method such as in Brebbia *et al.* (1984).

For a numerical method based on equations (3.1) and (3.2), we discretize the boundary Γ_{i} into *N*_{i} elements known as boundary elements such that If the starting and ending points (*r*,*z*) of the element (for *k*=1,2,…,*N*_{i}−1,*N*_{i}) are denoted by and , respectively, we choose the midpoint on as the boundary collocation point, such that
3.4

Across each boundary element , the temperature and heat flux (−*κ*∂*T*/∂*n*≡*q*) are assumed to be constant such that
3.5
where and and are temperature and heat flux at , respectively. Note that the approximations in equation (3.5) are known as the constant element approximation (París & Cañas 1997).

Substituting equation (3.5) into equations (3.1) and (3.2) leads to 3.6 and 3.7 respectively, where 3.8 The integrals in equation (3.8) may be evaluated numerically using, for example, the Gaussian quadrature.

The domain integral in equation (3.7) may be approximated into an equivalent boundary integral using the dual reciprocity method (Brebbia & Nardini 2000). To do so, *L*_{i} points denoted by for *i*=2, are selected at the interior of *R*_{i}. These points along with the *N*_{i} nodes selected on each boundary element serve as collocation nodes for the implementation of the dual reciprocity method.

Employing the steps outlined by Brebbia & Nardini (2000), the domain integral in equation (3.7) may be reduced into a boundary integral, giving
3.9
where for *m*=1,2,…,*N*_{i}+*L*_{i} and the coefficients and are defined implicitly by
3.10
and
3.11
respectively. The local interpolating function, , and the particular solution, , may be expressed as
3.12
and
3.13
respectively. Details on the functions in equations (3.12) and (3.13) are given by Wang *et al.* (2003).

Substituting equation (3.9) into (3.7) and letting (*ξ*_{i},*η*_{i}) in equations (3.6) and (3.7) be given by , for *n*=1,2,…, *N*_{i}−1,*N*_{i} (for *i*=1,3) and for *n*=1,2,…,*N*_{i}+*L*_{i}−1,*N*_{i}+*L*_{i} (for *n*=2), one obtains
3.14
and
3.15
where and are values of *λ*_{i}(*ξ*_{i},*η*_{i}) and *T*_{i}(*ξ*_{i},*η*_{i}) at , respectively.

Equations (3.14) and (3.15) along with equations (2.3)–(2.7) result in a system of linear algebraic equations that can be solved for the unknown temperature and heat flux at the boundaries and interior of each region *R*_{i}.

## 4. Results and discussion

### (a) Mesh-independent study

Before carrying out simulations on the temperature distribution inside the human skin during burns, it is important to ensure that the numerical scheme yields results that are mesh-insensitive. To do so, a mesh-independent study is carried out to determine the level of discretization that leads to convergent numerical results. Results from the mesh-independent study (not reported here) suggest that a total of 82 boundary elements and 64 interior points (selected inside the dermis, *R*_{2}) are sufficient to produce numerical solutions that are mesh-independent.

### (b) Temperature distribution using control values

Figure 2 shows the temperature profile inside the skin model obtained using values of parameters tabulated in table 1. The lowest temperature occurs at the bottom of the skin (depicted in black), which is at the body core temperature of 37^{°}C. The hottest surface (depicted in grey) is located at the surface where the heating disc is applied. The solid lines indicate the isotherms inside the human skin. Temperature at the surface of the epidermis is found to be approximately 55^{°}C. This may be due to the exposure of the human skin to the environment.

As the skin comes in contact with the heating disc, heat flows from the heating disc into the skin. The excess heat is transferred to the environment via convection and evaporation from the surface of the epidermis. Blood perfusion inside the dermis also serves as a medium to diffuse away the additional heat from the skin tissue.

Figure 3 plots the temperature distribution along the *z*-axis at *r*=0, 0.003 and *r*=0.0062 m. These values of *r* are selected to provide a more complete visualization of the temperature distribution inside the human skin. The temperature of the skin appears to increase uniformly inside the region occupied by the subcutaneous fat. Inside the dermis, a rapid increase at *r*=0 and *r*=0.003 m is observed. This may be due to the close proximity to the heating disc (*r*<0.005 m). At *r*=0.0062 m, the rate of temperature rise drops, which may be explained by the heat loss from the epidermis to ambient.

In order to validate the numerical results obtained using the boundary element model, comparisons are made with the results obtained using COMSOL Multiphysics (a partial differential equation calculator that utilizes the finite-element method). Table 2 summarizes the comparisons between the results at various points inside the skin model using the boundary element and finite-element models. An average difference of approximately 0.068^{°}C and an average percentage difference of less than 1 per cent indicate that the boundary element solution is accurate.

In the finite-element model with the grid-independent study, a total number of 6615 triangular elements have been used. On the other hand, the boundary element model requires a discretization that comprises only 146 nodes (82 boundary elements and 64 collocation nodes for the dual reciprocity method). It is apparent that the boundary element method has the advantage of reducing the computer memory needed.

Discrepancies between the results obtained using the boundary element method and the finite-element method as shown in table 2 may be due to two reasons. Firstly, in the numerical scheme developed in §3, constant element approximations have been used to describe the variations of temperature and heat flux across each boundary element. On the other hand, quadratic shape functions have been used in the finite-element model. The fact that the boundary element method manages to produce numerical results that are very close to the finite-element method despite using only constant element approximation highlights the advantages of the boundary element method. Secondly, when using the axisymmetric boundary element method, numerical integration has been employed to calculate the integration of the fundamental solution. Errors arising from the numerical integration may contribute to the overall errors inside the numerical scheme.

Investigations in the present study have been limited to a steady-state analysis. In clinical practice, this implies that the heating disc is placed onto the skin surface for an infinitely long time, such that the temperature increase in the skin reaches a steady state. In order to obtain a more reliable study of the temperature distribution inside the skin during burns, and ultimately for an analysis of the degree of thermal damage inflicted, a time-dependent analysis has to be carried out. This may be investigated in future studies.

## 5. Sensitivity analysis

The temperature distribution inside the human skin during contact with a hot solid has been successfully simulated in the preceding section. The results presented were obtained using values of thermal properties at the control level, i.e. values that we assume to represent the conditions of the normal human skin. Changes in the values of these thermal properties may lead to temperature distributions that are different from the one observed in §4.

In this section, a sensitivity analysis is carried out to investigate how changes in the various parameters of the human skin during burns may affect the skin temperature distribution. We also aim to identify the parameter(s) that have dominant effects on the human skin. Since the seriousness of burns and the thermal injury inflicted on the skin are dependent on the temperature increase, it may be advantageous if the clinical significance of each parameter is studied.

Eight parameters are considered here, namely the thermal conductivity of the epidermis, dermis and subcutaneous fat, blood perfusion rate of dermis, disc temperature, ambient convection coefficient, ambient temperature and evaporation rate of the skin surface. To investigate how each of these parameters may affect the skin temperature distribution, values of these parameters are varied by a certain amount from the control level (table 1). For each parameter investigated, values of the other parameters are maintained at its control.

Observations are made at various points along the *z*-axis at *r*=0 m, *r*=0.003 m and *r*=0.0062 m. As in figure 3, these plots are selected in order to obtain a more complete visualization of how the changes in each parameter affect the temperature distribution inside the skin.

### (a) Effects of epidermis thermal conductivity, *κ*_{1}

The control value of the epidermis thermal conductivity is chosen to be 0.21 W m^{−1} K^{−1}. In reality, factors such as individual variation and different parts of the human body may lead to different values of thermal conductivity. Assuming that these factors account for changes as much as ±50 per cent from the control, values of 0.105 and 0.42 W m^{−1} K^{−1}, which correspond to values at 50 and 150 per cent of the control, are investigated.

Figure 4 illustrates the results of this investigation. Variation of the epidermis thermal conductivity appears to have a very small effect on the temperature distribution inside the human skin. This may be due to the very thin structure of the epidermis, which causes any changes in the thermal conductivity to have an insignificant impact on the thermal response of the entire system. The results obtained here suggest that in examining the damage to the skin, the epidermis will most likely suffer the worst injury since the heat from the burn source is transferred directly across this region.

### (b) Effects of dermis thermal conductivity, *κ*_{2}

Using the same argument as that presented for the epidermis thermal conductivity, values of 50 and 150 per cent of the control are investigated for the dermis thermal conductivity. These values correspond to 0.185 and 0.74 W m^{−1} K^{−1}, respectively. The results obtained are presented in figure 5. Variation of the dermis thermal conductivity does not appear to have any significant effect on the temperature distribution inside the skin.

### (c) Effects of subcutaneous fat thermal conductivity, *κ*_{3}

Values of the subcutaneous fat thermal conductivity corresponding to 50 and 150 per cent of the control are given by 0.08 and 0.32 W m^{−1} K^{−1}, respectively. The results obtained are illustrated in figure 6. From figure 6, an increase in the subcutaneous fat thermal conductivity is found to produce a slight decrease in the temperature at the interface between the subcutaneous fat and the dermis layer (*z*=0.010 m). In other words, a larger fat thermal conductivity produces a cooler tissue temperature. A larger thermal conductivity permits the skin tissue to cool down faster, thus resulting in less damage to the tissue.

### (d) Effects of dermis blood perfusion rate, *ω*_{bl}

The control value of the blood perfusion rate inside the dermis is chosen to be 0.024 m^{3} s^{−1} m^{−3}. This value is assumed to represent the value at the high level. For investigation purposes, the mid value is selected to be 50 per cent of the control, i.e. 0.012 m^{3} s^{−1} m^{−3} while zero perfusion is chosen to be the low value. Note that, in reality, the blood perfusion inside the human skin may increase to 10 times its normal value during burn (Bowman 1985). This factor is not considered here since the present investigation is confined to a steady-state analysis.

Figure 7 shows the results obtained for the different values of blood perfusion rate investigated. The blood perfusion rate is found to have a negative effect on the temperature distribution inside the human skin, i.e. a larger blood perfusion rate produces a cooler skin temperature distribution. This agrees well with the thermoregulatory response of the skin tissue, as according to Ong & Ng (2005) and Ng & Lim (2008).

The negative effect of the blood perfusion rate may be explained by the role of blood perfusion. During burn, blood flow inside the dermis acts as a medium that transfers away excessive heat that is dissipated into the tissue from the heat source, in an attempt to cool the skin temperature. Thus, as the blood perfusion rate increases, the overall tissue temperature drops. The results observed in figure 7 suggest that the thermal damage of skin during burn may be reduced by an increased blood flow through the tissue.

### (e) Effects of disc temperature, *T*_{disc}

Values of 60 and 120^{°}C, which correspond to the low and high values, are chosen for the investigations of disc temperature. Results obtained from this analysis may provide us with an idea of how different heating temperatures may affect the thermal response of the skin tissues. The results are presented in figure 8.

Figure 8 reveals that the skin temperature at *r*=0 shows a small thermal variation when compared with the variations at *r*=0.003 and 0.0062 m. This is not surprising given that the source of heat lies in the regions along *r*=0 while at *r*=0.0062 m, cooling takes place at the surface of the skin that is exposed to ambient.

### (f) Effects of ambient convection coefficient, *h*_{amb}

We have chosen 7 W m^{−2} K^{−1} to represent the ambient convection coefficient at the control, which implies a free convection flow. This value is taken as the low level in this sensitivity analysis. Values of 25 and 250 W m^{−2} K^{−1}, which correspond to the transition from free to forced convection and forced convection, respectively, are chosen as the corresponding mid and high values (Incroprera & Dewitt 2002). Results are displayed in figure 9.

Changes in the temperature distribution at *r*=0 are minimal. On the other hand, one may observe significant changes in temperature at *r*=0.0062 m. When forced convection is induced on the skin surface, temperatures inside the skin appear to be greatly reduced. In clinical practice, this may imply that serious thermal damage can be avoided by inducing forced convection at the skin surface.

### (g) Effects of ambient temperature, *T*_{amb}

The control value for ambient temperature is 25^{°}C, which is the temperature of an air-conditioned room in Singapore. To investigate the effects of a lower and higher ambient temperature compared with control, values of 16 and 40^{°}C are considered. The value of 16^{°}C is chosen based on the lowest temperature that can be reached in a standard air-conditioned room while the 40^{°}C value represents an environment of extreme heat.

The results are presented in figure 10. The effects of the various ambient temperatures at both *r*=0 and at *r*=0.0068 m appear to be very small. This implies that reducing the ambient temperature may not help to avoid the risk of thermal damage on the skin of individuals who have suffered burns.

### (h) Effects of skin evaporation rate, *E*_{vap}

The values of evaporation rate represent the heat loss owing to the evaporation of sweat from the human skin. To investigate how lower and higher values may change the temperature distribution of the skin, values of 5 and 15 W m^{−2} are used. Results are given in figure 11. The different values of the evaporation rate show virtually no effect on the temperature distribution inside the human skin. This is only to be expected since heat loss owing to evaporation contributes only approximately 15 per cent of the total heat loss from the skin surface (Hedge 2008). Perhaps a larger evaporation rate may have some noticeable effects.

### (i) Discussion

Results from sensitivity analysis suggest that blood perfusion rate inside the dermis and ambient convection coefficient have a significant effect on reducing the temperature rise inside the skin during burn. It may also be inferred that the effects of these factors depend on the location of the analysis. For example, ambient convection coefficient has a strong effect on the temperature at a location that is further away from the heat source than points that are closer (figure 9). Although the sensitivity analysis has been performed for a steady-state problem, which may not be clinically practical, it is hoped that the information obtained may be useful to researchers as preliminary data for understanding the behaviour of the human skin towards heat injury. The Taguchi method can be carried out to improve the analysis.

## 6. Taguchi analysis

The Taguchi method is a statistical tool developed by Genichi Taguchi that is largely used for improving the quality of manufactured goods (Antony & Antony 2001). The Taguchi method is capable of addressing the number of different conditions needed to reliably ascertain whether low or high values are better in an experiment (Sudharsan & Ng 2000; Ng *et al.* 2008). Although the Taguchi method has been developed for the purpose of improving the quality of manufactured items, it has recently been applied in the field of numerical analysis (Sudharsan & Ng 2000; Ng *et al.* 2008). In this case, the Taguchi method is able to identify the combination of factors that produces the desired outcome (Ng *et al.* 2009).

The eight parameters that have been investigated in §5 are again investigated in this section. These eight factors, which are the thermal conductivity of the epidermis, dermis and subcutaneous fat, blood perfusion rate, disc temperature, ambient convection coefficient, ambient temperature and evaporation rate, are denoted by A, B, C, D, E, F, G and H, respectively. It is decided that a 2^{n} experiment is to be used, where *n* signifies the number of factors considered, eight in this case. Each factor has two levels, namely high and low. This results in a total of 256 (2^{8}) runs. The low and high values of each parameter are taken to be the same as those used in §5. These are summarized in table 3.

An important variable in the Taguchi method is the signal-to-noise ratio (SNR). Three types of SNRs are available, namely nominal the best, larger the better and smaller the better. Since we are interested in the largest temperature achieved inside the skin, it is decided that the larger the better SNR is used. Mathematically, this is given by
6.1
where *y* is the observation made, which is the average temperature across the skin model.

All the SNR values calculated using equation (6.1) are consolidated based on the high and low values. For instance, all the observations made for factor A at level 1 are consolidated for an average value. This applies to all other factors at each level 1 and level 2. The results for the average SNR are presented in table 4. Values of the effect for each parameter in table 4 are obtained by taking the difference of the average SNR between level 2 and level 1. A negative effect implies that the parameter at level 1 has a more dominant effect than that at level 2.

The largest magnitude of effect represents the parameter that has the most significant impact on the temperature distribution inside the human skin. From table 4, this is found to be the disc temperature, factor E, which complements the results obtained from the sensitivity analysis. From table 4, it is also found that the combination of factors at A_{2}B_{1}C_{2}D_{2}E_{1}F_{1}G_{2}H_{1} (subscripts 1 and 2 refer to the low and high levels, respectively) gives rise to a minimum temperature distribution inside the skin during burn.

## 7. Summary

We have successfully developed a model of the human skin to predict the temperature distribution during burn at a steady state. The present study, which includes a parametric study using sensitivity analysis and the Taguchi method, offers an alternative look at how various parameters affect the temperature distribution inside the skin and their order of importance and significance.

The Pennes bioheat equation has been used to describe the flow of heat flux inside the model. Results are obtained numerically using the dual reciprocity method, which show good agreement when compared with the solutions obtained using the more established finite-element method. The good agreement between the results obtained using the boundary element method and the finite-element method, together with the low computational demand of the boundary element method, suggests that the boundary element method may be used as a computationally inexpensive alternative for simulating the temperature distribution inside the human skin during burn.

Sensitivity analysis is carried out for eight different parameters that may have an impact on the temperature distribution inside the human skin. From the results obtained, it is found that the blood perfusion rate of the dermis, ambient convection coefficient and disc temperature are parameters that significantly affect the thermal response of the skin towards heat. Effects of the blood perfusion rate and ambient convection coefficient are negative, where larger values contribute to a drop in temperature. This implies that the blood perfusion rate and ambient convection coefficient are effective in reducing the temperature inside the skin during burn.

Analysis carried out in the present study may be further improved in future studies. From a practical point of view, a steady-state analysis, which has been performed in this paper, raises certain doubts, since it simulates the condition when the skin is heated by the heating disc for an infinitely long time. A more realistic model will include heating of the skin for a duration that would seem practical, and then allowing it to cool down (i.e. a time-dependent problem).

The thermal analysis carried out here is only the first step in the investigation of burn injury. By carrying out a time-dependent analysis, the transient temperature distribution obtained may be used to predict the degree of thermal damage inflicted onto the skin using the thermal damage model of Henriques & Moritz (1947). By doing so, one may attain a better understanding of the physiological response of the skin tissues towards heat.

## Footnotes

One contribution of 9 to a Theme Issue ‘Multi-scale biothermal and biomechanical behaviours of biological materials’.

- © 2010 The Royal Society