Determination of Laser Parameters in Thermomechanical Treatment of Skin Based on Response Surface Methodology

Abstract

An investigation was conducted to examine the photothermal and thermomechanical effects of short-pulse laser irradiation on normal tissues. This study analyzed the impact of short-pulse laser radiation on the heat-affected region within tissues, taking into consideration a set of laser variables, namely wavelength, intensity, beam size, and exposure time. The beam size ranged between 0.5 and 3 mm, and the intensity of the laser radiation ranged from 1 to 5 W/mm² at wavelengths of 532 and 800 nm. A three-layered, three-dimensional model was implemented and studied in a polar coordinate system (r = 10 mm, z = 12 mm) in COMSOL Multiphysics (version 5.4, COMSOL Inc., Stockholm, Sweden) to perform numerical simulations. The Pennes bioheat transfer model, Beer-Lambert, and Hooke’s law are integrated to simulate the coupled biophysics problem. Temperature and stress distributions resulting from laser radiation were produced and analyzed. The accuracy of the developed model was qualitatively verified by comparing temperature and mechanical variations following the variations of laser parameters with relevant studies. The results of Box-Behnken analysis showed that beam size (S) had no significant impact on the response variables, with p-values exceeding 0.05. Temperature (T_max) demonstrates sensitivity to both beam intensity (I) and exposure time (T), jointly contributing to 89.6% of the observed variation. Conversely, while beam size (S) has no significant effect on stress value (Smax), wavelength (W), beam intensity (I), and exposure time (T) collectively account for 71.6% of the observed variation in Smax. It is recommended to use this model to obtain the optimal values of the laser treatment corresponding to tissue with specified dimensions and properties.

1. Introduction

The application of lasers in the medical field has undergone substantial expansion, becoming integral components within medical systems and surgical procedures. Laser technology finds diverse applications across various medical domains, encompassing cancer diagnosis [1], cancer therapy [2], dermatology [3], ophthalmology (e.g., Laser Assisted In Situ Keratomileusis laser coagulation, and optical tomography) [4], and prostatectomy. Moreover, lasers play pivotal roles in cosmetic procedures such as unwanted hair removal and tattoo elimination [5]. Furthermore, lasers have demonstrated remarkable efficacy in the management of oral mucosal lesions, offering precise treatment with improved patient outcomes [6].

Ultra-short laser pulse durations, ranging from femtoseconds to a few picoseconds, offer significant advantages in high-quality and precise material processing. These durations primarily interact with electrons, minimizing heat conduction during the interaction with materials [7]. Consequently, material ablation occurs within a well-defined area with minimal mechanical and thermal damage to the target. In contrast, longer pulse durations, such as nanoseconds, lead to continuous heating of the target material, spreading laser pulse energy through heat conduction beyond the intended spot size. This can result in the boiling and evaporation of the irradiated material, creating an uncontrollable melt layer and potentially causing imprecise machining or marking [8,9].

Regarding interactions with living tissues, pulsed lasers induce both photothermal and photomechanical phenomena, posing risks of damage to affected tissues and surrounding healthy tissues. These challenges underscore the ongoing need for research aimed at enhancing the effectiveness of lasers in medical treatments [7,8]. Short-pulse lasers, a recent advancement in solid-state laser technology, operate at near-infrared wavelengths with pulse durations of less than 1 ps, offering promising avenues for improved medical applications [9].

The 532 nm and 800 nm lasers enjoy wide applications across diverse medical fields due to their specific optical properties and interactions with biological tissues [6]. In dermatology and skin treatments, the 532 nm laser finds common use in addressing vascular lesions such as port wine stains and facial telangiectasia due to its effective absorption by hemoglobin, resulting in vascular coagulation [10]. Conversely, the 800 nm laser is applied in hair removal and skin rejuvenation treatments, targeting melanin in hair follicles and pigmented lesions [11]. Notably, in ophthalmology, the 532 nm laser is preferred for photocoagulation procedures in diabetic retinopathy and other retinal diseases due to its efficient absorption by hemoglobin [12], while the 800 nm laser is utilized in photodynamic therapy for conditions such as age-related macular degeneration [13]. Specifically, both the 532 nm and 800 nm lasers play instrumental roles in laser surgery [14,15], with the former being effective for tissue ablation such as laser lithotripsy for urological conditions [16], while the latter is suitable for soft tissue surgeries, including cutting, vaporization, and coagulation, owing to its deeper tissue penetration and lower water absorption compared to shorter wavelengths [15].

Documented studies have compared continuous wave and pulsed lasers in hyperthermia treatment. Matthewson et al. [17] conducted a study comparing the effects of continuous (Model 212, CVI, Albuquerque, NM, USA) Nd-YAG laser (1064 nm) and pulsed laser (Model 2000, JK Lasers, Rugby, UK), 100 microseconds of duration, 10–40 Hz frequency) with an average power of 1 W on laser-induced hyperthermia concluding no difference in tissue damage regardless of pulsing rate in low-energy excitations. Cios et al. [18] reviewed the clinical and molecular effects of light irradiation across various wavelengths, including UV, blue, green, red, and IR, on dermal cells. Szczepanik-Kułak et al. [19] discussed morphea treatment using pulsed dye lasers (585 and 595 nm), fractional lasers (CO₂ and Er:YAG), Excimer lasers (308 nm), and Nd:YAG (1064 nm) lasers. Kuehlmann et al. [20] provided an updated overview of laser treatments for skin scars and burn wounds, covering indications, modes of operation, effectiveness, and side effects. Kaushik [21] outlined guidelines for fractional laser use in individuals with darker skin tones, suggesting it as a preferred treatment for various dermatological conditions (Fitzpatrick skin phototypes IV to VI). Finally, Eze and Kumar [22] introduced ablative, non-ablative, and fractional photothermolysis procedures for skin rejuvenation.

Shafirstein et al. [23] employed a finite element analysis (FEA) using the COMSOL package (Version 5.4) to evaluate the impact of transitional photodynamic degradation on wine stain laser treatment. The study involved the utilization of laser pulses with a frequency of up to three per millisecond to determine the depth of thermal penetration within the tissue layer through the application of diffusion approximation theory. Darif et al. [24] conducted a simulation to analyze the temperature changes and phases in materials subjected to a KrF-pulsed laser beam on a silicon surface. The pulse duration was 27 nanoseconds, and the impacts were evaluated at varying time intervals, investigating the impact of changing laser parameters. Jouvard et al. [25] and their associates performed a study that involved modeling the impact of laser irradiation on a metal surface. They utilized an Nd:YAG laser pulse with a duration of 180 ns, which resulted in the melting of the metal surface. The process of solidification was then monitored after the laser was no longer in use. Ganguly et al. [26] emphasized the influence of thermomechanical dynamics associated with laser therapy by utilizing a finite element model. The study utilized a three-layered model of the skin that was exposed to a focused Nd:YAG infrared laser beam. The results gained provided an understanding of how laser parameters, such as repetition rate, laser power, and pulse width, impact the ablation process. Liu and Chen [27] used bioheat transfer models to simulate the thermal behavior in the tumor and estimate the power dissipation requirement for the specified conditions during magnetic hyperthermia. Thermal damage was also evaluated based on the Arrhenius equation to build a modified thermal damage model with regeneration terms, and the equivalent thermal dose equation. Singh et al. [28] introduced a model of magnetic nanoparticle-assisted thermal therapy for minimally invasive, non-surgical cancer treatment. Pennes bio-heat and Arrhenius equations were implemented in COMSOL to model heat transfer and evaluate thermal damage in healthy and tumorous models. Kumar et al. [29] proposed a mathematical heat transfer model to predict the control temperature profile at the targeted position. This technique used a dual phase-lag (DPL) model of heat transfer in multilayer tissues with a modified Gaussian distribution heat source subjected to the most generalized boundary condition and interface at the adjacent layers. Wongchadakul et al. [30] suggested a thermomechanical model to simulate laser heating in shrinking tissue, exploring the impact of wavelength, laser irradiation intensity, and irradiation beam area on the process. Higher laser irradiation intensity resulted in increased temperature elevation and tissue shrinkage, with notable differences observed across varying wavelengths.

Skin tumors represent a significant health concern globally, necessitating efficient and minimally invasive treatment modalities. However, existing approaches often pose challenges due to their limited precision and potential damage to surrounding healthy tissue. By delving into the optimization of laser parameters and assessing the resultant thermal effects, our study seeks to address these critical gaps in current methodologies [31]. The incorporation of Box-Behnken optimization analysis allows us to discern optimal parameter values that not only ensure effective tumor ablation but also mitigate thermomechanical impacts on surrounding tissue [32,33]. Through this research, we aim to enhance our understanding of the intricate dynamics underlying skin-laser interactions, ultimately paving the way for safer and more efficacious clinical practices.

2. Materials and Methods

2.1. Model Design

This research considers particularly the characteristics of thermal transfer and mechanical distortion that occur within the layers of the skin during exposure to short-pulsed laser radiation (pulse duration less than 2.2 femtoseconds) under different treatment conditions (intensity 1–5 W/mm², beam size 0.5–3 mm, wavelengths: 532 and 800 nm, and exposure time 30–600 s). A finite element-based numerical simulation has been applied by COMSOL Multiphysics (version 5.4, COMSOL Inc., Stockholm, Sweden) to model the thermal changes and deformations of skin layers.

In this study, a Gaussian nature of the laser beam is assumed, and the Beer-Lambert law was used to describe the attenuation of the laser in the skin layers. The mechanical deformations of skin layers have been studied using equilibrium equations. A two-dimensional and three-layered skin model was used to determine thermal diffusion and mechanical deformations within the target skin for laser thermal treatment in a polar coordinate system (r = 10 mm, z = 12 mm). Using this model, it was possible to evaluate thermal deformations based on the thermal properties of the skin. It took an average of 5 min to compute the solution to the biophysics problem. The adopted model and properties were established based on the findings of previous studies [30,34].

2.2. Heat Transfer Simulation

A mathematical model for the heterogeneous skin tissue domain has been developed and solved after setting the boundary and the initial conditions. This was performed with the modeling of mechanical deformation, to study the interaction between the skin and the laser. The thermal distribution throughout the skin layers was obtained through the application of Pennes’ Bioheat Equation, which outlines the manner in which heat is transmitted throughout the skin and is expressed as follows:

$$\rho C\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+{\rho }_b{C}_b{\omega }_b\left(T_b-T\right)+Q_{met}+Q_{\mathrm{Laser}}$$

(1)

where ρ is tissue density, C is heat capacity, K is tissue thermal conductivity, T is tissue temperature, T_b is blood temperature, ρ_b is blood density, C_b is the specific heat capacity of the blood, and W_b is blood flow rate, Q_met is the metabolic heat source, and Q_laser heat source is associated with the laser source and related to thermal radiation in this analysis.

The boundary conditions were determined as instructed in previous studies [30]. The upper boundary of the skin layers was considered to be in the case of boundary convection:

$$-n\cdot (-k\nabla T)=h_{am}\left(T-T_{am}\right)$$

(2)

where T_am is the ambient temperature (°C), h_am is the convection coefficient of the air (W/m²·K). The outer section of the studied section is assumed by the internal body temperature, except for the upper surface. The absence of interlayer resistance between the epidermal, dermal, and subcutaneous tissue layers was assumed in the study; that is, continuous boundary conditions between these layers were considered.

This research considers particularly the characteristics of thermal transfer and mechanical distortion that occur within the layers of the skin during exposure to laser radiation under several different treatment conditions (intensity 1–5 W/mm², beam size 0.5–3 mm, wavelengths: 532 and 800 nm, and exposure time 30–600 s). A finite element-based numerical simulation has been applied by COMSOL Multiphysics (version 5.4, COMSOL Inc., Stockholm, Sweden) to model the thermal changes and deformations of skin layers.

The laser radiation was applied to the upper surface of the model, and the laser intensity changes along the skin depth were studied based on the Beer-Lambert law:

$$I\left(z\right)=I_0{e}^{\left(-\frac{r^2}{2{\sigma }^2}\right)}\times e^{\left(bz\right)}\times e^{-(a+b)z}$$

(3)

Skin layers are initially at an initial temperature of 36 °C and then suddenly become subject to the influence of laser radiation of constant intensity. The tissue absorption of laser radiation is described by the following equation:

$$Q_{\mathrm{Laser}}=a{I}_0{e}^{\left(-r^2/2{\sigma }^2\right)}\times e^{\left(bz\right)}\times e^{-(a+b)z}$$

(4)

where I determines the intensity of the laser radiation and I₀ is the intensity of the radiation on the surface of the skin, a is the absorption coefficient, b is the scattering coefficient, z is the depth of the skin, and σ is the width of the radiation exposure area.

In this study, isotropic material for skin tissue was assumed, and the model was simplified to a semi-static model, as indicated in previous studies [30,34]. The mechanical solid-state applied to the symmetric model designed using equilibrium equations, stress-strain relations, and strain-strain relations is described as follows:

$$\frac{\partial {\sigma }_{\Pi }}{\partial r}+\frac{\partial {\sigma }_{rz}}{\partial z}+\frac{{\sigma }_{rr}-{\sigma }_{\varphi \varphi }}{r}+F_r=0$$

(5)

$$\frac{\partial {\sigma }_{rz}}{\partial r}+\frac{\partial {\sigma }_{zz}}{\partial z}+\frac{{\sigma }_{rz}}{r}+F_z=0$$

(6)

$${\epsilon }_{rz}={\sigma }_{rz}(1+v)/E$$

(7)

$${\epsilon }_{rr}=\frac{\partial u_r}{\partial r},{\epsilon }_{zz}=\frac{\partial u_z}{\partial z},{\epsilon }_{\varphi \varphi }=\frac{u}{r}$$

(8)

$${\epsilon }_{rz}=\frac{1}{2}\left(\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}\right)$$

(9)

where σ is the stress, ε is the strain, F is the external load (equal to zero in the case of this study), ν is the Poisson’s ratio, E is Young’s modulus, and u is the mean displacement.

$${\epsilon }^{th}={\int }_{T_{ref}}^T\alpha dT$$

(10)

where T_ref is the reference temperature and α is the coefficient of thermal expansion of the tissue.

2.3. FEA Study

Free boundary conditions for all surfaces were also assumed when performing the mechanical deformation analysis. In this study, the mathematical solution of the bioheat equation and the equilibrium equation with their boundary conditions, was performed using the finite element technique using COMSOL Multiphysics. The model was meshed using tetrahedral elements using a mesh generator from COMSOL. The use of this type of element enhances accuracy and resolution, ensuring precise representation of complex geometries and thermal gradients [34]. The convergence sensitivity analysis showed that the mesh sizes are getting smaller, and a value of less than 1% change has appeared in the maximal temperature when the total number of mesh elements reaches 22,708 elements (Figure 1).

2.4. Optimization Study

The Box-Behnken approach is widely employed in industrial and medical domains to determine the optimal solution when multiple parameters affect a specific response variable. This method was utilized to optimize the Direct Laser Interference Patterning process by incorporating a Design of Experiments (DOE) approach.

Four independent variables were selected (wavelength, beam size, laser intensity, and exposure time). After reviewing previous studies related to the use of pulsed lasers in skin therapy, the ranges of the studied variables were established through the selection of the lower and upper bounds. The lower and upper levels of the variables under investigation are presented in

Table 1. The lower and upper levels of the studied variables.

Variables	Lower Level (−1)	Upper Level (+1)
I = Intensity (W/mm2)	1	5
S = Size (mm)	0.5	3
T = Time (s)	30	600
W = Wavelength (nm)	532	800

.

The optimization of the laser thermomechanical effect was performed using Minitab 19 software (Minitab, Inc. located in State College, PA, USA).

3. Results

Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show the effects of wavelength, laser intensity, beam size, and exposure time on thermal and mechanical distributions that appeared in skin layers.

3.1. Model Validation

As explained in the previous paragraphs, high temperatures can lead to thermal damage to skin tissue, causing pain due to thermal expansion. For this reason, it is necessary to control the amount of heat build-up and the thermal limits of the application spot. It mainly depends on the exposure time, the intensity of the laser radiation, and the radius of the laser beam.

Figure 2 shows the results of the circumferential and axial distributions of the apparent temperatures in the layers of the skin after irradiation of the region using two wavelengths (532 and 800 nm) and two intensity values (1 and 2 W/mm²).

Figure 2. The circumferential and axial thermal diffusion caused by the laser radiation at the laser radiation application center (r = 0, z = 0).

Figure 2. The circumferential and axial thermal diffusion caused by the laser radiation at the laser radiation application center (r = 0, z = 0).

On the one hand, it is obvious that an increase in the intensity of the laser radiation causes an increase in temperatures of no less than 20% in both distributions. On the other hand, the temperatures in the diagonal direction seem more stable, with a decrease of no more than 4% between the beginning and the end of the studied section. On the other hand, it is clear from the distribution of temperatures in the axial direction that the laser radiation was more focused on the irradiated focus and that the heat dissipation towards the depth occurs more clearly, as a 30% difference is observed between the center of the application area and the end of the section.

Figure 3 shows the distribution of heat and equivalent stresses at the mentioned wavelengths and intensities (1 and 2 W/mm²) and an exposure time of 30 s. The maximum temperature readings on the skin’s exterior are observed when utilizing a 532 nm wavelength due to the skin’s greater absorption coefficient at this wavelength compared to 800 nm. Upon exposure to laser radiation, the absorbed energy is transformed into heat energy, leading to a rise in temperature. On the other hand, when the exposure time is increased to 600 s, a transfer of the maximum temperature is noticed, corresponding to the wavelength of 800 nanometers, to deeper levels within the skin. This is attributed to the conductive role of the skin tissue, which allows the absorbed laser energy to spread to deeper layers (Figure 4).

The temperature changes with time are studied to accurately track the effect of exposure time. Figure 5 shows the temperature change with time, which allows the evaluation of exposure time. It is noticed that the temperature increases exponentially with time to reach a steady state after approximately 300 s. On the other hand, the increase in intensity by 1 W/mm² caused the temperature to increase by 27%.

Figure 6 shows the effects of laser beam radius and laser intensity, where a three-dimensional representation of the equivalent temperatures and stresses is shown at 600 s exposure time and 532 nm wavelength. Increasing beam radius seems to cause a slight increase in maximum temperature values (1%) while the positions of maximum values were varied. On the other hand, it seems clear that increasing the intensity of the laser beam leads to a noticeable increase in the site of the maximum temperature, which is located on the upper layer of the skin.

Figure 3. Distribution of temperature changes (°C) and equivalent stresses (MPa) at different wavelengths, intensities, and time 30 s.

Figure 3. Distribution of temperature changes (°C) and equivalent stresses (MPa) at different wavelengths, intensities, and time 30 s.

Figure 4. Distribution of temperature changes (°C) and equivalent stresses (MPa) of the studied model at different wavelengths, intensities, and time 600 s.

Figure 4. Distribution of temperature changes (°C) and equivalent stresses (MPa) of the studied model at different wavelengths, intensities, and time 600 s.

Figure 5. Temperature changes with time at intensities of 1 and 2 W/mm² and wavelength of 532 nm at the center of the studied model.

Figure 5. Temperature changes with time at intensities of 1 and 2 W/mm² and wavelength of 532 nm at the center of the studied model.

Figure 6. Distribution of the equivalent temperatures (°C) and equivalent stresses (MPa) across the model layers at a time of 600 s, intensities of 1 and 2 W/mm², and wavelength of 532 nm.

Figure 6. Distribution of the equivalent temperatures (°C) and equivalent stresses (MPa) across the model layers at a time of 600 s, intensities of 1 and 2 W/mm², and wavelength of 532 nm.

It is also clear from the previous Figure 8 and Figure 9 that the wavelengths 532 and 800 nm generally show similar patterns of thermal diffusion, due to the small differences between the values of the thermal absorption coefficients for the layers of the skin. In general, thermal laser treatment is based on two steps: During the early stages of thermal laser treatment, it seems that the tissue was heated directly within the depth of light absorption and that the heat did not diffuse deeply into the tissue; however, with time, it appears that the heat diffused deeper into the tissue. In addition, blood perfusion plays a cooling role that may prevent the temperature from increasing above the permissible limits.

On the other hand, simulation results suggest that the maximum temperature is located within the central exposure area, while it seems that the temperature tends to decrease when moving away from the exposure area due to the surrounding biological tissues when using all the studied wavelengths and specific intensities. Temperatures can increase within the deeper layers due to the diffusion of heat from the central spot to the surrounding tissues, which is a cooler area compared to the exposure area.

It is noteworthy that the longitudinal heat distribution plays a crucial role in the diffusion of heat across the various tissue layers. Additionally, it should be acknowledged that the size of the heat focus spot decreases as the wavelength increases. The longitudinal temperature gradient field demonstrates greater effectiveness compared to the transverse direction, owing to the upper surface of the skin being subjected to a convective thermal restriction, leading to the continual conveyance of heat towards the surrounding tissues. It should also be mentioned that laser intensity levels can lead to an increase in energy absorption, which causes an increase in temperature and thermal diffusion at all wavelengths.

Figure 7. Distribution of the equivalent temperatures (°C) and equivalent stresses (MPa) across the model layers at a time of 600 s, intensities of 1 and 2 W/mm², and wavelength of 800 nm.

Figure 7. Distribution of the equivalent temperatures (°C) and equivalent stresses (MPa) across the model layers at a time of 600 s, intensities of 1 and 2 W/mm², and wavelength of 800 nm.

3.2. Distribution of Equivalent Stresses

The simulation results are shown in Figure 2, Figure 3, Figure 5 and Figure 7, and they demonstrate the location of the maximum equivalent stress values on the skin’s upper surface when exposed to the laser in all scenarios, regardless of whether the maximum temperature is situated on the surface or within the tissue. Figure 8 shows the equivalent stress variations (at z = 0 and r = 0) with time and a wavelength of 800 nm, laser intensities of 1 and 2 W/mm², and a beam size of 1 mm. An exponential increase in stress was observed with time, with a difference of not less than 50% between the two studied strengths. Moreover, it is noticed that equivalent stresses reach a stable state after a time of about 400 s, which confirms the dependence of the mechanical effect on the thermal effect associated with laser radiation.

Figure 8. Variations of equivalent stresses (at z = 0 and r = 0) with time at 800 nm wavelength, 1 and 2 W/mm² laser intensities, and 1 mm beam width.

Figure 8. Variations of equivalent stresses (at z = 0 and r = 0) with time at 800 nm wavelength, 1 and 2 W/mm² laser intensities, and 1 mm beam width.

When studying the equivalent stresses according to the radial direction, it is noticed that the maximum value appears in the central area and a decrease in the values with decreasing distance from the surface, while the maximum values fade in the longitudinal direction when the depth increases beyond the distance of 3 mm. In addition, the differences between 1 and 2 W/mm² intensities are more pronounced according to the diagonal direction than according to the longitudinal or axial direction.

3.3. Distribution of the Total Displacement

Figure 9 shows the three-dimensional distribution of the overall displacements in the skin model after exposure to wavelengths of 532 and 800 nm, intensities of 1 and 2 W/mm², and bandwidths of 1 and 2 mm at an exposure time of 600 s. It is observed that maximum displacement is located at the edges of the skin surface at the wavelengths 532 and 800 nm, but this region of maximum displacement goes deeper when the level of the intensity of the laser beam increases. When the width of the laser beam is increased to 2 mm at the same intensity, it is noticed that the maximum value moves towards the depth, while changing the intensity does not change the position of the maximum value. It is observed that the maximum value of the total displacement increases at a beam width of 2 mm compared to the value of 1 mm.

Figure 9. Distribution of model discplacements in skin model layers at wavelengths 532 and 800 nm, laser intensities of 1 and 2 W/mm², and beam widths of 1 and 2 mm.

Figure 9. Distribution of model discplacements in skin model layers at wavelengths 532 and 800 nm, laser intensities of 1 and 2 W/mm², and beam widths of 1 and 2 mm.

3.4. Optimization of Thermomechanical Effect

The effect of processing parameters (wavelength W, beam size S, laser intensity I, and exposure time T) on temperature and stress distributions was studied in the skin using Box-Behnken design optimization (

Table 2. The Box-Behnken design of thermomechanical optimization (parameters, levels, and responses).

Run	Intensity I [W/mm2]	Size S [mm]	Time T [s]	Wavelength (λ) W [nm]	Temperature Tmax [°C]	Stress Smax [MPa]
1	3	3	30	532	42.56	0.5802
2	3	1.75	315	532	66.22	0.5708
3	1	0.5	315	800	45.52	0.908
4	1	3	315	532	45.55	0.1715
5	5	0.5	315	532	87.22	0.9674
6	1	0.5	315	532	45.55	0.1715
7	1	1.75	30	800	38.29	0.174
8	3	1.75	315	800	66.08	2.93
9	5	1.75	600	800	92.92	5.53
10	1	1.75	600	532	46.64	0.1732
11	5	1.75	30	532	46.83	0.9834
12	5	0.5	315	800	86.99	4.97
13	3	1.75	315	532	66.22	0.5708
14	3	0.5	600	800	69.74	3.28
15	1	3	315	800	45.52	0.908
16	3	1.75	315	800	66.08	2.93
17	3	1.75	315	800	66.08	2.93
18	1	1.75	30	532	38.29	0.1769
19	3	1.75	315	532	66.22	0.5708
20	1	1.75	600	800	46.60	1.02
21	5	3	315	532	87.22	0.9674
22	5	3	315	800	86.99	4.97
23	3	0.5	600	532	69.89	0.565
24	3	3	600	800	69.74	3.28
25	3	0.5	30	532	42.56	0.5802
26	3	3	600	532	69.89	0.565
27	5	1.75	600	532	93.19	0.9625
28	5	1.75	30	800	46.83	0.96
29	3	3	30	800	42.56	0.57
30	3	0.5	30	800	42.56	0.57

). The optimal parameter levels corresponding to the minimum values for each temperature were also determined T_max (≤43 °C) and stress S_max (≤0.2 MPa).

Table 2 shows the experiment array consisting of 30 runs according to the Box-Behnken design to evaluate the effect of the following parameters: wavelength (W), intensity (I), beam diameter expressed as beam size (S), and skin exposure time to the laser (T). On the change in each of the two responses studied (maximum temperature T_max and stress S_max).

Table 3. Analysis of variance (ANOVA) to estimate the parameters’ effect for each response (T_max, S_max) and determine p-values (at α = 5%).

	Source	DF	Adj SS	Adj MS	F-Value	p-Value	Contrubtion	R-sq
Tmax	I	2	4769.00	2384.50	57.24	0.000	50.1%	90.27%
	S	2	0.05	0.02	0.00	0.999	…
	T	2	3718.45	1859.23	44.63	0.000	39.5%
	W	1	0.08	0.08	0.00	0.965	…
	Error	22	916.40	41.65			9.7%
	Total	29	9414.06
Smax	I	2	17.2383	8.6191	9.32	0.001	24.1%	71.62%
	S	2	0.0000	0.0000	0.00	1.000	…
	T	2	9.1414	4.5707	4.94	0.017	12.8%
	W	1	24.9403	24.9403	26.96	0.000	34.7%
	Error	22	20.3488	0.9249			28.3%
	Total	29	71.6918

shows the results of the ANOVA analysis of variance between the averages of the studied response under the influence of the four influential parameters using IBM SPSS Statistics version 26.0 (SPSS Inc., Chicago, IL, USA), clarifying the presence or absence of statistically significant differences depending on the probability value (p-value) at a confidence interval (α = 0.05). The results of the analysis of variance showed that there was no statistical significance for the effect of both beam size S (p = 0.999 > 0.05) and wavelength (W) (p = 0.965 > 0.05). At temperature T_max, in contrast, both the beam intensity (I) and the exposure time (T) have statistical significance, with a combined contribution of the two parameters estimated at 89.6%. Figure 10a illustrates the average impact of the mentioned parameters on (T_max). The figure indicates that the response remains unaffected by S and W. In contrast, the results of the analysis of variance showed that there is no statistical significance for the effect of the beam size S (p = 1.00 > 0.05) only on the stress value S_max. In contrast, the wavelength (W), the beam intensity (I), and the exposure time (T) have statistical significance. The combined contribution of the three parameters is estimated at 71.6%. Figure 10b also shows the average effect of the mentioned parameters on (S_max).

The response surface for both skin temperature T_max (Figure 10a) and stress state S_max Figure 10b can be presented in terms of beam intensity (I) and exposure time (T). The value of both T_max and S_max can also be extrapolated as a function of beam intensity (I) and exposure time (T) in Figure 11c,d, respectively. Figure 11c,d show the possibility of obtaining minimum values for the skin temperature T_max and its stress state S_max at low values of beam intensity (I) and exposure time (T).

Equations (11)–(14) are the regression equations as a function of the processing parameters (I, S, T) for predicting the value (S_max)₅₃₂, (S_max)₈₀₀, (T_max)₅₃₂, and (T_max)₅₃₂, respectively.

(S_max)₅₃₂ = 0.409 − 0.059 I − 0.006 S + 0.00145 T + 0.0003 I² + 0.002 S² − 0.000006 T² + 0.0000 I·S + 0.000813 I·T

(11)

(S_max)₈₀₀ = −1.191 + 0.581 I − 0.006 S + 0.00622 T + 0.0003 I² + 0.002 S² − 0.000006 T² + 0.0000 I·S + 0.000813 I·T

(12)

(T_max)₅₃₂ = 28.93 + 3.28 I − 0.18 S + 0.0760 T + 0.023 I² + 0.051 S² − 0.000124 T² + 0.000 I·S + 0.01662 I·T

(13)

(T_max)₈₀₀ = 29.03 + 3.24 I − 0.18 S + 0.0757 T + 0.023 I² + 0.051 S² − 0.000124 T² + 0.000 I·S + 0.01662 I·T

(14)

The determination factors are (R-sq)_Smax = 93.48% and (R-sq)_Tmax = 97.89%.

By applying the response optimization feature in the Box-Behnken approach, it is possible to query the levels of parameters corresponding to a specific goal after analyzing and comparing the responses corresponding to the distribution of levels in the experiment array. In this study, the appropriate goal is to determine the minimum limits for each of the two responses: skin temperature (T_max)_min and stress (S_max)_min. The results of the analysis indicated that the minimum values for both T_max and S_max can be obtained according to the following levels: (I = 1, S = 1.8, T = 30, W = 800). These levels were tested in the simulation model (FEA), and the response results were: (T_max)_min = 38.29 [°C], (S_max)_min= 0.174 [MPa].

Figure 10. Main effects of parameter levels on: (a): skin temperature Tmax, and (b) stress S_max.

Figure 10. Main effects of parameter levels on: (a): skin temperature Tmax, and (b) stress S_max.

Figure 11. Box-Behnken optimization results (a): Response surface S_max as a function of beam intensity I and exposure time T. (b): Response surface T_max as a function of beam intensity I and exposure time T. (c): The effect of changing the values of beam intensity I on changing the values of temperature T_max and stress S_max. (d): The effect of changing the values of exposure time T on changing the values of temperature T_max and stress S_max.

Figure 11. Box-Behnken optimization results (a): Response surface S_max as a function of beam intensity I and exposure time T. (b): Response surface T_max as a function of beam intensity I and exposure time T. (c): The effect of changing the values of beam intensity I on changing the values of temperature T_max and stress S_max. (d): The effect of changing the values of exposure time T on changing the values of temperature T_max and stress S_max.

4. Discussion

In this study, we investigated the thermo-optical and optomechanical effects of short-pulse laser irradiation on a layered skin model, aligning our methodology with previous research for robust comparisons. Our analysis focused on determining the heat-affected region resulting from laser radiation across various laser variables, including wavelength, radiation intensity, laser beam size, and exposure time. To capture the intricate dynamics, we conducted numerical simulations employing a finite element model within the COMSOL Multiphysics program. This choice was motivated by the program’s capabilities in studying thermal diffusion, coupled with static and kinetic mechanical analyses.

The accuracy of our model was rigorously verified by examining temperature distributions, mechanical stresses, and deformations. Notably, our results underscored the pivotal roles of exposure time and laser intensity in achieving the requisite treatment temperature. Furthermore, we observed a balanced interaction between intensity and exposure time in attaining the optimal temperature. These findings align closely with the observations of Wongchadakul et al. [30].

Our analysis of short laser pulses demonstrated their ability to deliver high energy to target sites efficiently, producing a well-defined heat-affected area. We delved into the thermal and mechanical effects over time intervals ranging from 30 to 600 s, focusing particularly on the thermomechanical effect, which emerged as predominant in laser-tissue reactions.

The study’s parametric analysis encompassed the frequency, beam size, intensity, and exposure time of pulse radiation, elucidating the photomechanical interaction. The coupling feature in COMSOL Multiphysics facilitated simulations, revealing that equivalent stress peaked around 2 MPa, gradually diminishing with distance from the irradiation point. These findings align closely with the observations of Wongchadakul et al. [30] and underscore the significance of temperature-stress relationships, especially in the axial direction.

Utilizing a second-order polynomial model, we quantitatively evaluated the effects of laser intensity, beam size, wavelength, and exposure time on temperature prediction. Our modeling efforts revealed that a 532 nm laser with an intensity of 1 W/mm², a beam size of 3 mm, and an exposure time of 309 s could produce a maximum temperature of 43 °C.

Our study highlights the need for further experimental investigations to elucidate variable blood perfusion rates and assess thermal damage, particularly in tumor tissues with high blood perfusion rates. We suggest employing established thermal damage equations, such as the first-order Arrhenius equation [35], modified Arrhenius equation [36], or temperature-dependent time-delay equation [37], to refine predictions and accommodate tissue-specific variations. Moreover, our investigation could be expanded by employing controlled micro/nano-scale periodic surface structures in conjunction with femtosecond pulsed lasers to enhance tumor ablation [38].

Looking ahead, our thermomechanical model can be extended to analyze tumor cell necrosis using quantitative models under imaging guidance and track drug delivery employing shorter laser pulses. Incorporating thermal damage details and Von Mises stresses from temperature-time history would enhance computational predictions, supporting pre-clinical and clinical assessments. Moreover, integrating accurate blood perfusion information from imaging techniques can refine predictions derived from the Pennes Bioheat Equation, fostering advancements in understanding micro- and nano-scale heat transfer in biological tissues.

5. Conclusions

In conclusion, our study delved into the photothermal and thermomechanical effects resulting from short-pulse laser irradiation on normal tissues. By analyzing a range of laser variables, including wavelength, intensity, beam size, and exposure time, we explored their influence on the heat-affected region within tissues. Utilizing a three-layered, three-dimensional model implemented in a polar coordinate system, numerical simulations were conducted using COMSOL Multiphysics.

Integration of the Pennes bioheat transfer model, Beer-Lambert law, and Hooke’s law facilitated the simulation of coupled biophysics phenomena, enabling the examination of temperature and stress distributions following laser radiation. Qualitative verification of our model’s accuracy was achieved through comparisons of temperature and mechanical variations with relevant studies.

Box-Behnken analysis revealed that beam size (S) exhibited negligible impact on response variables, with p-values exceeding 0.05. Temperature (T_max) demonstrated sensitivity to both beam intensity (I) and exposure time (T), collectively contributing to 89.6% of the observed variation. Conversely, while beam size (S) showed no significant effect on stress value (S_max), wavelength (W), beam intensity (I), and exposure time (T) collectively accounted for 71.6% of observed variation in equivalent stress (Smax).

Our research envisions extending the photomechanical model to analyze tumor cell necrosis and drug delivery using shorter laser pulses, enhancing computational predictions, and refining blood perfusion data from imaging.

In summary, this research advances our understanding of micro- and nano-scale heat transfer in biological tissues, offers promising applications in laser-based medical interventions, and provides deeper insights into thermal and mechanical interactions within living organisms.