Thermo-vibrational analyses of skin tissue subjected to laser heating source in thermal therapy

Laser-induced thermal therapy, due to its applications in various clinical treatments, has become an efficient alternative, especially for skin ablation. In this work, the two-dimensional thermomechanical response of skin tissue subjected to different types of thermal loading is investigated. Considering the thermoelastic coupling term, the two-dimensional differential equation of heat conduction in the skin tissue based on the Cattaneo–Vernotte heat conduction law is presented. The two-dimensional differential equation of the tissue displacement coupled with the two-dimensional hyperbolic heat conduction equation in the tissue is solved simultaneously to analyze the thermal and mechanical response of the skin tissue. The existence of mixed complicated boundary conditions makes the problem so complex and intricate. The Galerkin-based reduced-order model has been utilized to solve the two-sided coupled differential equations of vibration and heat transfer in the tissue with accompanying complicated boundary conditions. The effect of various types of heating sources such as thermal shock, single and repetitive pulses, repeating sequence stairs, ramp-type, and harmonic-type heating, on the thermomechanical response of the tissue is investigated. The temperature distribution in the tissue along depth and radial direction is also presented. The transient temperature and displacement response of tissue considering different relaxation times are studied, and the results are discussed in detail.


Mina Ghanbari 1* & Ghader Rezazadeh 2,3
Laser-induced thermal therapy, due to its applications in various clinical treatments, has become an efficient alternative, especially for skin ablation. In this work, the two-dimensional thermomechanical response of skin tissue subjected to different types of thermal loading is investigated. Considering the thermoelastic coupling term, the two-dimensional differential equation of heat conduction in the skin tissue based on the Cattaneo-Vernotte heat conduction law is presented. The two-dimensional differential equation of the tissue displacement coupled with the two-dimensional hyperbolic heat conduction equation in the tissue is solved simultaneously to analyze the thermal and mechanical response of the skin tissue. The existence of mixed complicated boundary conditions makes the problem so complex and intricate. The Galerkin-based reduced-order model has been utilized to solve the two-sided coupled differential equations of vibration and heat transfer in the tissue with accompanying complicated boundary conditions. The effect of various types of heating sources such as thermal shock, single and repetitive pulses, repeating sequence stairs, ramp-type, and harmonictype heating, on the thermomechanical response of the tissue is investigated. The temperature distribution in the tissue along depth and radial direction is also presented. The transient temperature and displacement response of tissue considering different relaxation times are studied, and the results are discussed in detail.
In modern clinical treatments, thermal therapy is considered one of the most efficient existing alternatives. In thermal therapy, electromagnetic (EM) energy, ultrasonic waves, and other devices based on thermal conduction have been used as heating sources 1 . Thermal therapy aims to modify tissue temperature in a targeted region over time to compel a desired biological response. In the majority of designed thermotherapies, thermal therapy is delivered to a target tissue volume with minimal impact on intervening or surrounding tissues. Thermal therapy emphasizes three techniques, namely, hyperthermia and thermal ablation. Hyperthermia can be categorized into long-term low-temperature hyperthermia in which the tissue temperature rises to 40-41 °C for 6-72 h and moderate-temperature hyperthermia in which tissue experiences 42-45 °C for 15-60 min. In thermal ablation, also known as high-temperature hyperthermia, the tissue temperature rises to higher than 50 °C for more than 6 min 1 . Laser-induced thermal therapy (LITT) plays a significant role in many oncology services as an alternative to conventional surgical interventions, especially for patients who are not good candidates for surgery. In other words, LITT is a percutaneous tumor-ablation procedure that delivers therapy due to the interstitially placed high-power lasers in the tumor. It delivers thermal therapy into the tumor cavity under MRI guidance, causing tumor ablation by compelling coagulation necrosis 2,3 . Reported studies show that it is also a very good option for the treatment of pediatric brain tumors 4 and small palpable invasive breast carcinomas 5 .
Skin ablation utilizing energy-based devices has attracted increasing interest in the last few years. Skin ablation is an effective process not only for cosmetic purposes such as resurfacing, treating scars, or antiaging but also for therapeutic applications. In skin ablation, as the energy is divided into fractions, deep dermal penetration of the energy is achieved with minimal effect on the epidermis. It ensures a rapid recovery time compared with traditional ablative lasers 6 .
Several developed ablative and non-ablative laser devices have been provided physicians with a wide palette of treatment options. In some studies, employing ablative lasers for the treatment of photoaged skin has shown efficacious results 7,8 . However, these painful procedures have several side effects, such as pigment alterations, infection, scarring, long-lasting erythema, and significant downtime. Ablative fractionated CO2 or erbium: yag lasers are commonly applied techniques. The literature on non-ablative fractional lasers (NAFLs) and ablative www.nature.com/scientificreports/ fractional lasers (AFLs) was reviewed by Tierney et al. 6 . This review supported the use of NAFL and AFL as an effective and safe treatment for photoaging. They showed that for the treatment of photoaging, fractionated resurfacing has important advantages over ablative laser resurfacing treatments. Bipolar radiofrequency is the other alternative method for the treatment of skin laxity, wrinkling with a low risk of scarring or persistent pigmentation 9 . It has been shown that bipolar fractionated RF treatment is considered a rapid wound healing response in which treatment intervals of at least 14 days should be recommended to allow the performance of the remodeling process. In new technologies, high-intensity focused ultrasound (HIFU) is used to treat UVinduced hyperpigmentation. This method was commonly indicated for skin laxity. Vachiramon et al. reported the efficacy and safety of high-intensity focused ultrasound for UVB-induced hyperpigmentation in human subjects 10 . In the last few years, a new technology based on thermomechanical principles has been developed to present a novel treatment modality. Fractional treatment of aging skin with Tixel was reported by Elman et al. 11 . A novel device of Tixel, which is based on thermomechanical ablation technology, synthesizes temperature control and sophisticated motion. By applying Tixel as a novel technology, average treatment pain, downtime, and erythema clearance could be improved significantly. Therefore, Tixel could be used safely for nonablative and ablative resurfacing, such as scarring, the incidence of bleeding, or postinflammatory hyperpigmentation. Sintov and Hofmann assessed the effect of Tixel on the skin penetrance of three hydrophilic molecular models 12 .
In the present work, no significant damage to the skin tissue or dermal coagulation was reported. The underlying wound-healing processes after skin ablation with thermomechanical ablation were investigated by Kokolakis et al. 13 . It was concluded that the wound-healing process after TMA is much faster, and the recovery time minimizes significantly compared to other ablative techniques. The investigation of coagulation characteristics and ablation of a new CO laser and a high-power Tm:fiber laser were presented to assess their potential application for fractional ablation of the skin 14 . The use of a CO laser showed approximately two times larger coagulation zones than the CO 2 laser. The practicality of using noninvasive methods for objective skin assessment was evaluated following skin rejuvenation treatment 15 . After laser treatment, considerable improvement in facial skin aesthetics was recorded in brown spots, UV spots, and pores after 3 weeks, without significant changes in the tissue at the molecular level, as assessed by micro biopsy. Thermal therapies and physiological studies such as hyperthermia, skin and tumor ablation, cryosurgery, frostbite, skin burns, and body thermal regulation and response to environmental conditions address the temperature rise in living tissues rising from the absorbed heat of the laser plus. Therefore, temperature distribution and the process of heating transfer and in living tissues play an important role in the effectiveness of thermal therapy methods. Presenting exact mathematical modeling of the heat transfer phenomenon is difficult due to the complexity of the nature of heat transfer in living tissue. Simplifications and assumptions must be made to make the problem obedient, in addition to capturing the essential features of the process. As the determination of temperature distribution in blood perfused tissue is very significant, heat transfer in bio tissues has been studied numerously.
Several bioheat equations have been reported for modeling heat transfer in living tissue. The first bioheat equation was reported by Penne based on simplifying assumptions 16 . The assumptions of the presented model concerned four central factors: equilibrium site, blood perfusion, vascular architecture, and blood temperature. All pre-arteriole and post-venule heat transfers between blood and tissue were neglected. The flow of blood in the small capillaries was assumed to be isotropic. The effect of blood flow directionality and the local vascular geometry were also not considered in Penne's model. Despite the shortcomings of Penne's equation, it has enjoyed significant success in many applications, such as thermal stimulation of the whole body, cryosurgery, blood perfusion, hyperthermia therapy, and measurements. Chen and Holmes modified Penne's perfusion term, in which vascular geometry and blood flow directionality were taken into account 17 . Jiji et al. presented a threetemperature model for peripheral tissue where the temperature of arteries, veins, and tissue temperature was defined in the model for analyzing heat transfer 18 . Due to the difficulty of the three-temperature model, the three coupled deep layers were reduced to a single equation for the tissue temperature by applying simplifications 19 . The extracted differential equations for the bioheat equation in the discussed methods were all parabolic equations. They were obtained based on the classical Fourier's law, where heat pulses are supposed to propagate at infinite speed. An extended heat conduction model was introduced by Cattaneo 20 and Vernotte 21 . Applying Cattaneo-Vernotte (C-V) for extracting the differential heat conduction equation causes the partial differential equation to transform from parabolic to hyperbolic. In several types of research, the classical heat conduction equation has been modified to ensure finite speed pulse propagation to be compatible with physiological considerations and physical reality in a transient process 22 . To take into account the microstructure interaction effect in the heat transfer process, Tzou introduced a generalized correlation between heat flux and temperature 23,24 . As the heat conduction model requires two-phase lags of the temperature gradient and heat flux, the heat transfer model is referenced as a dual-phase lag (DPL) model. In this model, the time delay parameter was a new indicator of bioheat efficiency in living tissue.
Numerous studies have been reported on the use of the modified heat conduction equation in analyzing the temperature distribution in tissue. Analytical estimation of the temperature and thermal damage in living tissue based on a hyperbolic bioheat model was investigated by Alzahrani and Abbas 25 . The governing partial differential equation subject to laser irradiation was solved analytically in the Laplace domain. Comparing the analytical results with the existing experimental results showed the effectiveness of the mathematical model for biological heat transfer. In another work, the non-Fourier effect of laser-mediated thermal behaviors in bio tissues was investigated 26 . Optical transmission and energy deposition were obtained by applying the Monte Carlo method. The dual-phase-lag model was applied and solved by employing a three-level finite difference method. Transient heating within skin tissue due to time-dependent thermal therapy was presented in the context of a memorydependent heat transport law 27  www.nature.com/scientificreports/ heat source velocity and the memory-dependent derivative on the thermal injuries and temperature of skin tissues were precisely investigated. Theoretical investigation of the temperature distribution in three-dimensional biological skin tissue was reported, in which the skin tissue was irradiated by multifiber lasers 28 . The results showed that the arrangement layout, spot size, and interval distance of the laser beams affected the obtained irradiated zone. Several types of research have been reported on the dual-phase-lag bioheat transfer model [29][30][31] .
In one of the studies, to describe the interaction of a multipulse heat source and the skin, the dual-phase lag (DPL) bioheat transfer model and Henrique's burn assessment model were employed 32 . The influences of the biological parameters on the temperature distribution were discussed in detail. The thermomechanical behavior taking place in instantaneously heated skin tissue via an analytical approach was explored 33 . The generalized thermoelastic model involving a dual-phase-lag model of bioheat transfer was applied in a multilayer skin structure. Variable mechanical and thermal properties with spatial location and temperature were presented. Ezzat presented analytical thermomechanical responses of viscoelastic skin tissue 34 . The influences of variable thermal conductivity and volume material properties on the heat transfer of bioheat and the mechanical heat-induced response in a human skin plane were examined. Although several works have been presented for modeling heat transfer in bio tissues, the coupling thermoelastic term resulting from stretching or contracting of an elastic living tissue has not been considered in the heat transfer equation. Moreover, the vibration of elastic tissue in exposure to a thermal source, which is important, especially in giving thermal shock to the tissue, is rarely studied 35,36 . In this work, a two-dimensional heat conduction equation of the skin tissue coupled with the two-dimensional equation governing displacement of the tissue was presented based on a modified Fourier law. Very complicated boundary conditions appear in this problem that should be satisfied. The two-side coupled equations of heat transfer and tissue movement are solved simultaneously to analyze the transient temperature and vibrating response of the tissue. The effect of several types of laser heat sources on the response of the tissue is investigated in more detail. The effect of thermoelastic terms on the temperature response of the tissue is also studied. Figure 1 shows the biological skin tissue exposed to laser heating in cylindrical coordinates. When biological tissue receives laser irradiation, it will experience localized heating, optical energy absorption, and thermal expansion. By supposing an axisymmetric temperature distribution in the skin tissue, a two-dimensional distribution is considered along with the r-and z-directions in this study. Tissues can also vibrate along with the radial (r) and the skin depth directions (z) by absorbing optical energy. The two-sided coupled governing equations of the two-dimensional heat conduction and two-dimensional vibration of the elastic media in cylindrical coordinates are extracted. They are solved simultaneously to obtain the temperature response of the bio tissue exposed to laser heat. It is assumed that the radius of the laser beam is so small that the laser beam is assumed to be a localized heating source at r = 0.

Governing equations
According to the anatomical structure of skin tissue, it is generally considered as a multilayered structure consist of epidermis, dermis, fat layers. Different layers of skin tissue have different nonlinear thermomechanical behaviors. Moreover, the physical properties of each layer are also age-dependent and can vary in different ages of humans 37,38 . Few types of research have been reported the nonlinear elastic behavior to identify the mechanical parameters of human skin. Delalleau et al. defined an identification method related to a new numerical model to account for the nonlinear behavior of human skin in vivo 39 . In that work, a relationship between mechanical and physiological aspects was first presented. The orientation of the dermis collagen fibers in the direction of the stress during an experiment made changes in the elasticity of the skin and showed that skin tissue becomes stiffer. This leads to a three-phase behavior law and validates the nonlinear mechanical approach. In another research Finite element modeling of human skin was introduced using an anisotropic, nonlinear elastic constitutive model 40 . It was shown that the nonlinear mechanical stress-strain response of skin is originated from the collagen network in the skin. It was hypothesized that the force-stretch response of collagen is governed by the entropic of long-chain molecules. A constitutive model derived from the statistical mechanics of long-chain molecules, corresponding to the " collagen network in the skin", demonstrated the mechanical response of the www.nature.com/scientificreports/ skin. In this work, for simplicity, a linear single-layer skin model is used where the skin is treated as a single homogeneous layer with constant thermal and physical properties.

Two-dimensional displacement equation of the tissue. The balanced equation of an elastic media
is expressed as 41-43 : where ρ is the tissue density, b i is the body force per unit volume of the tissue, S ij refers to stress components and U i denotes the tissue displacement components. The balance equation in cylindrical coordinates in the r-and z-directions is expressed as [41][42][43] : where u and v and w are displacement components of tissue displacement in the r-, z-, and θ-directions. The constitutive equation of an elastic media in the presence of the temperature distribution is defined as: In Eq. (5), δ ij is the Kronecker delta, and T is the temperature. and µ are the Lame constants related to the modulus of elasticity E and Poisson's ratio ν as 41 : where, ε ij = 1 2 u i,j + u j,i is the strain tensor. In relation (5), B is defined as: where α T is the thermal expansion, E is the elasticity modulus of the tissue and υ is Poisson's ratio. The components of the strain tensor in cylindrical coordinates are expressed as 41 : Substituting strain components (8) into constitutive equation (5), the components of the stress tensor in cylindrical coordinates can be expressed as: T 0 is the reference temperature.
Substituting relation (9) into balance equations (2) and (3), the two-dimensional vibration of tissue in the r-and z-directions takes the following form: www.nature.com/scientificreports/ Two-dimensional hyperbolic heat transfer equation. Due to the importance of temperature distribution in bio tissue in many medical therapies, such as brain tumor invasive breast carcinomas, hyperthermia, frostbite, and skin burns, predicting the accurate temperature response of bio tissue to laser irradiation seems to be very important. Due to the complexity of the nature of heat transfer in living tissue, presenting a mathematical model for heat transfer needs logical simplifications. Several simplifications and assumptions must be made to make the problem tractable while capturing the essential features of the process. The local form of the energy balance equation in index notation is as follows: where σ is the stress, ε is the strain, u is the specific internal energy per unit volume, Q is known as the heat supply and denotes the rate at which heat per unit volume is produced by internal sources, and V is the volume. The rate at which heat is conducted into the body per unit area per unit time across the element of the surface is called the heat flux vector and is represented by the symbol q i . Fourier's law states that the time rate of heat transfer through a material is proportional to the negative gradient of temperature and to the area at right angles through which the heat flows. It specifies a linear relationship between the temperature gradient and heat flux and is expressed as: In Eq. (13), T is temperature and k is the thermal conductivity W m k . In Fourier's law, the diffusion of heat gives rise to infinite speeds of heat propagation, which is incompatible with physical reality and physiological considerations in a transient process. Rapid thermal energy deposition is indicated in the medium, i.e., the occurrence of any local temperature disturbance leads to an immediate perturbation in temperature at each point in the medium. Applying Fourier's law in the energy equation results in Pennes' bioheat equation, which is commonly used for modeling heat transfer in biological systems. Much attention has been assigned to developing the classical heat conduction equation to ensure finite-speed pulse propagation. The generalized Fourier's law is called the Cattaneo-Vernotte heat conduction law. Utilizing it, the governing partial differential equation is transformed from parabolic to hyperbolic type. The constitutive equation of the Cattaneo-Vernotte model is characterized as follows: τ is the relaxation time, assumed to be a non-negative constant.
Combining equations (12) and (14) results in the following heat conduction equation 46 : In Eq. (15), c is the specific heat capacity of the tissue. The thermoelastic coupling term in the heat conduction equation is so small and usually neglected in practice 44 . It is important in many applications, such as thermal shocks 45,46 , ultrafast laser heating in the thermal processing of materials 47,48 , modeling of the vibration of resonant microelectromechanical systems (MEMS) [49][50][51][52] , and dynamic crack propagation 53 . Therefore, the thermoelastic effect should be taken into account in the mentioned application. By considering the thermoelastic effect in Eq. (15), it is modified to 46 : Q is the volumetric heat generated in the tissue by blood perfusion, metabolism, and laser pulses as: In Eq. (17), Q L is the volumetric heat generated from a laser heat source 25 : . I 0 is the laser intensity, µ a is the absorption coefficient, and e 1 , e 2 , O 1 and O 2 are functions of diffuse reflectance (R d ) and are presented by Gardner et al. 25 . f 2 (t) is a time-dependent function and can be considered in any form as a repetitive pulse function, step function, harmonic function, and so on. δ is the penetration depth and is characterized as follows 25 : g is the anisotropy factor, and µ s is the scattering coefficient. (17) is the volumetric heat generated by blood perfusion and is defined as: www.nature.com/scientificreports/ where ρ b and c b are the density and the specific heat capacity of the blood, ω b is the blood perfusion, and T b is the artery temperature.Q M is the heat generated by the metabolic process due to various physiological processes occurring in the rest of the body and is neglected in this study. Substituting relations (8) and (17) in relation (16), the two-dimensional hyperbolic heat conduction equation is expressed as: By considering the following nondimensional parameters: where α = k/ρc m 2 /s is the heat diffusivity of the tissue. Substituting the parameters (22) into Eqs. (10)(11), and (21) yields: with the following coefficients: In Eq. (25) the source term is expressed as: Nondimensional initial conditions for the displacement components of skin tissue and temperature are expressed as: The nondimensional boundary conditions for u and v along r direction are as follows: The nondimensional boundary conditions for u and v along z direction are as follows:

Numerical method
In this study, the Galerkin-based reduced-order model was utilized to solve the coupled two-dimensional partial differential equations of hyperbolic heat transfer and vibrating bio-tissue equations. Utilizing the Galerkin method, a continuous operator problem, such as a boundary value differential equation is converted to a discrete problem. It is done by applying linear constraints determined by finite sets of basic functions. This method is the equivalent of applying the method of variation of parameters to function space by converting the equation to a weak formulation. It is based on the weak formulation of an equation and limits the possible solutions to a smaller space than the original one, which can be solved more easily. A linear combination of a set of prescribed basis or shape functions is considered for the unknown function. Satisfying the boundary conditions of the problem is dependent on selecting proper shape functions. The accuracy of the solutions depends significantly on the number and type of shape functions. In this work, complex boundary conditions, especially those related to stresses on the surface of the tissue, appear. To satisfy the complicated boundary conditions, the following approximate solutions are considered for Eqs. (23)- (25).  www.nature.com/scientificreports/ Applying the Galerkin integrals as: gi K (8) hj + A i K (9) hj

Numerical results
For numerical calculation, values of the heat and physical properties of the skin tissue, blood, and the parameters of the laser heat source are listed in Table 1.
For the modeling and its solution to be more understandable, a clear flow chart of the modeling and the numerical method is presented in Fig. 2. The coupling terms in equations are shown in red font in the chart.   Single-pulse heating source. Considering the single pulse heating source shown in Fig. 3a, the temperature response of the skin tissue is presented and shown in Fig. 3b.
The results indicate that in a laser single pulse heating source with high step time, the tissue can reach higher temperatures than a heating source with low step time. It is also observed that the step time does not affect the gradient of temperature in the heating process but increases the temperature gradient in the cooling process. The displacement response of tissue in the depth direction (z-direction) is shown in Fig. 4a, and in Fig. 4b, it is compared to the displacement response of tissue in the radial direction (r-direction).
A higher step time increases the maximum displacement of the tissue and reaches the maximum temperature. It also increases the displacement gradient of tissue in the cooling process, while it does not affect the gradient of displacement in the heating process. Figure 5 shows the temperature of the skin tissue versus time considering different values of relaxation time. It was mentioned that in the generalized Fourier law, an inertial term is considered for the heat transfer phenomenon, which is concerned with the relaxation time. An increase in the relaxation time leads to an increase in the inertial effects and therefore causes the reached maximum temperature to decrease.
Repetitive pulse heating source. Considering the repetitive pulse heating source shown in Fig. 6a, the displacement response of the skin tissue is shown in Fig. 6b. The results indicate that increasing relaxation time decreases the temperature gradient at the beginning of the heating process. In other words, relaxation time decreases the first maximum overshoot and increases the pick time in the response of the skin tissue. Relaxation time also reduces the amplitude of temperature vibrations of the skin in the whole heating process. In this work, the heat conduction phenomenon in the skin tissue is modeled based on the non-Furies conduction equation   www.nature.com/scientificreports/ which is attributed to the Cattaneo-Vernotte model. In the presented model, τ is the relaxation time. It is a function of the material. It is defined as the time necessary for storage of the thermal energy required for the propagative transfer to an element in thermal contact. In other words, the term τ ∂q ∂t in the non-Fourie's heat conduction equation behaves as an inertial term in the conduction phenomenon which restricts the infinite propagation speed in the bio tissue. So, by increasing the relaxation time, the effect of the thermal inertia term increases. Therefore, the maximum temperature that the skin can be reached during thermal radiation decreases. Inertia term also creates a delay in the thermal and vibrating response of the tissue that leads to the increment of peak time in the time-domain thermomechanical response of the tissue.
The effect of the period and width of the repetitive pulses on the thermomechanical response of the tissue was investigated. The indicated results in Fig. 7a,b show that when the period of the pulses reduces, the amplitude of temperature vibrations decreases; therefore, the temperature response of tissue will be similar to the response subjected to the step heating source. Moreover, by increasing the width of laser pulses, the tissue can reach high temperatures during the heating process.
Repeating sequence stair. Different types of repeating sequence stair heating sources are shown in Fig. 8a,b,c in which the equal energy is delivered to the tissue in all cases. In Fig. 8a, the sequence stairs are applied in the cooling stage, while in Fig. 8b, the sequence stairs are applied in the heating stage, and in Fig. 8c, symmetric sequence stairs are applied in the heating and cooling stages of the skin heating. The non-dimensional temperature of the skin tissue versus time is shown in Fig. 8d). When sequence stairs are applied in the heating stage, the tissue can experience a higher maximum temperature than when they are applied in the cooling stage. In other words, heating the skin gradually by sequence stairs and cooling it suddenly, tissue can experience higher temperature rise compared to the case when it is heated suddenly and cooled gradually. When the tissue encounters thermal shock whether, in the heating or cooling process (cases (a), and (b)), it delivers or losses the www.nature.com/scientificreports/ energy suddenly. Therefore, the high percentage of the delivered energy is converted to mechanical energy which leads to the vibration of the tissue rather than the increment of its temperature. However, in case (c) where the tissue is heated and cooled gradually via sequence stairs, most of the delivered energy is used to increase the internal energy of tissue rather than its displacement. Therefore, the maximum temperature, the tissue can be reached via both heating and cooling with sequence stairs (case (c)) is high compared with cases (a), and (b). In Fig. 8e) the time-domain displacement response of the tissue along radial and depth direction is presented. The displacement of tissue in the depth direction is more considerable than the radial direction.
Step heating source. Figure 9 shows the temperature response of the skin tissue, considering a step heating source. At very low relaxation times, the skin tissue does not experience any overshoot in the thermal response. Increasing the relaxation time may lead the tissue to encounter an overshoot during its transient response. The nondimensional displacement of tissue along the z-direction at t = 0.1 for different values of r is shown in Fig. 10a,b. The results indicate that by moving away from the center of the tissue where the pulse heating source is applied, the displacement of tissue in the r-and z-directions decreases and eventually vanishes.
Repetitive ramp heating source. Figure 11a shows the respective ramp heating source with different repeating times, and Fig. 11b) represents the transient response of the tissue. Reducing the repeating time decreases the vibrating amplitude of temperature, therefore causing the skin to experience a low-temperature rise. In Fig. 11c the vibration response of the tissue considering different relaxation times is presented. In lower relaxation time, the tissue encounters with higher overshoot in the time-domain response of it.
Harmonic heating source. The transient thermal response of the skin tissue is studied supposing a harmonic heating source. Figure 12a shows the temperature response of tissue to the harmonic heating source considering different value frequency ratios (ω r ) . The frequency ratio is defined as the ratio of the excitation frequency to the characteristic frequency of the tissue ( α h 2 ). The frequency ratio displaces the time at which the tissue reaches the maximum temperature. Figure 12b Table 1. The blood properties and laser parameters based on which the results are calculated are also listed in Table 1. In Some experimental studies, it has been shown that there are great similarities between human and pig skins, especially in the vascular organization. An experimental study was conducted on laser heating of pig skins considering a single pulse with various exposure duration times and laser powers 54 . The reported experimental data are shown in Fig. 13a). Figure 13b) indicates the numerical results calculated in this work considering a single laser pulse with a laser power of 122 kW m 2 . Neglecting the effect of the thermoelastic coupling term, the maximum temperature the skin reaches several exposure durations is listed in Table 2. A good agreement is observed between the results of the presented numerical solution and those obtained from the experimental study.
The results of the presented coupled thermoelastic model are compared with the results obtained based on the finite element method 55 . Therefore, by considering the initial and boundary conditions, the material properties are the same as those used in Ref. 55 , the maximum temperature of the skin tissue in two different points in the skin tissue is compared with the obtained values in Ref. 55 . The results are presented in Table 3. The obtained maximum temperature of the tissue shows agreement with the results obtained based on FEM. It shows that the present numerical method is effective and accurate.

Conclusion
In the present study, the thermomechanical behavior of the skin tissue exposed to different types of heating sources was presented. 2-D coupled differential equations of the hyperbolic heat transfer and 2-D dynamic displacement of tissue were extracted, considering the thermoelastic coupling term. The mixed accompanying  www.nature.com/scientificreports/ boundary conditions (displacement and force) governing the displacement equation that has made the solution of the equations so complicated were imposed using a heuristic method. A Galerkin-based reduced-order model was utilized to solve the coupled equations with mixed boundary conditions. The transient temperature and displacement response of the tissue were analyzed under different types of thermal heating loading. It was shown that when a single pulse heating source with a high step time is applied on the skin, the tissue experiences a higher temperature rise. No changes in the temperature gradient were observed in the heating stage of the process. However, the temperature gradient in the cooling process was considerable. In the case of repetitive pulses, a decrease in the period of the pulses led to a reduction in the amplitude of the temperature vibrations. It was observed that when a sequence stair is applied in the heating stage, the tissue can reach a higher maximum temperature compared with the case when it is applied in the cooling stage. Increasing the relaxation time in the