#### Abstract

A semianalytical method based on the separation of variables is presented to solve the problem of heat transfer in multilayer skin tissues. This method is used to predict the temperature field and thermal damage of skin tissues, which are subjected to the time-varying laser heating and fluid cooling. The effect of convective heat transfer coefficient of the skin surface on temperature field and thermal damage is studied. The temperature field and thermal damage of skin tissue with different laser heating cycles are calculated. The relationship between the laser heating period and the temperature fluctuation is obtained. The effects of the peak time of laser heating on temperature and the thermal damage of skin tissues are investigated. The relationship between the peak time of laser and the temperature fluctuation of skin tissue is drawn. The results calculated in this paper can provide some guidance for hyperthermia of skin tissue. The presented semianalytical method can be applied to solve the general heat transfer problems of multilayer structures.

#### 1. Introduction

Thermotherapy has been widely used in the treatment of skin tissue diseases. Laser, microwave, and radio frequency are often used to irradiate the surface of skin tissue during thermotherapy [1–3]. Accurate prediction of skin tissue temperature and thermal damage during hyperthermia plays an important role in thermotherapy [4]. More and more researchers have studied the heat transfer and thermal damage of the skin tissue. Important research can be summarized as follows.

Deng and Liu gave several closed-form analytical solutions to the bioheat transfer problems which are based on Green’s function method [5]. Xu et al. studied the heat transfer of skin tissue by using pence model; they established a mathematical model of heat transfer of three-layer skin tissue and solved the model by using finitely different method [6]. Liu et al. established DPL model for heat transfer of multilayer skin tissue; the Laplace transform method was used to solve this model [7]. Tian et al. established a skin-fabric five-layer heat transfer model; a semianalytical method was employed to solve this model [8]. Ahmadikia et al. used analytical methods based on Laplace transform to study the distribution of temperature field in skin tissue under constant or transient heat flux conditions; parabolic and hyperbolic skin tissue heat transfer equations were studied [9]. In the study of Kumar et al., a dual-phase-lag model of bioheat transfer had been studied by using Gaussian distribution source term with the most generalized boundary condition during hyperthermia treatment [10]. Lin solved the analytical solution of heat transfer of single skin layer which was based on the separation of variable methods; Fourier and non-Fourier cases were studied [11]. Hooshmand et al. developed an analytical solution method for solving the problem of heat transfer in biological tissue, generalized dual-phase-lagging (DPL) model was studied, and the calculation results were compared with the classical DPL model [12]. Kumar et al. dealt with the study of heat transfer and thermal damage in triple layer skin tissue using fractional bioheat model. An implicit finite difference scheme is obtained by a new method [13]. Kumar and Rai studied the thermal characteristics in skin tissue using time fractional dual-phase-lag bioheat transfer (DPLBHT) model subjected to Dirichlet boundary condition. They solved this problem using the finite element Legendre wavelet Galerkin method. They treated skin tissue as a single layer [14]. Li et al. predicted the temperature of different biological tissues using Fourier and non-Fourier heat transfer models. And the temperature response of different tissues was measured. The experimental results were compared with those predicted by the model [15]. Hobiny and Abbas used the Laplace transform method to solve the problem of biological tissues subjected to moving heat source [16]. Ghanmi and Abbas derived an analytical solution for the fractional transient heating within the skin tissue during hyperthermia [17]. Sur et al. used Laplace transform to study the thermal behavior of skin tissue when it was subjected to moving heat sources [18]. Sharma and Kumar established a nonlinear dual-phase-lag heat transfer model of biological tissue. The finite element Runge–Kutta method was used to solve this model. The temperature field and thermal damage of skin tissue were studied using the numerical solution of the model [19].

Skin tissue is the largest organ of the human body, including epidermis, dermis, and subcutaneous tissue [20]. Thermophysical parameters of different layers of skin tissue are different. So it is difficult to solve the heat transfer model of skin tissue by analytical method [9, 11, 21, 22]. Numerical methods were often used to solve this problem [23]. Several papers used analytical methods, but all of them treated skin tissue as a single layer [12]. In this paper, a semianalytical method is developed to solve heat transfer in multilayer skin tissues. A special method of separation of variables is used to derive the analytical solution. In the derivation process, the time domain is firstly divided into several small time periods. Then the temperature function at the interface of different layers of skin tissue is assumed to be linear in a small period of time. After that, the analytical solution of separated variables of skin tissue heat transfer is deduced in a small period of time. According to the condition that the temperature at the interface is equal and the heat flux at the interface is continuous, the temperature at the interface is calculated. Temperature field of skin tissue during laser heating and fluid cooling is obtained by the semianalytical solution. The effects of several parameters such as convective heat transfer coefficient, period, and peak time of laser heating on temperature field are studied. These parameters can be used to control the temperature of the skin tissue and the thermal damage in the appropriate range.

#### 2. Mathematical Model

Human skin tissue consists of epidermis, dermis, and subcutaneous. The geometric model of skin tissue is shown in Figure 1. The corresponding mathematical model can be established by combining the boundary conditions. The governing equation can be expressed as (1) [2, 24–27]:

The boundary conditions can be expressed as follows:

At *x* = 0,

At *x* = *L*_{3},

The initial condition is

The corresponding dimensionless governing equation is

The dimensionless boundary conditions are as follows:

At = 0,

At = 1,

The dimensionless initial condition is

#### 3. Solution Method

If the value of is small enough, in time domain [0 ], can be expressed as [28]

Equations (7) and (9)–(11) show that the boundary conditions on both sides of each layer of skin tissue are known in a short time interval [0 ]. A extended shifting variable method is developed by Lee and Lin, which can be used to solve the temperature field of each layer of skin tissue [29].

Assume

Here, (*ξ*), *i* = 1, 2, 3; *m* = 1, 2 are the shifting functions to be specified and (*ξ*, *τ*) is the transformed function. Substituting (12) into (5), the two following subsystems are obtained. The first subsystem is expressed in terms of the transforming variables as follows; the transform governing equation is

If *i* = 1, substitute (12) into (6) and (10):

In order to homogenize the boundary conditions when separating variables, (17) and (18) need to be satisfied:

So,

Assume

If *i* = 2, 3, substitute (12) into (10) and (11):

In order to homogenize the boundary conditions when separating variables, (24) and (25) need to be satisfied:

So,

Assume

The characteristic equation of (13) is

Assume

If *i* = 1, substitute equation (30) into (17) and (18):

Equation (31) is a second-order homogeneous ordinary differential equation. Its characteristic function is as follows:where is the positive solution of the following equation:

If *i* = 2, 3, substitute (30) into (24) and (25):where is the positive solution of the following equation:

Obviously, satisfies the following equation:

Assume

Substitute (40) into (13) and integrate it in []:

The solution of (41) is

#### 4. Thermal Damage

In addition to predicting the temperature field of skin tissue during thermotherapy, it is often necessary to predict its thermal damage. Formula (47) is developed by Henriques and Moritz which can be used to calculate thermal damage of skin tissue [1]:where represents the thermal damage factor of skin tissue. represents the activation energy. represents the temperature of the skin tissue. *R* represents the universal gas constant. According to the thermal damage factors, the thermal damage is divided into three levels in the literatures. , the first level of thermal damage; , the second level of thermal damage; , the third level of thermal damage.

#### 5. Numerical Results

When the surface of skin tissue is irradiated by pulsed laser, the variation of laser energy with time can be shown in Figure 2. Assume the initial temperature of the skin tissue *T* (*x*, 0) = 37°C. At *x* = *L*_{3}, *T* (*L*_{3}, 0) = 37°C. The thermophysical parameters of skin tissue can be shown in Table 1 [4].

Figure 3 shows the temperature history at the skin surface, ED interface, and DS interface in a laser heating cycle. It is found that the temperature rises first and then decreases with time at *x* = 0, *x* = 0.1 mm. If *h* = 7 W/m^{2}°C, the temperature peaks appear at *t* = 2 s and *t* = 2.4 s, respectively. If *h* = 700 W/m^{2}°C, the temperature peaks appear at *t* = 1.4 s and *t* = 1.7 s, respectively. So we can speculate that the larger the value of *x*, the shorter the time required for temperature peak. If *h* = 7 W/m^{2}°C, at *x* = 1.6 mm, the temperature first remains constant and then continues to rise. When *t* = 10 s, the peak value of temperature still does not appear. If *h* = 700 W/m^{2}°C, at *x* = 1.6 mm, because this location is far away from the heating source, temperature varies a little with time. It can be inferred from the comparison of solid and dotted lines that the larger the value of *h*, the lower the temperature. Investigating its reason we can know that the larger the value of *h* is, the more heat is taken away by water cooling on the surface of skin tissue, so the lower the temperature of skin tissue is.

Figure 4 shows the temperature distribution in the *x* direction at different times. As shown in the picture, if *h* = 7 W/m^{2}°C, *t* = 1 s, the temperature of the skin tissue surface is the highest, which is 65°C. The temperature distribution of skin tissue decreases in the *x*-axis. Comparing the temperature distribution curves of *t* = 0.5 s and *t* = 1 s, it is found that the temperature of *T* = 1 s is higher than that of *t* = 1 s. This is because the heat flux of laser heating keeps rising in 0 to 1 s. Because the heat flux of laser heating begins to decrease after *t* = 1 s, the temperature distribution at *t* = 5 s and 10 s is much lower than that at *t* = 1 s. Comparing the temperature distribution of *h* = 7 W/m^{2}°C and *h* = 700 W/m^{2}°C, it is found that, at the same time, the temperature distribution of skin tissue at *h* = 700 W/m^{2}°C is significantly lower than that at *h* = 7 W/m^{2}°C. Because the cooling effect is greater than the heating effect, the temperature of skin tissue surface is lower than 37°C at *t* = 5 s and *t* = 10 s.

Figure 5 demonstrates the relationship between heat flux and heating time at different locations of skin tissue. As can be seen from the figure, when *h* = 7 W/m^{2}°C, the heat flux density on the surface of skin tissue (*x* = 0) increases first and then decreases, which is in accordance with the trend of surface heating heat flux. At the skin-dermis interface (*x* = 0.1 mm), the heat flux is slightly lower than that at the skin surface. At the epidermal-subcutaneous interface (*x* = 1.6 mm), the heat flux increases first and then decreases, but the range of variation is very small. On the left side of the skin tissue (*x* = 6 mm), the heat flux is 0 and remains unchanged. When *h* = 700 W/m^{2}°C, the heat flux of heating is less than that of cooling in the initial period of time. The heat flux of skin surface (*x* = 0) and dermal-epidermal interface (*x* = 0.1 mm) changed from negative to positive, continued to increase, and then gradually decreased. Because the boundary condition we set is that the right side of the skin tissue is adiabatic, on the left side of the skin tissue (*x* = 6 mm), the heat flux is 0 and remains unchanged too.

Figure 6 demonstrates the relationship between thermal damage and heating time at different locations of skin tissue. As can be seen from the figure, thermal damage of skin tissue increases with time accumulation. The surface of skin tissue has the greatest thermal damage. If *h* = 7 W/m^{2}°C, the thermal damage value of skin tissue surface is close to 100, reaching the second-degree thermal damage, and the thermal damage value at the junction of epidermis and dermis exceeds 0.53, reaching the first-degree thermal damage. Because the greater the value of *h*, the better the heat dissipation and the smaller the thermal damage of skin tissue, the thermal damage at *h* = 700 W/m^{2}°C is six orders lower than that at *h* = 7 W/m^{2}°C. It can be seen from the figure that although the heat flux of laser heating increases first and then decreases, the heat damage increases continuously. As the heat flux of laser heating decreases, the trend of increasing skin tissue thermal damage slows down. In order to analyze the effect of *R* and *H* on thermal damage, a variance analysis was carried out. Table 2 shows the calculated results. By calculation, *R*^{2} (*x*) = 0.163812, *R*^{2} (*h*) = 0.318886, *R*^{2}_{(interaction)} = 0.318877. Therefore, it can be inferred that the influence level of *x* on thermal damage is lower than the average and the influence level of *h* on thermal damage is above average. These two factors and their interaction have a great influence on thermal damage.

Figure 7 demonstrates the effect of heating cycle Tc on the temperature history at skin surface. As can be seen from the figure, because heating and cooling are the same in the first five seconds, so the temperature curve coincides completely. If *T*_{C} = 5 s, the first pulse of laser heating is completed after 5 seconds, the heat flux of heating no longer decreases, and the second laser heating pulse starts. Therefore, the temperature on the surface of skin tissue starts to rise when *t* = 5 s, and the rising process is fluctuant, and the fluctuating period is *T*_{C} = 5 s. After five seconds, the shorter the heating cycle is, the faster the temperature rises. Temperature fluctuates with heating period. After 30 seconds, the amplitude of temperature fluctuation tends to be balanced. The shorter the heating period *T*_{C} is, the greater the temperature fluctuation amplitude *T*_{s,max} is.

**(a)**

**(b)**

Figure 8 shows the effect of heating cycle *T*_{C} on the temperature history at the interface between epidermis and dermis. Comparing with Figure 7, it can be seen that the temperature fluctuation at the interface of epidermis and dermis is smaller than that at the skin surface. It is the same as the result shown in Figure 6, The shorter the heating period *T*_{C} is, the greater the temperature fluctuation amplitude *T*_{s,max} is. When *T*_{C} = 5 s, the amplitude of fluctuation reaches the maximum value; at the junction of epidermis and dermis, *T*_{ed1,max} = 82.5°C, it is 9.5°C lower than the surface temperature of skin tissue. When *T*_{C} = 15 s, the amplitude of fluctuation reaches minimum value; at the junction of epidermis and dermis, *T*_{ed3,max} = 65.2°C, it is 7.9°C lower than the surface temperature of skin tissue.

**(a)**

**(b)**

Figure 9 shows the effect of heating cycle *T*_{C} on the temperature history at the interface between dermis and subcutaneous. Because the location is far from the heat source, the range of temperature fluctuation is much smaller than the skin surface and the interface between epidermis and dermis. It takes a long time for temperature fluctuation to reach a steady state. As can be seen from the figure, when *T*_{C} values are different, the time required for temperature fluctuation to reach stability is 75 s, 100 s, and 150 s, respectively. It can be inferred that the greater the value of temperature *T*_{C} is, the less time it takes for temperature fluctuation to reach stability. It can also be seen from the figure that the larger *T*_{C} value is, the larger temperature fluctuation range is.

**(a)**

**(b)**

Figure 10 shows the effect of heating cycles on the thermal damage at different locations of skin tissue. Because the surface of skin tissue is close to the heat source, the thermal damage on the surface of skin tissue (*x* = 0) is much greater than the interface between dermis and epidermis (*x* = 0.1). With the accumulation of heating time, the thermal damage of skin surface in all three cases exceeds 10^{4}, reaching the third-degree burn. The thermal damage with *T*_{C} = 5 s is four orders higher than the thermal damage with *T*_{C} = 10 s. It can be inferred from the figure that the smaller the heating period is, the more serious the thermal damage of skin tissue is. The thermal damage at the interface between epidermis (*x* = 0.1 mm) and dermis exceeds 10^{4}, reaching the third-degree burn when *T*_{C} = 5 s or *T*_{C} = 10 s. When *T*_{C} = 15 s, the thermal damage value at the interface between epidermis and dermis exceeds 100 to reach second-degree burn. Because the dermis-subcutaneous interface (*x* = 1.6 mm) is far away from the heat source, no burn occurs with *T*_{C} = 10 s and *T*_{C} = 15 s, and the first-degree burn occurs with *T*_{C} = 5 s. In order to study the influence of *T*_{C} and *x* on thermal damage, variance analysis was carried out. Table 3 shows the calculation results. By calculation, *R*^{2} (*x*) = 0.165575, *R*^{2} (*T*_{C}) = 0.167528, and *R*^{2}_{(Interaction)} = 0.3308943. Therefore, it can be inferred that the influence level of *x* and *h* on thermal damage is lower than the average, respectively. These two factors and their interaction have a significant effect on thermal damage.

Figure 11 indicates the relationship between laser heating cycle and temperature fluctuation amplitude. *T*_{s,max} and *T*_{s,min} represent the maximum and minimum values of temperature fluctuation amplitude of skin tissue surface, respectively. *T*_{ed,max} and *T*_{ed,min} represent the maximum and minimum values of temperature fluctuation amplitude of the interface between epidermis and dermis, respectively. *T*_{ds,max} and *T*_{ed,min} represent the maximum and minimum values of temperature fluctuation amplitude at the interface between dermis and subcutaneous, respectively. It can be seen from the figure that the amplitude of temperature fluctuation increases first and then decreases with the increase of laser heating period. The range of temperature fluctuation increases with the increase of laser heating period. Skin tissue surface temperature fluctuation range is the largest, the epidermis and dermis interface is the second, and the dermal and subcutaneous interface is the smallest. With the increase of laser period, the minimum amplitude of temperature fluctuation becomes closer and closer; when *T*_{C} > 16 s, the lowest temperature of skin tissue is equal at all locations.

Figure 12 shows the effect of time *t*_{p} on skin surface temperature history. It can be seen from the figure that the skin tissue surface temperature fluctuates and rises at the beginning time and then tends to fluctuate steadily. If *t*_{p} = 3 s, the fluctuation range of skin surface temperature is 80–107°C. If *t*_{p} = 1 s, the fluctuation range of skin surface temperature is 45–75°C. If *t*_{p} = 0.5 s, the fluctuation range of skin surface temperature is 42–63°C. We can infer that the larger the *t*_{p} value is, the larger the maximum value of skin tissue surface temperature fluctuation is and the smaller the minimum value of skin tissue surface temperature fluctuation is.

**(a)**

**(b)**

Figure 13 shows the effect of peak time *t*_{p} on the temperature history at the interface between epidermis and dermis. Comparing with Figure 12, we can see that the temperature of this interface is lower than that of the skin surface. The trend of fluctuation is consistent with Figure 12.

**(a)**

**(b)**

Figure 14 shows the effect of peak time *t*_{p} on the temperature history at the interface between dermis and subcutaneous. Because this location is far away from the heat source, the temperature of this location is much lower than that of the skin tissue surface and the interface between the epidermis and dermis. The temperature at this position fluctuates with the heat flux heated by laser.

**(a)**

**(b)**

Figure 15 indicates the relation between peak time *t*_{p} and temperature fluctuation Amplitude. The meaning of *T*_{s,max}, *T*_{s,min}, *T*_{ed,max}, and *T*_{ed,min} is the same as shown above_{}. It can be seen from the figure that, with the increase of peak time *t*_{p}, the amplitude of temperature fluctuation becomes larger and larger at all positions of skin tissue. The temperature fluctuation range of skin surface and epidermis-dermis interface increases with the increase of peak time *t*_{p}. The temperature fluctuation range of dermal-subcutaneous interface increases first and then remains constant with the increase of *t*_{p}. The larger the *t*_{p} value is, the smaller the range of temperature fluctuation is.

Figure 16 indicates the effect of peak time *t*_{p} on the thermal damage history at different locations of skin tissue. It can be seen from the figure that the greater the *t*_{p} value, the more severe the thermal damage of the skin tissue. With *t*_{p} = 3 s thermal damage value of skin surface and epidermal-dermal interface exceeds 10^{4}, reaching the third-degree burn, and thermal damage value of dermal-subcutaneous interface exceeds 1 to reach a second-degree burn. If *t*_{p} = 0.5, the thermal damage on the skin surface (*x* = 0) exceeds 100 to reach a second-degree burn. The thermal damage on the skin surface (*x* = 0.1 mm) exceeds 1 to reach a first-degree burn. Because the dermal-subcutaneous interface is far away from the laser heating heat source, there is no burn. If *t*_{p} = 1, the thermal damage on the surface of skin tissue exceeds 10^{6}, reaching the third-degree burn. The thermal damage on the epidermis-dermis interface of skin tissue exceeds 10^{6}, reaching the second-degree burn. The thermal damage on the dermal-subcutaneous interface of skin tissue is less than 0.53, lower than the first-level burn standard. In order to study the influence of *t*_{p} and *x* on thermal damage, variance analysis was carried out. Table 4 shows the calculation results. By calculation, *R*^{2} (*x*) = 0.165721, *R*^{2} (*t*_{p}) = 0.167244, and *R*^{2}_{(Interaction)} = 0.3325. Therefore, it can be inferred that the influence level of *x* and *t*_{p} on thermal damage is lower than the average, respectively. These two factors and their interaction have a significant effect on thermal damage.

#### 6. Discussion

It is of great significance to accurately predict the temperature field and control the thermal damage during skin tissue hyperthermia. In this paper, a special separation variable method is used to derive the semianalytical solution of temperature field in skin tissue hyperthermia. Compared with the traditional finite element method, finite difference method, Laplace transform method, and so on, this method has clearer physical meaning and higher calculation efficiency. The temperature field and thermal damage of skin tissue during laser hyperthermia are predicted by the analytical solution derived in this paper. Several parameters that can affect the temperature field and thermal damage of skin tissue during heating have been studied. Through analysis we can find that the convective heat transfer coefficient between the skin tissue surface and the cooling fluid has a greater effect on the temperature field and thermal damage. We can change the injection speed of the cooling fluid to change the convection heat transfer coefficient and then control the temperature and thermal damage of the skin tissue in an appropriate range. In addition, the heating period and peak time of the laser also have a greater impact on the skin tissue temperature field and thermal damage during hyperthermia. Adjusting the parameters such as laser heating period or peak time can control the skin tissue temperature to fluctuate within a stable range. This operation can achieve the desired effect of hyperthermia without thermal damage to skin tissue.

#### 7. Conclusion

In this paper, a semianalytical method is presented to solve the heat transfer problem of multilayer skin tissues. The temperature field and thermal damage of skin tissue under time-varying laser heating and fluid cooling are obtained by the presented method. The effects of several parameters on skin tissue temperature and thermal damage are investigated. The main phenomena are revealed as follows:(a)The smaller the convective heat transfer coefficient *h* is, the faster the temperature of skin tissue rises(b)The larger the convective heat transfer coefficient *h* is, the smaller the heat flux of skin tissue is(c)The smaller the heating period is, the more serious the thermal damage of skin tissue is(d)The range of temperature fluctuation decreases with the increase of laser heating period *T*_{C}(e)The greater the *t*_{p} value, the more severe the thermal damage of the skin tissue

In this paper, the Fourier model of heat transfer in biological tissue is used to study the heat transfer in skin tissue. In recent years, some literatures have used non-Fourier model to simulate the heat transfer in biological tissue. Our future research direction is to use the method of this paper to solve the non-Fourier heat transfer model of biological tissue and then to study the heat transfer of skin tissue.

#### Nomenclature

i: | Layer number of skin tissue: i = 1, epidermis; i = 2, dermis; i = 3, subcutaneous |

: | Thermal conductivity of i-th layer |

: | Temperature |

: | Density of i-th layer |

: | Specific heat capacity of i-th layer |

: | Time |

: | Thickness of i-th layer |

: | x-axis coordinates, |

: | Blood perfusion rate |

: | Density of epidermis blood |

: | Specific heat capacity of blood |

: | Metabolic heat generation in the skin of i-th layer |

: | Heat source due to external heating i-th layer |

: | Heat convective coefficient |

: | Temperature of cooling fluid |

: | Temperatures of arterial blood |

: | Heat flux function |

: | Dimensionless |

: | Dimensionless T |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless t |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless |

: | Dimensionless . |

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors are grateful to the financial body, Department of Education of Guangdong Province (associated grant number: 2017KTSCX218).