Theoretical analysis on thermal treatment of skin with repetitive pulses.

Thermal ablation is an efficient method of medical treatment, such as cancer therapy, wound closure, laser cutting, freckle removal and other treatments. In order to guarantee the curative effect and the safety of the patients, the thermal response of the tissue which is subjected to the heat source need to be carefully monitored. However, it is too difficult to achieve real-time monitoring on the full-field temperature. In the present study, efforts were made to build up a theoretical model for the prediction of the thermal response in the human skin. The Dual-Phase-Lag (DPL) bio-heat transfer model and the Henrique's burn assessment model were employed to describe the interaction of multi-pulse heat source and the skin. The repeated multi-pulse laser is a common heat source in the thermal treatment and the thermal responses of the skin would be complicated under the common effects of the non-Fourier effects and the multi-pulse source. The Green's function approach was used to solve the governing equations analytically. The closed-form solution for the temperature distribution of the skin was obtained and the thermal damage was estimated based on the temperature results. The influences of the biological parameters (the phase lags of the heat flux and the temperature gradient) and the heat source parameters (the pulse number and the duty ratio) on the temperature distribution, the burn degree and the irreversible burn depth of the irradiated region were discussed.

Introduction

Thermal ablation is an efficient method of medical treatment, such as cancer therapy, wound closure, laser cutting, freckle removal and other treatments. In thermal ablation, human tissues are exposed to different kinds of heat producers, including microwave, radio-frequency, infrared radiation, magnetically excitable thermo-seeds and other thermal sources[1]. However, the high temperature could hurt the healthy cells by loosening the cell membranes and denaturing the proteins[2]. In order to guarantee the curative effects and the patients’ safety, the temperature distribution in human tissue should be strictly controlled. Given that real-time monitoring on the full-field temperature cannot be achieved at present, it is extremely significant to make precise prediction of the temperature distribution before thermal treatment. The mechanism of the thermal responses in human tissue during the hyperthermia need to be investigated in detail.

Among the existing literatures, Pennes[3] was the first to study the heat conduction in the living biological tissues by proposing a bio-heat transfer model on the basis of Fourier’s law. Since then, the Pennes model has been widely applied by researchers. For example, Yue[4] developed a one-dimensional steady-state bioheat transfer model of living tissues in cylindrical coordinates and derived the analytical solution with the usage of Bessel’s equation. Shih[5] investigated the thermal response of a semi-infinite biological tissue due to a sinusoidal heat flux at the skin. Laplace transform method was employed to obtain the semi-analytical solution for the governing equations. Fu[6] explored the thermal response inside the skin in the process of thermal treatment on the basis of Pennes model. The influences of the tumor shape, including the location, geometry and size, on the temperature distribution were taken into account. Singh[7] investigated the thermal effects on the breast tumor during the radio-frequency ablation. The effects of temperature-dependent parameters of the multi-layer breast model, such as electrical and thermal conductivities, were considered. Bhowmik[8] carried out a numerical study on the thermal response of skin with different fat thickness. The Pennes model and the numerical simulation method were employed to derive the temperature distribution of the skin tissue. Marchandise[9] probed the feedback from human nociceptive system on the thermal stimuli with finite element simulation and Fourier’s law employed. The influences of the water-content rate in both hairy skin and glabrous skin were discussed emphatically.

All the above researches were based on the traditional Fourier’s law which assumed that thermal propagation velocity is infinite. However, this assumption was found not suitable for some issues in which the non-Fourier effects shows up, such as heat transfer problems in low-temperature environment and in the biological tissue[10-12]. To make the thermal conduction model more accurate, Cattaneo[13] and Vernotte[14] established a thermal wave model (C-V model) with a finite thermal propagation velocity by involving the heat flux relaxation time τ. Ulteriorly, Tzou[15] promoted this heat transfer model by introducing another relaxation time τ for the phase lag of temperature gradient. The improved model is called dual-phase-lag (DPL) bio-heat conduction model which allows either the temperature gradient to precede heat flux vector or the heat flux vector to precede temperature gradient. Researchers[16,17] investigated the heat transfer problems in spherical coordinate system with DPL model and numerical method employed. The influences of the non-Fourier effects induced by the phase lags of the heat flux and the temperature gradient on the thermal behavior were considered. Afrin[18] explored the non-equilibrium heat transfer in living biological tissues by considering the arterial and venous bloods. They found that the phase lags of the heat flux and the temperature gradient depend on the properties of artery, vein and tissue, blood perfusion rate and convective heat transfer rate. Nobrega[19] studied the heat conduction in the skin subjected to the pulse laser heating and fluid cooling. The analytical solutions for the one-dimensional skin model were presented with three different heat conduction models employed. Jaunich[20] analyzed the temperature distributions in the skin which is subjected to a short pulse laser. The results from the hyperbolic model and the parabolic model were compared.

It is of great significance to explore the mechanism of thermal damage to make sure of the reliability and safety of thermal treatments. For this purpose, Henriques[21] proposed an expression for the protein denaturation process based on the first order approximation of the Arrhenius equation, which has been widely used in the burn degree prediction. Lu and co-authors[22,23] developed a computational approach to examine the heat transfer process, heat-induced mechanical response, as well as the associated pain level with numerical simulation and the finite difference method employed. Liu[24] analyzed the thermal response for estimating thermal damage in laser irradiated biological tissue with DPL model employed. Afrin[25] investigated the thermal damage in the living biological tissue induced by laser irradiation with a generalized DPL model employed. The results derived from the generalized and classical DPL model were compared and the influences of the laser parameters were discussed. Wang[26] investigated the thermal response of the skin which is subjected to the pulse boundary heat flux with Laplace transform technique employed. Venugopalan[27] explored the thermodynamic response of biological tissue by experimental and analytical methods.

The complication of the transient heat conduction equations with non-Fourier effects under consideration makes the theoretical results very difficult to derive. Green’s function method can be beneficial to solve the complex differential equations of mathematical physics. Zur[28] investigated the free vibrations of functionally graded circular plates with elastic supports on the basis of the classical plate theory. Green’s function method was employed to solve the boundary value problem. Chao[29] presented the thermoelastic behavior of two circular inclusions in an infinite plane which is subjected to a point heat source. The solution obtained in the study can be treated as the Green’s functions for the crack problem associated with two circular inclusions. Seremet[30,31] inspected the steady-state thermoelastic problems for domains described in cylindrical coordinate system and Cartesian coordinate system. A new method was proposed to derive the Green’s function for incompressible Lamé equations and the Green’s function method was used to solve the governing equations. Zhao[32] investigated the coupled thermoelastic vibration of the Euler–Bernoulli beam with cracks which was subjected to heat flux. The explicit expressions of the temperature and displacement response in the beam were obtained with the usage of Green’s function method. Ma[33] has solved some heat transfer problems by using Green’s function method, such as the three-dimensional living biological tissue subjected to a scanning laser beam, the DPL heat conduction in a finite medium subjected to a moving feat source, and the bi-layered circular plate irradiated by laser pulse.

In clinical medicine, the continuous laser or pulsed laser is selected for different disease. In addition, it is reported that both single pulse laser and repetitive pulse laser can be used[32]. At present, the thermal response of human tissue under continuous laser or single pulse laser has been investigated widely. However, there is few reports about the response induced by repetitive laser. In this case, the thermal conduction equations should be solved respectively in different time segmentations because of the discontinuity of the thermal source, which makes it difficult to predict the temperature distribution during the thermal ablation. To address this issue, the present study will focus on the heat transfer conduction procedure in a three-dimensional model in vivo skin which is subjected to a repetitive pulse laser. The governing equations will be established on the basis of DPL model. Exact solution of temperature distribution will be obtained by employing the Green’s function method, which is meaningful for the further exploration on the physical mechanism of the thermal treatment. The thermal damage will be evaluated and the influences of non-Fourier effects will be discussed.

Mathematical models

The heat conduction model employed in the present study is the DPL model, which is modified from the Fourier’s law with two relaxation times under consideration. This heat conduction model can be expressed as[34]:where, τ and τ represent the phase lags of the heat flux and the temperature gradient, q the heat flux, k the heat conductivity, T the temperature, and t the spatial and temporal coordinates, respectively.

The thermal conduction properties of the skin are quite different from that in the classic homogeneous materials. The thermal response will be influenced by the convection between blood and the skin, blood perfusion in vascular beds and metabolic heat generation. Pennes[3] considered the influences as an average distribution in the biological tissue and established a bio-heat conduction model (Pennes model) which has been widely used for decades of years. With the combination of the Pennes model and the DPL model[15], proposed a modified bio-heat transfer model, which can be written as:where, ρ and C are the density and specific heat of the skin. ρ, C, w represent the density, specific heat and perfusion of blood, respectively. T is the temperature of arterial blood, Q for the metabolic heat generation and Q for the external heat source.

In the present study, the thermal transfer process in a cuboid skin model which is subjected to laser pulses is investigated. As is shown in Fig. 1, the length, width and height of the model are l, b and h, respectively. The rectangular laser pulses irradiate the cuboid at the center of the top surface. The governing equation on the basis of DPL bio-heat conduction model can be expressed as:where, θ  = T – T is the temperature rise in the skin, T for the core temperature of human body. Q represents the power density of the laser source, which can be written as:φ(z), Q(x), Q(y) and Q(t) are expressions for the spatial and the temporal distributions of the laser energy, which are given bywhere, µ is the absorption coefficient of the skin, H(*) for the Heaviside function, R for the optical reflectivity, I for the power density of the input laser energy, R for the geometric dimension of the laser spot. to( and te( are the onset moment and the ending moment of the ith pulse and n is the number of the pulse.

Figure 1

Illustration of the 3D skin model.

The model is considered as a cuboid domain in human body. So temperatures of the side surfaces and the bottom surface are assumed equal to the core temperature T. The top surface is considered as the skin so the natural convection between the top surface and the air is under consideration. The boundary conditions can be expressed as:where, f = − h (T − T). T is the ambient temperature and h is the convective heat transfer coefficient.

The initial conditions can be expressed as:

The laser pulse which irradiates the skin is treated as the volume heat generation in the substrate. Figure 2 shows the spatial distribution of the input energy power. The energy input is in the exponential decay distribution along depth direction. On the top surface of the skin, the heating energy evenly distributes in a square area. In the time domain, the laser energy is distributed in the form of repeated multi-pulse, as is shown in Fig. 3. It is assumed that the period of each laser pulse is t with the duration time being t. It is clear that t = to( − te(. The duty ratio is defined as r = t / t, which is an important parameter in determining the thermal response of the skin.

Figure 2

Spatial distribution of the heat source power.

Figure 3

Temporal distribution of the heat source power.

It is universally accepted that the denaturation of protein can be used to measure the degree of thermal damage. The denaturation rate defined by[21] is given as:where, A represents the frequency factor, R for the universal gas constant and Eα for the activation energy of the denaturation reaction. The thermal damage is determined by temperature and exposure duration, which can be evaluated as[21]:

Exact solution for governing equations

Solution in terms of Green’s function

In 1985, Frankel[35] presented the Green’s function approach for hyperbolic heat conduction in a one-dimensional medium. In the following, we will extend this procedure into DPL model in a three-dimensional finite medium and develop the general solution to Eqs. (3)–(10). The results will be used to solve the present heat conduction problem.

Denote an operator L as following:where, , .

Then Eq. (3) can be expressed as:where,

According to the procedure presented by[35], the following integration is defined:where is the corresponding Green’s function for Eq. (3), which is written to represent the “effect/cause” relationship. Here we denote and to show the point. The introduction of ε > 0 is to invoke causality at a later time in the analysis[35]. The operator L’ in Eq. (14) with respect to the cause variables x’, y’, z’ and τ can be written as:where, .

Integration of the right part of Eq. (16) by parts yields the following expression:where,where L’* is the formal adjoint operator of the operator L’, which can be written as[35]:

The Green’s function for Eq. (3) is chosen to satisfy the following auxiliary problem:

And the following homogeneous boundary conditions and additional requirement should be satisfied:

Equation (27) represents the causality principle, which is merely a statement that no effect can be experienced prior to a cause . Substituting Eqs. (19)–(23) into Eq. (18) and taking the limit as in Eqs. (19)–(23) yields

So the temperature distribution of the skin model in the present study can be expressed with the usage of Green’s function as following:

Determination of the Green’s function

The Green’s function can be derived from Eqs. (25)–(27). The associated equation can be written as:

The mode functions satisfying the boundary conditions are[36]:where, , , . According to the boundary conditions, ζ can be derived by the following equation:

By operating on Eq. (30) with , and incorporating the homogeneous boundary conditions of Eq. (26), one can obtain:where,

The inversion formula iswhere , , .

With Eq. (27) under consideration, the solution for Eq. (33) can be written as:where, , .

Substitution of Eq. (37) into Eq. (36) yields the Green’s function for the governing equations as:

Solution for the nonhomogeneous equation

With the boundary conditions and the initial conditions under consideration, the solution for Eq. (3) can be expressed as:

Substitute Eqs. (4)–(8), (15) and (38) into Eq. (39), and the solution for the nonhomogeneous Eq. (3) can be derived as:where,

Numerical results and discussions

The parameters of the heat source and the skin in the present study are shown in Tables 1 and 2, respectively[19,37,38]. The core temperature of human body is usually set as 37 °C, and the temperature of arterial blood and the initial temperature of the vivo skin can also be set as T = 37 °C and T = 37 °C. The ambient temperature is described as T = 20 °C .

Table 1

Parameters of laser beam.

Parameters	Values
I0, parameter of laser power(W/m3)	3 × 104
Ra, energy reflectivity coefficient	0.024
R0, size parameter of laser spot (m)	2 × 10–3

Table 2

Thermal parameters of the skin.

Parameters	Values
Skin specific heat (J/kg∙K)	2348
Skin density (kg/m3)	911
Blood specific heat (J/kg∙K)	3770
Blood density (kg/m3)	1060
Thermal conductivity (W/m∙K)	0.235
Metabolic heat Generation (W/m3)	368.1
Length (m)	0.02
Width (m)	0.02
Thickness (m)	0.01
Blood perfusion (s−1)	0.0016

The burn assessment model employed proposed by[21] is used, which is based on the description of protein thermal denaturation rate. And the parameters of the model are shown in Table 3[19]. It is accepted that: represents the first-degre burn, for the second degree-burn (irreversible burn) and for the third-degree burn[15,39].

Table 3

Parameters for thermal damage prediction.

Temperature range (oC)	Eα/R (K)	A, frequency factor (s−1)
T ≤ 55	7.5 × 104	3.1 × 1098
T > 55	3.54 × 104	5.0 × 1045

The exact values of the two phase lags for the skin are still under study. However, it is widely agreed that the values fall to the order of magnitude ~ s[40-42]. So the phase lags will be chosen based on this range.

Degeneration to the Fourier’s law

The DPL heat transfer model can degenerate to the Pennes model by setting τ = τ = 0. The expression for the distribution of the temperature rise based on the Pennes model can be expressed as:where, the is the Green’s function for the Pennes model, which can be written as:

Substituting Eq. (42) into Eq. (41) can yield the following formulas:where,

A finite element simulation is carried out to validate the analytical solution derived in the present study (by ABAQUS-6.14, Dassault Systemes Simulia Corp.). For ease of comparing the numerical and analytical solutions, the value of blood perfusion (w) is set to be zero and the analytical solution is reduced to the Pennes model. The Abaqus software is employed for the finite element analysis. The skin model is meshed into 160,000 elements, which are set to be 8-node linear heat transfer bricks. The grids is densified at the center point of the top surface, as is shown in Fig. 4. With the usage of Fortran, the multi-pulse heat source is applied as a body heat generation with the expression of the heat source given by Eqs. (4)–(8). The simulation is set to be transient and the time period (20 s) is divided into two stages: the heating stage (0 ~ 3 s) and the non-heating stage (3 ~ 20 s). The temperature increment of the two stages are 0.01 s and 0.1 s, respectively.

Figure 4

The FEM model of the skin (created by Abaqus 6.14, https://www.3ds.com/zh/products-services/simulia/products/abaqus/).

Figure 5 shows the comparison of the results derived by the FEM and analytical solution. The heat source involves three pulses (n = 3) and the duty ratio is set to be 0.3 (r = 0.3, t = 1 s and t = 0.3 s). The boundary condition on the top surface of the model is natural convection condition. And the temperature on the other surfaces is equal to the core temperature T. The value of the convective heat transfer coefficient is set to be h = 10 W/(m2 K) with the assumption of natural convection between the skin and the air. It is shown in Fig. 4 that the results derived by FEM and analytical solution are in good agreement.

Figure 5

Comparison between FEM and analytical solutions.

Results based on DPL model

In the following, the thermal response of the skin irradiated by repetitive laser pulses are investigated based on the DPL model.

Figure 6 shows the influences of the heat flux phase lag τ and temperature gradient phase lag τ on the temperature distribution induced by repetitive laser pulses. The pulse number is set to be n = 10 and the duty ratio is r = 0.5 (t = 1 s, t = 0.5 s). It is shown in Fig. 6 that the temperature variation can be divided into two stages: the heating stage and the non-heating stage. The heating stage is the first ten seconds during which the laser pulses irradiate on the skin and the temperature rises in a serrated form[43,44]. The non-heating stage is the latter ten seconds when the heat source is taken off and the temperature decreases due to the heat transfer and the convective heat dissipation. Figure 6 (a) and (b) show the temperature variation with time at the center point of the irradiated region (x = l/2, y = b/2, z = 0). A great τ leads to a tardy heat conduction progress, which means the temperature variation shows up before the heat transfer starts. On the opposite, a greater τ leads to a faster heat balance progress, which means that more heat flux will be transferred from high-temperature region to the low-temperature region before the temperature variation appears. It is shown in Fig. 6a that when τ takes a great value, the temperature rises rapidly during the first three pulse periods and then the growth tendency slows down. During the last three pulse periods, the tendency of temperature growth nearly disappears especially when τ = 10 s. In the non-heating stage, a greater τ leads to a faster cooling rate because of the fluctuation characteristic caused by τ. If τ is great enough, the temperature can be even lower than the initial temperature. Figure 6b demonstrates that the temperature rise induced by every laser pulse tends to be the same and the tendency of temperature growth tends to be linear when τ increasing from 0.1 s to 10 s. The larger the value of τ is, the more smoothly the temperature decreases during the non-heating stage. Figure 6c,d show the peak temperature variations (appearing at the moment t = 10 s) on τ and τ. It can be found that the peak temperature of the skin induced by the multi-pulse heat source increases with the increment of τ and the decrement of τ. The phase lag of heat flux strengthens the influence of external heat source and leads to intense temperature variation by hindering the heat transfer process. On the contrary, the phase lag of the temperature leads to gentle temperature change by promoting the thermal balance process.

Figure 6

Influences of the phase lags on the temperature variation induced by repetitive laser pulses.

Figure 7a shows the influences of the pulse number on the temperature variation of the skin. In the present study, the total input energy is assumed to be constant and not change with the pulse number. It is shown in Fig. 7a that the amplitude of temperature rise is the maximum when n = 1 and it decreases obviously with the increment of n. Figure 7b shows the influences of duty ratio on the temperature variation. The vibration amplitude of temperature increases with the decrement of duty ratio. However, the overall level of the temperature changes little. When r = 1, the multi-pulse heat source regenerate to a constant heat source and the temperature changes smoothly. Figure 7c shows the peak temperature variation with the pulse number. The temperature reaches the maximum value when n = 1 and then reduces rapidly. With the increasing of n, the reduction tendency keeps attenuating and the value of the peak temperature tends to be stable. Figure 7d shows the peak temperature variation with the duty ratio. The increment of the duty ratio leads to the decrement of the peak temperature by decreasing the amplitude of the temperature vibration. The pulse number and the duty ratio do not affect the total input energy of the laser beam, but the temporal distribution of the laser energy is decided by the two parameters.

Figure 7

Influences of the laser parameters on the temperature variation.

During the laser therapy, the skin is subjected to the laser pulses so that the temperature rise and the thermal damage occur in the irradiated region and the tumor cells will be killed. In order to protect the healthy part and kill the tumor efficiently, it is significant to estimate the burn degree in the irradiated skin. In the present study, the Henrique’s model[21], which is based on the protein denaturation rate estimation, is adopted to assess the burn degree. Figure 8 shows the thermal damage variation with time at the center point of the irradiated area (x = l/2, y = b/2, z = 0). The parameters of the heat source are set to be I = 3 × 104 W/m3, r = 0.5, n = 10. The ordinate is set to be logarithmic to show the thermal damage more distinctly. Figure 8a shows the influences of the heat flux phase lag τ. The phase lag of heat flux postpones the thermal balance progress and aggravates the heat accumulation. Consequently, more intense temperature rise will be induced, which has been shown in Fig. 6, and more severe burn is induced in the irradiated region. The value of the burn index Ω increases sharply during the heating stage (the former ten seconds) and then keeps constant after the heat source disappears. It is suggested that the thermal damage is irreversible when Ω reaches 1. Figure 8b shows that the influences of the temperature gradient phase lag τ is opposite to that of τ. A great τ will promote the heat equilibrium process and reduce the heat accumulation in the irradiated region. As a result, the burn degree decreases with the increasing of τ. It is shown in Fig. 8b that in the condition of τ = 1 s, no irreversible burn occurs in the skin when the temperature gradient phase lag increases to τ = 10 s.

Figure 8

Influences of the phase lags on the thermal damage induced by the multi-pulse laser.

The parameters of the heat source, including the pulse number and the duty ratio, show crucial effects on the thermal response of the skin. As a result, the thermal damage in the irradiated region is also under the influences of these parameters. Figure 9a,b show the temporal distribution of the burn at the center point of the irradiated area. As is shown in Fig. 9a, the burn degree decreases with the increment of the pulse number n. This tendency is tremendous especially when n ≤ 5. When n keeps increasing, the burn degree tends to that caused by a constant heat source. Figure 9b shows the influences of the duty ratio on the thermal damage. Under the premise of the constant energy input, the increment of the duty ratio results in the decrement of the burn degree. The obvious difference caused by duty ratio shows up at the end of the heating stage (i.e., the moment t = 10 s). Figure 9c,d show the thermal damage distribution along the depth direction under different pulse numbers and duty ratios, respectively. The burn occurs mainly in the small region close to the top surface of the skin. The thermal damage maximizes at the top surface and decreases with the increasing of the duty ratio and the pulse number. However, the comparison of Fig. 9a,c and b,d shows that the influences on the thermal damage caused by the pulse number are much more remarkable than that caused by the duty ratio. The area of region where is irreversibly burnt increases obviously when the value of n decreases from 5 to 1.

Figure 9

Influences of the heat source parameters on the thermal damage induced by the multi-pulse laser.

According to the Henrique’s model, the thermal damage is regarded to be irreversible when Ω reaches 1. In order to describe the irreversibly burnt region more clearly, the irreversible burn depth at the center point of the irradiated region has been calculated and the influences of the duty ratio and the pulse number are as shown in Fig. 10. Figure 10a,b show the time history of the depth. The irreversible thermal damage appears at the moment t = 7 s, extending downward for about eight seconds and then trends to be stable at the moment t = 15 s. The decrement of the pulse number leads to an increasing burn depth and takes a shorter time for the irreversible burn to appear, especially when n ≤ 5. The decreasing of the duty ratio also benefits the thermal damage. However, the appearing moment of the irreversible thermal damage changes little when the value of r increases from 0.1 to 1. Figure 10c,d show the irreversible burn depth variation with the pulse number and the duty ratio, respectively. The burn depth reduces sharply when the pulse number increases from 1 to 5 and then tends to be stable with the continuous increment of pulse number. That is to say, on the premise of the constant input energy, single-pulse heat source causes the much deeper burn than multi-pulse does. The increasing of the duty ratio also leads to a decrement of the irreversible burn depth. However, the influence of the pulse number is more important than that caused by the duty ratio.

Figure 10

Influences of heat source parameters on the depth of the irreversible burn induced by the multi-pulse laser.

Conclusions

The temperature response and thermal damage of skin were derived when it is subjected to repetitive pulse laser based on the DPL biological heat transfer model and the Henrique’s burn assessment model.

In the present study, the influences of the two phase lags and the heat source parameters on the temperature distribution and the thermal damage in the skin were researched. When τ is great, the temperature rises rapidly in the initial period of the heating stage and the rise tendency of the temperature gradually disappears as time passes. In addition, the increment of τ also leads to the increasing of the maximum temperature and the thermal damage in the irradiated region. On the contrary, the increasing of τ results in the decrement of the maximum temperature and the burn degree.

The duty ratio and the pulse number obviously affect the thermal response of the skin. The quantity of the total input energy is set to be constant. That is to say, the values of n and r do not affect the quantity of the energy input but the energy temporal distribution. The peak temperature maximizes when n = 1 and decreases obviously with the increment of npulse. As a result, the burn degree and the depth of the irreversible thermal damage decreases with the increment of n. The influence of the pulse number is the most obvious when n = 1, keeping decreasing with the increment of n and gradually disappearing when n is great enough (n ≥ 50 in the calculation for example). The duty ratio has similar effects on the skin thermal response with the pulse number does. The increment of r causes the decrement of the peak temperature, the burn degree and the depth of irreversible burn. However, the influences of the pulse number is more obvious than the duty ratio.

The validation of thermal damage parameter seems difficult because the quantitative description and measure is hard. However the temperature result validation is possible.