Computational Modelling of Dual-Phase Lag Bioheat Process in Cooling of Cutaneous Tissue exposed to High Heating for Burn Injury Prediction

Heat transfer in biological systems is critical in analytic and therapeutic burn applications. Timely burn evaluation and appropriate clinical management are critical to ameliorate the treatment outcome of burn patients. To apply appropriate burn treatment, it is necessary to understand the thermal parameters of the skin. The paper aims to model the non-Fourier bioheat process in the human skin using a multi-domain trivariate spectral collocation method to determine skin burn injury with non-ideal properties of tissue, blood perfusion and metabolism. The skin tissue internal water evaporation during direct heating is considered. Parametric studies on the effects of skin tissue properties, initial temperature, blood perfusion rate and heat transfer parameters for the thermal response and exposure time of triple-layer cutaneous tissues are carried out. The study shows that the initial tissue temperature, the thermal conductivity of the epidermis and dermis, relaxation and thermalisation time and convective heat transfer coefficient are critical parameters necessary for skin burn injury baseline examination. The thermal conductivity and blood perfusion rate also exhibit negligible effects on the burn injury threshold of the cutaneous tissue. The present study is aimed to assist burn evaluation for reliable experimentation, design and optimisation of thermal therapy delivery.


Introduction
Skin burns from short exposure to high heat flux such as exposure to a hot substance, flash fire or laser irradiation occur ubiquitously in our daily lives under different domains including domestic, industrial, and extreme situations such as in military combat. Under these conditions, the heating process is very short typically within the time duration of less than 5s and the heat flux incident of the skin surface can be as high as 83.2kW/m 2 [1]. Burn injury under such heating condition would require burn evaluation and different therapeutic procedures, which might result in lifelong psychological and physical scars [2]. For accurate burn evaluation from such instantaneous thermal injuries, several researchers have intensified research on the heat transfer analysis on living biological tissues.
Heat transfer analysis in living biological tissues is critical for different diagnostic and therapeutic applications, which require temperature changes. Generally, analysing heat transfer processes in biological tissues is arduous due to the complexity of biological structures and composition and involves several processes including thermal conduction in solid tissues, metabolic heat generation, convective heat transfer of blood circulation, perfusion through capillary tubes with the tissues, evaporation, and external heat interaction with the ambient environment. To model the thermal properties of living biological tissues, several bioheat transfer equations have been formulated in the literature including Pennes [3], Weinbaum and Jiji [4], Wulff [5] and Nakayama and Kuwahara t T [6] among others. However, Arkin et al in [7], compared various bioheat transfer models and concludes that the Pennes bioheat model is by far the most practical due to its simplicity for fast prediction of transient temperature profiles. Pennes bioheat model is based on the classical Fourier's law with the assumption that propagation speed of thermal disturbance is infinite which in physical reality is impractical. This is because energy propagation occurs in finite speed as thermodynamic processes in living tissues maintain steady state condition after a certain time lag period [8]. To solve the paradox in the Pennes bioheat equation, Morse and Feshbach [9], Cattaneo [10] and Vernotte [11] independently proposed the linear extension of Fourier's law to establish the finite speed of thermal signal propagation in heat diffusion theory. The resulting thermal wave (TW) model is based on singlephase-lag due to the wave-like behaviour of heat transport in biological structures regarded as C-V constitutive relation as: where the relaxation time, q  is added to capture the micro-scale response in time. Also, since living tissues are highly non-homogeneous and non-uniformity and require relaxation time to accumulate enough energy to transfer to the nearest medium. The TW bioheat model gained its popularity because of its established experimental validations and capability to achieve accurate predictions than the classical Fourier law model. However, the TW bioheat model is restricted to micro-scale response in time and not the micro-scale response in space [12][13][14]. The limitations of the TW bioheat model is overcome by a corresponding model termed dual-phase lag (DPL) model of bioheat transfer, which considers the micro-scale response in space using phase lag in heat flux and temperature gradient [12]. Recently, Rukolaine in [15] explained key unphysical effects of DPL model of heat conduction. The DPL bioheat model is effective to examine the effects of microstructural interactions in the fast transient process of heat transport to analyse the thermal behaviour of skin tissue under various skin conditions [16][17][18].
Additionally, the DPL model characterises microstructural interactions in heat transport and is developed with the first-order Taylor series expansion. Different methodologies have been employed to analyse the DPL bioheat model for single, double, and triple-layer thermal models [19][20][21][22][23][24][25][26]. These methods examine the difference in physiological and thermal properties of the skin. However, since phase change occurs in biological tissues over a wide range, the boundary conditions at the interface between the two adjacent layers would result in highly nonlinear mathematical models. Therefore, for an accurate understanding of the temperature distribution in the human skin when exposed to heat, the use of numerical methods is effective to reduce dependency on timeconsuming and costly experimental analysis. Several researchers have employed different numerical schemes including finite element method [27][28][29], boundary element method [30][31][32], Monte Carlo [33] and finite difference method (FDM) [34] due to adaptation of numerical schemes to adapt to complex geometries and governing equations. The purpose of this paper is to present the non-Fourier thermal modelling of triple-layer cutaneous tissue for the prediction of skin burn injury with non-ideal properties of tissue, metabolism, and blood perfusion. The study proposed a hybrid method using the multi-domain trivariate spectral collocation method to investigate the temperature distribution within cutaneous tissue subject to instantaneous heating of high heat flux. The skin tissue internal water evaporation during direct heating is considered.

Modelling of Human Skin Exposed to Burn
The human skin helps protect the body from microbes and the elements as well as support regulating body temperature under various thermal conditions. The human skin is constructed of two distinct layers: epidermis and dermis, with the hypodermis or subcutaneous layer, which is not part of the skin, but included in the skin structure, as it is attached to the dermis by collagen and elastic fibres, as illustrated in Figure 1. The present work considers the effect of microstructural interactions in the transient process of heat transfer to deal with the paradox of the classical Fourier's model and account for the limitations of the thermal wave model. Under the present consideration, the temperature gradient at any point in the material at time t + τT corresponds to the heat flux vector at the same point at time t + τq, which can be expressed mathematically as: The above model covers a wide range of space and time for physical observations. Rewriting Eq. (3) gives   From the local energy balance, the equation for conservation of energy can be written as: Rearranging Eq. (5) By substituting Eq. (6) in (4) we arrived at Expanding Eq. (7) produces From the above DPL model in Eq. (8), if τq and τT are set to zero, then the Pennes' model is recovered.
However, if τT is set to zero and the first-order Taylor expansion is used only in time, then the thermal wave or hyperbolic model is recovered as: Further, if τT is set to zero and the first-order Taylor expansion is used both in time and space, then Tzou' model in [12] is recovered as: The three-dimensional DPL model in the Cartesian coordinates gives by taking the thermal conductivity as constant, the three-dimensional DPL model in the Cartesian coordinate gives The two-dimensional DPL model in the Cartesian coordinates gives Taking the thermal conductivity as constant, the three-dimensional DPL model in the Cartesians co-ordinates produces The one-dimensional DPL model with varying thermal conductivity is expressed as Then the one-dimensional DPL model with constant thermal conductivity is expressed as:

Modelling Evaporation during the High Heating
During heating of tissues to high temperature, the skin tissue moisture, normally around 70-75% evaporates, absorbing latent heat. This is because at 100°C, during ablation, water vaporisation in the tissues occurs leading to dehydration of the tissues. As the temperature increases (>100°C), the continuous vaporisation and dehydration of the tissue initiate the carbonisation process of the tissues [35][36][37]. These biological processes are fundamental for the analysis of skin burn at high temperature and without due consideration of these processes, results obtained from the models would differ significantly from experimental results. To this end, the vaporization terms are included in the above models to accommodate the phase change due to evaporation. Under such condition, the enthalpy-based model can be adopted [38]. Therefore, without loss of generality, in this work, a simple method to incorporate the simple water-related processes into the thermal models is employed to improve the ablation models at high temperatures. In the simple model, the rate of vaporisation is modelled as established in [39] as  and w  represents the latent heat of vaporisation for water (2260kJ/Kg) and tissue water density respectively.
The tissue water density is a function of temperature and Eq. (18) can be written using Chain rule as: (19) where 0 W T   < for evaporation

Two-dimensional DPL Model for the Triple Layer Skin Tissue
Consider the triple-layer skin tissue in Figure 2, the two-dimensional triple-layer DPL thermal model can be written as The two-dimensional, three-layer DPL thermal model can be written as: 2  2  2  3  2  ,  ,  , , which can be written as   2  2  2  2  3  2  2  2  2  2  2  2  2  ,  ,  , ,  T  T  T  T  T  T  T  T  T  T  T  c is the effective heat capacity The present study considers a skin insult scenario which consists of direct contact to the with a heated disk of 1cm diameter that is maintained at a defined constant temperature for a specified time. The surface area peripheral of the disc is assumed to experience convective heat transfer with an ambient air temperature of 25 o C during the burn process. After completion of the insult, the disc was removed and the entire skin surface is cooled by natural convection of the surrounding air [31,32]. The initial, boundary and interlayer conditions are stated as:

Initial conditions
At the initial condition, the temperature of the skin tissue is equal to the blood. Therefore,

Initial conditions
At the initial condition, the temperature of the skin tissue is equal to the blood and is expressed as: (28a) (28b) (28c)  The developed skin burns model is used to determine the required time to generate the different degrees of burns on the human body based on the understanding of the classification in Table 1.

Results and Discussion
The above-developed solutions are simulated and the effects of various thermal and flow parameters on temperature for the assessment of burn injury are investigated. The effects of thermal conductivity of the triple layer of the skin for the prediction of burn injury threshold are presented in Figures 3 -5. In Figure 3, it is shown that as thermal conductivity of the epidermis increases by a factor of 2, the injury baseline shifts towards the lefthand side (LHS) of the base value (0.210W/mK). This is because, as the thermal conductivity of the epidermis increases under continuous temperature exposure, the thermal resistance reduces resulting in increased heat penetration of the tissue. Under this condition, the time required to reach a second-degree burn injury situation becomes minimal, which implies that the higher the thermal conductivity of the skin layer tissue, the lower the degree of burn injury. Moreover, the opposite trend is observed as the thermal conductivity decreases by a factor of 2, which results in the injury baseline to shift towards the right-hand side (RHS) of the base value (0.210W/mK). This is because as thermal conductivity reduces, the thermal resistance increases resulting in reduced heat penetration of the tissue. Consequently, under such burns condition, more time would be required to reach a second-degree burn injury. Figure 4 shows that as the thermal conductivity of the dermis increases by a factor of 2, the injury baseline shifts towards the RHS of the base value, i.e., 0.370W/mK. This phenomenon occurs, as an increase in the dermis thermal conductivity causes heat at the epidermis-dermis interface to be readily transferred to the deeper tissue and consequently resulting in prolonged time required to reach the injury baseline. However, an opposite trend is observed when the thermal conductivity of the dermis is decreased by a factor of 2, resulting in the injury baseline to shift towards the LHS of the base value (0.370W/mK). Moreover, the above finding implies that the thermal conductivity of the dermis reduces as the temperature of exposure heat decreases. From Figure 5, the thermal conductivity of the hypodermis exhibits no significant effect on the prediction of the injury baseline. Figure 6 illustrates the effect of blood perfusion rate on the prediction baseline. Under constant volumetric capacity of the blood, the effects of no perfusion (0ml/100g/min), half of the maximum dilatation (75ml/100g/min) and maximum blood flow (150ml/100g/min) of the skin vessels. In Figure 6, the blood perfusion rate exhibits negligible effects on the dermis and invariably does not affect the burn injury prediction in the dermis layer. The effects of initial skin tissue temperature on the burn injury baseline prediction are presented in Figure  7. From Figure 7, the initial skin tissue temperature exhibits significant effects on skin exposed to burn injury. This is because the warmth of the skin tissue increases as the temperature experienced throughout the heating and cooling periods increases resulting in an increased degree of burn injury. The effect of convective heat transfer coefficient for the prediction of burn injury baseline is illustrated in Figure 8. In Figure 8, it is shown that the heat transfer coefficient exhibits significant effects on the burn injury baseline for minimal exposure time. However, the effect of the convective heat transfer coefficient becomes negligible, with prolonged exposure time.   This is because, under prolonged exposure of the skin layers to burn, the burn injury is dominated by conductive heat transfer in the tissue rather than the convective heat transfer at the surface of the skin. Moreover, as the convective heat transfer coefficient increases by a factor of 2, the injury baseline shifts towards the LHS of the base value (1250W/m 2 K). Under such condition, minimal time is required to reach the second-degree injury baseline. However, when the convective heat transfer coefficient is decreased by a factor of 2, the injury baseline shifts towards the RHS of the base value (1250W/m 2 K), which shows that as convective heat transfer coefficient decreases, more time is needed for the skin burn to reach a second-degree burn injury. Figure 9 shows the initial tissue temperature profiles for the three-layer skin tissue whilst the injury baseline using a uniform tissue temperature of 34°C and the initial temperature profile for two and three-layer slabs is presented in Figure 10. Figures 9 and 10 shows that the selection of initial temperature 34°C for all the three layers is highly reasonable. It is worth noting that uniform initial temperature can be applied for all three layers for a more accurate analysis of the cutaneous tissue burn injury. Moreover, the skin injury layer baseline range of combinations of different values of thermal properties is presented in Figure 11. In Figure 11, the effect of a worst-case combination of parameters is demonstrated where the upper baseline represents a combination of low thermal properties of skin components, low heat transfer coefficient and low initial tissue temperature. The lower baseline represents a combination of high values of these parameters.

Effects of Relaxation time and Thermalisation time
The effects of relaxation time and thermalisation time on the skin temperature during the cooling process is presented in Figure 12. The simulation shows that the peak temperature predicted by the classical bioheat equation is higher than the DPL during the heating process. However, the temperatures from all considered cases converge during the cooling process. This is because the finite speed of wave propagation and the peak temperature predicted by the DPL occurs with time lags. The time lag causes longer thermal dissipation, i.e., cooling by the heat conduction of tissue and by the blood perfusion period of peak temperature. The predicted higher temperature by the classical bioheat model than the DPL model implies that the accumulation of the thermal damage is overestimated since the level of the accumulative thermal damage depends primarily on the peak temperature. Moreover, the relaxation and thermalisation time, as shown in Figure 12, enhances heat diffusion, which enables a reduction in the thermal damage of the skin tissues. This effect is significant as the heat transfer process approaches the peak temperature during heating. However, the curves are close and almost similar values are predicted during cooling for the relaxation time and thermalisation time considered. a. Maximum temperature distribution in the skin tissue before cooling b. Maximum temperature distribution in the skin tissue after 4sec of cooling c. Maximum temperature distribution in the skin tissue after 7sec of cooling Figure 13. Temperature distribution of the skin tissue before and after cooling The temperature distribution of the skin tissue at depth 0.006m and length 0.024m for the period of 8sec is presented in Figure 13. In Figure 13, the dimensionless length and depth of the tissue is shown and the negative values in the dimensionless depth indicate depth zero i.e., below the skin surface. The dimensional length ranges from -1 to +1 to depict symmetry as stated in the boundary conditions. The dimensional length ranges from -1 to +1 to depict symmetry as stated in the boundary conditions. Further, the results of the present study are compared with the different established works in the literature including Henriques [42], Fugitt et al. [43], Stoll and Greene [44], Takata [45] and Wu [46] as presented in Figures 14-19. Figure 14 demonstrates the comparison of the various results for the effects of surface temperature on burn injury. Figures 15 -17 demonstrates the required time to reach the first, second and third-degree burn injuries when the skin tissue temperature is maintained at 50°C, 70°C, and 90°C respectively. The effects of tissue depth and surface temperature on the degree of burn and burn injury are illustrated in Figures 18 and 19. It is shown that the thermal damage value or the value of the burn injury is minimal for low tissue temperature but rapidly increases with temperature above 50°C.In Figures 14-19, it is shown that the present study quantitatively agrees with previous studies for low-level damage at a temperature below 50 o C. However, there is significant variation among the skin burn damage models results for the high exposure temperatures. In fact, at a temperature above 50 o C, only Henriques [42] and Takata [45] corresponds reasonably for all temperatures evaluated whilst in agreement with the results of the present study. Figure 20 demonstrate the validation of the present study with the experimental results of Takata [45]. From the comparison, the present study shows a significant agreement with the experimental results in [45], which demonstrates the reliability of the obtained results of the study.

Conclusion
The non-Fourier modelling of DPL bioheat process in human skin exposed to high heating for burn evaluation and prediction is presented. The effects of cutaneous tissue properties, initial temperature, blood perfusion rate, relaxation and thermalisation time and heat transfer parameters for the thermal response and exposure time of triple-layer cutaneous tissue is evaluated. The present work employs a multi-domain trivariate spectral collocation method to solve the aforementioned effects of the DPL bioheat process. The study shows that an increase in the thermal conductivity of the epidermis decreases the thermal resistance and readily causes an increased heat penetration of the cutaneous tissue. The finding implies that a high thermal conductivity of the tissue lowers the degree of burn injury. The study also shows that increase in dermis thermal conductivity results in a prolonged time to reach the injury threshold and the blood perfusion rate exhibits no net effect on the prediction of burn injury in the dermis layer. However, the initial skin tissue temperature exhibits significant effects on burn injury exposure.
The relaxation and thermalisation time play fundamental roles in the analysis of burn injury. The relaxation time and thermalisation time enhances heat diffusion and reduces the thermal damage in the skin tissues. In addition, the heat transfer coefficient plays significant effects on burn injury for small exposure time. The effect of the convective heat transfer coefficient becomes minimal for prolonged exposure time. However, for accurate skin burn evaluation, it is convenient to use uniform initial temperature for the triple-layer cutaneous layers. The most dominating factors in the burn injury follow the order: relaxation and thermalisation time, initial tissue temperature followed by the epidermis layer thermal conductivity, dermis layer thermal conductivity and convective heat transfer coefficient. The present work would help in the quantification of skin burn, and for clinicians and biomedical engineers to experiment, design, characterise and optimise strategies for delivering thermal therapies.