Theoretical Study on Cryogen Spray Cooling in Laser Treatment of Ota’s Nevus: Comparison and Optimization of R134a, R404A and R32

Abstract

Cryogen spray cooling (CSC) could be applied clinically for the laser therapy of Ota’s nevus, a dermal hyperplastic pigmented disease with a morbidity rate of 0.1–0.6% in the Asian population. An accurate, efficient, complete simulation system that considers the entire spray cooling process, including cryogen flow in the tube nozzle, spray dynamics and internal phase change heat transfer (cold injury) in skin tissue, was established to determine suitable cryogen and cooling parameters. The optimum spray distances for R134a, R404A and R32 were determined to be 66.0, 43.1 and 22.5 mm, respectively. The corresponding maximum surface heat fluxes were 363.5, 459.9, and 603.6 kW∙m⁻², respectively. The maximum surface heat flux of R32 with small spray distance was 1.66 times as large as that of R134a, indicating the potentially good cooling performance and precise targeted cooling of R32 during the laser therapy of Ota’s nevus. The cooling durations that caused cold injury of skin tissue were 2.3, 1.4, and 1.1 s for R134a, R404A, and R32, respectively. The interval between CSC and laser irradiation was optimized to 90–162 ms for R134a, R404A and R32, in consideration of the cooling effect, depth, uniformity, and risk of cold injury.

1. Introduction

Ota’s nevus is a dermal disorder of hyperpigmentation with a morbidity rate of 0.1%–0.6% in the Asian population [1]. It mostly distributes on the face and neck, and is caused by the abnormal migration of primitive melanocytes in skin tissue at the embryonic stage. As shown in Figure 1, the human skin is a multi-layered structure comprised of the epidermis, dermis and subcutaneous tissue (from the outermost to the innermost layer). Normal melanin with a diameter of 0.1–0.2 μm is buried in the epidermal layer, and its content determines skin color [2]. Melanin lesions are buried in the dermis layer and can be seen as irregular spindles with diameters of 0.25–0.65 μm on histopathologic slides [3]. Laser surgery with 532 nm/1064 nm Nd:YAG and 755 nm alexandrite lasers is implemented in clinical settings due to selective photothermolysis [4]. Laser heating can thermally damage melanin lesions, accompanied with the explosive vaporization of melanosome. The surrounding tissues survive maximally. However, the high laser absorption of epidermal normal melanin often causes adverse side effects, such as epidermal thermal injury, hypopigmentation and bleeding [5]. Effective skin cooling technologies must be developed urgently to solve these problems.

Cryogen spray cooling (CSC) has a high cooling capacity and spatial and temporal selectivity. It allows for the rapid evaporation of the liquid film (completely evaporation time ≤30 ms [6]) on the skin surface due to the large temperature difference between liquid film and environmental air, and is nontoxic to the human body. Thus, CSC can facilitate the pre-cooling of the epidermis before laser irradiation. CSC has been successfully used in the laser therapy of port-wine stain (PWS) birthmarks to avoid thermal injury on the epidermis caused by the competitive laser absorption between epidermal melanin and dermal blood. However, the potential application of CSC in the laser therapy of Ota’s nevus has not yet been investigated. The CSC parameters used in PWS cannot be directly applied to Ota’s nevus owing to the large differences in tissue composition and structure, as well as the special requirement in terms of larger cooling capacity. The potential clinical application of CSC in the laser therapy of Ota’s nevus aims to avoid the competitive laser absorption between epidermal and dermal melanin, which cannot be realized by changing the laser wavelength. The selection of cryogen for safe and clinical application should follow these rules [7]: (1) high thermal conductivity for rapid heat removal from skin tissue; (2) low boiling point and high latent heat for enough cooling and the fast evaporation of the liquid film on the skin surface; (3) no adverse health effects or danger to the human body, as well as environmental effects. In the laser treatment of PWS, R134a with a normal boiling point (T_b) of −26.1 °C was approved by the U.S. Food and Drug Administration (FDA) for clinical utilization, but its insufficient cooling capacity is not applicable for darkly pigmented human skin [8]. R404A with a low T_b (−46.5 °C) may be a substitute for R134a, because it produces a powerful cooling effect with an increment of 17.8% in surface heat flux [9,10]. The nontoxicity and safety in the human body of R404A were approved by the American Society of Heating Refrigerating and Airconditioning Engineers (ASHRAE) [11]. Furthermore, the histological observation investigation conducted by Dai et al. [8] demonstrated that the utilization of R404A with spurt durations up to 300 ms did not produce obvious morphological changes to dark-pigmented human skin. For R32, animal studies indicated that the acute inhalation toxicity is very low and no adverse human health effects were reported, as documented by the International Occupational Safety and Health Information Centre (CIS) [12]. Although R32 is flammable, the local R32 vapor concentration is around 0.24–4.76 g/m³ in a supposed space 1 m³, under the clinical spurt duration (5–100 ms) [13] and mass flow rate (47.63 g/s) parameters, which is far less than the lower flammability limit (equated mass concentration, 306 g/m³) [14], indicating the safety in terms of flammability. Moreover, the suction cup joined with a vacuum hand pump [15] can be implemented to collect the cryogen vapor to ensure its safe application for the human body. R134a and R404A have zero ozone-depleting potential, but their extremely high global warming potentials (GWP; 1430 for R134a, 3850 for R404A) necessitate the substitution of a new cryogen. In the laser surgery of PWS, most of the laser energy is absorbed by hemoglobin, whereas a small amount is absorbed by epidermal melanin; this can be realized by selecting an appropriate laser wavelength. However, the laser absorption is similar between normal epidermal melanin (µ_a = 750–800 m⁻¹, λ = 450–800 nm [16]) and abnormal dermal melanin (µ_a = 700–800 m⁻¹, λ = 450–1050 nm [17]) in Ota’s nevus, resulting in a high temperature increase in the epidermal melanin and surrounding tissue. Therefore, a large cooling capacity is required for the CSC-assisted laser therapy of Ota’s nevus. R32 with a low boiling point (−51.9 °C) and low GWP (675) may be used as a potential cryogen in the laser therapy of Ota’s nevus. However, little work has focused on R32′s spray dynamics, cooling effect and parameter optimization for practical applications. Although the spray distances for R134a and R404A spray cooling have been optimized for the laser therapy of PWS [18], the spray distance and other parameters should be re-optimized for Ota’s nevus due to the large difference in tissue composition and structure.

CSC is a complicated process that merges the gas–liquid flow in the capillary tube, the atomization, evaporation and breakup of spray droplets, the boiling heat transfer on the skin surface, and the phase-change heat transfer (cold injury) inside bio-tissue. Experimental investigations cannot effectively reveal spray cooling mechanisms, but numerical simulation is a powerful tool for this task. Nilpueng and Wongwises [19] proposed the homogeneous flow model to simulate the internal cryogen flow in a straight-tube nozzle by dividing the flow state into subcooled, metastable, and two-phase regions. Many theoretical models targeted at spray dynamics have been developed to simulate the atomization, evaporation and breakup of spray droplets. For example, Zeng and Lee [20] proposed an atomization model assuming that the droplet primary breakup is affected by air drag force and bubble growth on the droplet surface. Witlox et al. [21,22] established the relationship between droplet diameter and superheat degree under mechanical transition and atomization by using a piecewise linear function. Other researchers focused on the initial droplet diameter distribution [23]. Witlox et al. [21] improved the Rosin–Rammler distribution model to simulate spray dynamics and found that after atomization, evaporation and secondary breakup occur when a droplet flies. Ranz and Marshall [24] identified the classical convective heat and mass transfer correlation for water droplet evaporation. Afterward, Spalding [21], R-Y [24] and N-G [25] models were developed for consideration of the blowing effect. Zhou et al. [26] proposed an N-G-R-M model by comprehensively considering the convection and blowing effect to predict droplet evaporation with different evaporation rates accurately. Droplets break up due to the shear force of gas. Four available droplet breakup models, namely, Taylor analogy breakup (TAB) [27], wave breakup [28], Kelvin–Helmholtz [29] and stochastic secondary droplet (SSD) [30], can be used to simulate droplet breakup under different Weber numbers. Aiming at the simulation of bio-tissue phase-change heat transfer upon CSC, Li et al. [31] also constructed a multi-scale heat transfer model based on cellular frostbite theory to simulate cell dehydration and intracellular ice formation. This model has been proven successful in predicting safe spray durations of CSC to prevent cold injury of skin tissue. However, little work has focused on CSC for the laser therapy of Ota’s nevus.

The aforementioned investigations only focused on a single physical process during spray cooling. In the present work, an accurate, efficient, complete simulation system that considers the entire process of cryogen spray was constructed by considering the processes of gas–liquid flow through a nozzle tube, the atomization, evaporation and breakup of spray droplets, and phase-change heat transfer (cold injury) inside bio-tissue, using well established mathematical models for the literature. Besides, the coupling of the heat transfer model on the skin surface constructed a bridge between the micro-spray droplet characteristics and surface boiling heat transfer, which can prevent the direct simulation of the complicated droplet impact cooling process, thus highly reducing the computation load. With the complete simulation system, the simulation results were compared with experiments in terms of droplet temperature and diameter, as well as the surface temperature of skin tissue. Then, the spray characteristics of R134a, R404A and R32 spray cooling were compared. Afterward, the spray distances of CSC with three cryogens were optimized to determine the best cooling capability and investigate the potential application of the new cryogen R32 in clinics. In addition, the interval between CSC and laser irradiation was determined by comprehensively considering cooling performance, cooling depth, cooling uniformity and risk of cold injury to provide theoretical guidance for the safe and precise clinical application of CSC in the laser treatment of Ota’s nevus.

2. Mathematical Modeling for Cryogen Spray Cooling

To establish the complete simulation system, a 1D homogenous flow model was employed to solve the two-phase cryogen flow in a tube [19]. The discrete phase model (DPM) was employed to simulate spray atomization, droplet evaporation and breakup under the Euler–Lagrange framework. Cell dehydration and intracellular ice formation under CSC were also simulated using the multi-scale bio-tissue heat transfer model based on cellular frostbite theory [31]. Given that the present study mainly concerns the simulation of the spray field and heat transfer inside skin tissue, the details of the 1D homogenous flow model [19] are not discussed in this paper.

2.1. Discrete Phase Model for Spray Dynamics

2.1.1. Droplet Motion, Breakup and Evaporation

The DPM method tracks individual droplet motion under the Lagrangian framework to simulate droplet motion, breakup and evaporation. In consideration of the effects of gravity and gas drag force on a discrete droplet, the motion of the droplet can be calculated by

$$\frac{d{u}_d}{dt}=f_D\left(u_a-u_d\right)+\frac{g\left({\rho }_d-{\rho }_a\right)}{{\rho }_d}$$

(1)

where u_d and u_a are the velocities of the droplet and ambient gas, g is the acceleration of gravity, ρ_d and ρ_a are the densities of the droplet and ambient gas, nad f_D(u_a−u_d) is the drag force of the droplet per unit mass from surrounding gas. f_D is the drag force coefficient introduced in Morsi and Alexander’s work [32].

Violent flashing atomization often occurs at the nozzle exit during CSC, because of the sudden release of high-pressure cryogen spraying into a low-pressure environment. It is very difficult to directly simulate such a complicated flash, including bubble formation on the droplet surface and the droplet breakup induced by the bubble explosion, coupled with the phase-change heat and mass transfer. Furthermore, the present work mainly focuses on the fully-developed spay after flashing atomization at the nozzle exit, because the clinical spray distance is often selected to be larger than 10 mm to obtain the good cooling performance. Thus, a piecewise linear model proposed by Witlox et al. [22] concerning the effect of superheat degree on different atomization forms (mechanical, transition and flashing atomization) was adopted in the present simulation to obtain the droplet diameter distribution after flashing atomization occurs. The later model validation indicates the rationality of such treatment, which also avoids the complicated simulation of flash. The droplet temperature after flashing atomization is assumed to be the boiling point in an atmospheric environment. The stochastic secondary droplet (SSD) model [30] was applied to simulate droplet breakup in the fully-developed spray induced by the gaseous shear force, due to the model’s advantage of small computation load. In the SSD model, the breakup of numerous droplets with a high Weber number is regarded as a discrete random event, in which the breakup degree is related to the size of the parent droplet.

For flashing spray, the cryogen droplet temperature decreased to the boiling point at the nozzle outlet. Droplet evaporation dominated during droplet flight with a high velocity in the fully-developed spray. The convection/diffusion controlled model proposed by Miller et al. [33,34] was used to calculate the mass transfer during droplet evaporation, as follows:

$$\frac{d{m}_d}{dt}=k_c{A}_d{\rho }_{d}ln\left(1+B_m\right)$$

(2)

where m_d and A_d denote droplet mass and surface area. B_m is the Spalding mass transfer number. The mass transfer coefficient, k_c, was computed by the N-G-R-M model that considers the convective and blowing effect of air at the droplet interface, proposed by Zhou et al. [26]. The heat transfer between droplet and ambient gas was simulated with the lumped parameter method on the basis of the assumption of uniform density and temperature inside the droplet. The energy conservation equation is

$$m_d{c}_{p,d}\frac{d{T}_d}{dt}=h{A}_d\left(T_a-T_d\right)-\frac{d{m}_d}{dt}{h}_{fg}+A_d{\epsilon }_d\sigma \left(T_r^4-T_d^4\right)$$

(3)

where T_d and T_a are the droplet and gas temperature, c_p,d and h_fg are the specific heat and latent heat of the droplet, and ε_d and T_r are the spectral emissivity and radiation temperature of the droplet. ε_d is assumed to be 0.96, which is the same as that of water because of the similar optical properties, since there are no measured data for cryogens. σ is the Stefan–Boltzmann constant, i.e., σ = 5.67 × 10⁻⁸ W·m⁻²·K⁻⁴. The key parameter for droplet evaporation, i.e., convective heat transfer coefficient (h), can be determined by the Nusselt number correlation proposed by Sergei S. Sazhin [34],

$$Nu=\frac{\ln \left(1+B_t\right)}{B_t}\left(2+0.552R{e}^{1/2}+P{r}^{1/3}\right)$$

(4)

where Nu = hd/λ_a, Re = ρ_dd|u_d−u_a|/η, Pr = c_p,dη/λ_a, B_t = c_p,a(T_a−T_d)/h_fg. η and h_fg are the viscosity and latent heat of droplet. λ_a and c_p,a are the heat conductivity coefficient and specific heat of gas. The Nusselt number correlation (Equation (4)) was obtained by Sergei S. Sazhin [34] using dimensional analysis, and is applicable to the convective heat transfer for moving droplet evaporation with high velocity, such as cryogen droplet evaporation during CSC. Its accuracy has been validated by Zhou et al. [26]. However, the temperature-dependent thermal properties of the droplet and the surrounding air should be used to accurately implement the Nusselt number correlation in the current study.

2.1.2. Governing Equations for Continuous Gas

The simulation of the gas phase was conducted under the Eulerian framework by considering mass, momentum and energy conservation. The mass conservation equation is expressed as

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho u\right)=S_m$$

(5)

S_m is the mass source term induced by droplet evaporation. The momentum conservation equation can be obtained by solving the Navier–Stokes equation

$$\frac{\partial }{\partial t}\left(\rho u\right)+\nabla \cdot \left(\rho uu\right)=-\nabla p+\nabla \cdot \overline{\overline{\tau }}+\rho g+F$$

(6)

where F is the momentum source term described by the change of droplet momentum in the control volume. p is the pressure, and $\overline{\overline{\tau }}$ is the tensor of stress. The energy conservation equation is

$$\frac{\partial }{\partial t}\left(\rho E\right)+\nabla \cdot \left(u\right(\rho E+p\left)\right)=\nabla \cdot \left(k_{eff}\nabla T_a-{\displaystyle \sum _i{h}_i{J}_i+(\overline{\overline{\tau }}\cdot u)}\right)+\rho g+S_h$$

(7)

where k_eff pertains to effective thermal conductivity, and S_h is the energy source term. h_i and J_i are the enthalpy and diffusion rate of the ith component in the gas phase, respectively. The two-way coupling of mass, momentum, and energy transfer between discrete droplets and the continuous gas were established with the source terms in Equations (5)–(7).

2.2. Multi-Scale Heat Transfer Model for Skin Tissue Frostbite

When the tissue cell temperature is below 0 °C, water in the cell starts to freeze. The cell state can be regarded as unfrozen, partially frozen, or completely frozen. In the partially frozen state, extracellular ice forms first, yielding to the water transported on the cell surface because of the concentration difference. The cell death induced by intracellular ice formation (IIF) occurs with a sufficiently high cooling rate. Meanwhile, if the cooling rate is relatively slow, cell dehydration is induced because of water transport on the cell surface. Differing between traditional engineering materials, the macro-scale heat transfer and the micro-scale cellular-level biophysical events (intracellular ice formation and cell dehydration) must be considered simultaneously when solving the heat and mass transfer dynamics to evaluate the cold injury of skin tissue under CSC, and this task can be executed with the multi-scale heat transfer model developed by Li et al. [31]. Macro-scale heat transfer equations were implemented to calculate the temperature distribution in the skin tissue, considering the heat source term induced by intracellular ice formation and cell dehydration,

$${\rho }_j{c}_j\frac{\partial T}{\partial t}=\nabla \cdot ({\lambda }_j\nabla T)+a\left[(1-\mathrm{PIIF})\frac{2\pi R_{ec}\rho L}{\Delta x^2}\frac{d{R}_{ec}}{dt}+\frac{(V_c-V_{ns})\rho L}{\Delta x^3}\frac{d\mathrm{PIIF}}{dt}\right]$$

(8)

where ρ, c and λ are the density, specific heat and heat conduction of tissue coefficient. The subscript j denotes the unfrozen (uf), completely frozen (cf) and partially frozen states (pf), respectively. V_c and V_ns are the cell volume and total volume of Krogh unit. PIIF is the probability of IIF, and L is the latent heat of water. ∆x and R_ec are the Krogh unit dimension and equivalent radius of the extra-cellular space. a equals 0 for the unfrozen and completely frozen states, and 1 for the partially frozen sate.

The micro-scale mass transfer equation can be implemented in the Krogh unit that represents a specific space with a similar tissue structure, containing certain numbers of cells and surrounded by interstitial fluid. The mass transfer inside the Krogh unit was considered by solving the radius change of the cylindrical extracellular space because of the water transport on the cell surface. The ice formation was coupled with the macro-heat transfer equation as the latent heat source term. PIIF is calculated by the following equation [35],

$$\mathrm{PIIF}=1+\exp \left[{\displaystyle \underset{T_f}{\overset{T}{\int }}A\mathsf{\Omega }{\left(\frac{T}{T_{f,0}}\right)}^{1/2}\frac{{\eta }_{c,0}}{{\eta }_c}\frac{A_c}{A_{c,0}}\exp \left(\frac{-k_c{(T/T_{f,0})}^{1/4}}{{(T-T_f)}^2{T}^3}\right)}dt\right]$$

(9)

where T_f is the water freezing point, and A_c is the valid superficial cell area for water transport. k_c and Ω are the thermodynamic and kinetic parameter for crystal nucleation on the cell surface. η_c is the intracellular solution viscosity, determined by Bischof and Rubinsky’s work [36]. By solving the multi-scale heat transfer equations, the temperature distribution of skin tissue, the dynamic change in R_ec, and the PIIF can be obtained and used to evaluate cold injury during CSC. Cold injury occurs when the PIIF is 1 or R_ec reaches its maximum value R_ec,max, which is determined by the equivalent radius of Krogh unit volume without the water component [36].

2.3. Simulation Conditions and Model Validation

Figure 2 indicates that the simulation considered a 3D axisymmetric cylinder domain with a radius of 25 mm and a height of 103 mm, consisting of spray (Z_s = 100 mm) and skin tissue (Z_t = 3 mm) domains. The origin of the coordinates was located at the nozzle exit, and the cryogen was spurted onto the skin surface vertically. In accordance with the histological anatomy of Ota’s nevus [37], the skin tissue domain consisted of four planar layers, namely, the epidermis, dermis without melanin, dermis with melanin, and dermis without the melanin layer. The Krogh unit with a dimension of (∆x × ∆x × ∆x) was used to represent the specific space containing certain numbers of cells and surrounding interstitial fluid (∆x = 22 μm), determined by the average distance between two adjacent sinusoids in the scanning electron micrographs of the tissue [36]. In accordance with the Asian skin structure, the epidermis thickness was assumed to be 60 μm. The total thickness of the dermis layer was 2940 μm. Melanin volume fractions of 15% and 35% were uniformly distributed in the epidermis and dermis melanin layers, and no melanin existed in the other layers.

The experimental investigation of heat transfer performance during R404A spray cooling by Zhou et al. [38] indicated that a spurt duration longer than 50 ms has an insignificant effect on cooling capacity. As such, a spurt duration of 50 ms was selected in the current work. Three cryogens of R134a, R404A and R32 were compared in the present work. The physical properties of three cryogens at 25 °C and 1 atm are compared in

Table 1. Physical properties of R134a, R404A and R32 (25 °C, 1 atm).

Cryogens	ρd (kg/m3)	cp,d (J/kg∙K)	λd (W/m∙K)	Tb (°C)	psat (MPa)	σ (N/m)	η (Pa∙s)	hfg (J/kg)
R134a	1206.7	1424.6	0.0811	−26.1	0.67	0.0081	0.00019	216,000
R404A	1044.1	1542.3	0.0646	−46.5	1.25	0.0046	0.00013	230,260
R32	961.0	1936.7	0.1259	−51.9	1.69	0.0068	0.00011	270,910

, and the thermal properties of the skin components are documented in

Table 2. Thermal properties of skin components [39].

Properties	States	Epidermis	Dermis
ρ (kg·m−3)	Unfrozen	1120	1090
	Completely frozen	1120	921
c (J·kg−1·K−1)	Unfrozen	3200	3500
	Completely frozen	3200	1230
λ (W·m−1·K−1)	Unfrozen	0.34	0.41
	Completely frozen	0.34	2.25

.

The constant flow velocity boundary condition was implemented at the nozzle inlet. The pressure outlet boundary was employed at the surface of the spray domain. Direct simulation considering the complex processes of spray droplet impact cooling, heat transfer inside the liquid film, and interactions between droplet and thin liquid film on the skin surface is very difficult. Alternatively, the heat transfer model for the skin surface (Z_t = 0) was applied to simulate the integrated cooling effect during CSC, as follows:

$$-\lambda {\frac{\partial T}{\partial z}|}_{Z_t=0}=q(\varphi ,R,t)$$

(10)

The surface heat flux was calculated with the general correlation proposed by Tian et al. [10] as follows:

$$\frac{B{i}^{\ast }}{B{i}_{\max }^{\ast }}=\{\begin{array}{ll}\frac{F{o}_{\mathrm{s}}}{F{o}_{\mathrm{s},}{}_{\max }}& 0\le \frac{F{o}_{\mathrm{s}}}{F{o}_{\mathrm{s},}{}_{\max }}\le 1\\ 1.265{\mathrm{e}}^{(-\frac{F{o}_{\mathrm{s}}}{1.9F{o}_{\mathrm{s},}{}_{\max }})}+0.25& 1<\frac{F{o}_{\mathrm{s}}}{F{o}_{\mathrm{s},}{}_{\max }}\le 14\end{array}$$

(11)

where Bi* is the spray Biot number that measures the relationship between surface convective heat transfer resistance and internal thermal resistance of the skin during CSC [10] (Bi* = h*δ/λ, h* = q/(T₀ − T_d)). T₀, δ and λ denote the initial temperature, thickness, and heat conductivity coefficient, respectively. T_d is the spray droplet temperature (Fo_s = at/δ²). a is the thermal diffusivity.

Note that the heat transfer inside the skin tissue is mainly the conduction of cooling capacity under CSC, and the blood perfusion can be ignored because of the short spurt duration of 50 ms. Furthermore, the heat transfer at the domain boundaries inside the skin can be ignored, since the computation lengths in the axial and radial direction are more than five times as large as the cooling depth (Z_t,d ~ 300 µm) and cooling radius (r_s ~ 5 mm at clinical spray distance). Thus, the adiabatic boundary condition was implemented for the other sides of the skin tissue domain shown in Figure 2. The initial temperature of the skin tissue and the air in the spray domain were assumed to be 37 °C and 25 °C, respectively. The initial temperatures of the cryogens were assumed to be their boiling points (T_b), as shown in Table 1.

The flowchart of simulation was drawn in Figure 3, wherein the complete simulation system consists of two parts, the spray field and the phase-change heat transfer inside the skin tissue. The solution procedure includes the following steps:

(1)

Calculate the initial spray parameters (droplet initial temperature, velocity, diameter and distribution) via the 1D homogenous flow model [19] and the piecewise linear model [22];

(2)

Compute the droplet velocity, diameter and temperature via the droplet movement, breakup and evaporation models [30,32,33,34], under the Lagrangian framework using Equations (1)–(4);

(3)

Calculate the gas concentration, velocity and temperature distribution by solving the mass, momentum and energy governing equations (Equations (5)–(7)) under the Eulerian framework with the implicit time scheme and SIMPLE algorithm;

(4)

Calculate the heat and mass transfer between the droplet and the surrounding gas using the two-way coupling method with the source terms S_m, F and S_h;

(5)

Calculate the heat flux on the skin surface with the surface heat transfer model (Equations (10) and (11)) [10];

(6)

Compute the temperature distribution by solving the macro-scale bio-heat conduction equation (Equation (8)) [31] coupled with the latent heat (L) released from the cellular water freezing;

(7)

Compute the probability of intracellular ice formation (PIIF) and the radius of the cylindrical extra-cellular space (R_ec) by solving the micro-scale mass transfer equation (Equation (9)) [35].

Since the number of spray droplets decreases along the radial direction, the uniform radiant-type tetrahedral mesh was used to discretize the computational domain. Two sets of grid sizes were employed for the spray field and skin tissue because of the large differences in domain size. The mesh sizes with an azimuth angle (Δ$\varnothing$), radial size (Δr) and cell height (Δz) of 9° × 0.8 mm × 2 mm, 6° × 0.8 mm × 1 mm and 3° × 0.4 mm × 1 mm, respectively, were used to test the mesh sensitivity for the spray domain. Figure 4 illustrates that no apparent distinction in droplet diameter distribution can be seen when further refining the grid from 6° × 0.8 mm × 1 mm to 3° × 0.4 mm × 1 mm. Thus, the 6° × 0.8 mm × 1 mm grid was selected for the spray domain.

As shown in Figure 5, no significant distinction in tissue temperature can be seen with cell size smaller than 6° × 80 µm × 20 µm. Thus, the cell size for skin tissue can be determined to be 6° × 80 µm × 20 µm. The cell number is 5,625,000 (60 × 625 × 150). The cell size of 6° × 80 µm × 20 µm is small enough to discretize the computational domain at different skin layers.

Figure 6 shows the comparison of axial temperatures in the skin tissue with uniform (6° × 80 µm × 20 µm) and nonuniform mesh systems, which were comprised of four different grid subsystems with identical Δ$\varnothing$ and Δr, but different cell heights Δz of 3 µm, 5 µm, 10 µm and 20 µm for the epidermis, dermis without melanin, dermis with melanin, and dermis without melanin layers, respectively. Figure 6 indicates that no significant distinction in axial tissue temperature was found using the uniform versus nonuniform mesh systems. To save computational load, the uniform mesh system was employed for the current simulation.

R134a spray cooling was first simulated to validate the accuracy of the present model. A straight tube nozzle (inner diameter, 1.0 mm and length, 60 mm) was implemented to produce spray droplets. The simulated spray pattern (Figure 7a) was drawn in a 2D diagram to validate the accuracy of the present algorithm.

Figure 7 illustrates that good agreement in the spray pattern was found between the experiment of Zhou et al. [40] and the simulated results of the present model. The simulated droplet temperature and diameter, and the surface temperature of the skin tissue, match well with those in the experiment, proving the accuracy of the present model (Figure 8).

As shown in Figure 9, most of relative errors in the simulated results were less than 6%, which further validates the accuracy of the present model. The maximum relative errors for simulated droplet temperature and skin surface temperature were 11.9% and 10.6%. A slightly larger maximum relative error (19.1%) in droplet diameter occurred near the nozzle exit, mainly because of the measurement uncertainty of PDPA due to the disturbance of the incompletely broken liquid block in the thick spray near the nozzle exit, and the absence of an accurate flashing atomization model for cryogens with low boiling points and high volatilities.

3. Results and Discussions

The spray fields with R134a, R404A and R32 were compared in this study. Afterward, the spray distances with three cryogens were optimized, and the potential application of the new cryogen R32 in clinics was investigated. The spray durations to prevent cold injury were determined using the multi-scale heat transfer model for skin tissue frostbite, and the intervals between the CSC results with the three cryogens and laser irradiation were optimized by comprehensively considering cooling performance, cooling depth, cooling uniformity and risk of cold injury.

3.1. Comparison of Spray Characteristics Using the Three Cryogens

As depicted in Figure 10, R32 produced the lowest droplet temperature because of its lowest boiling point and highest latent heat. However, little differences were found in the spray pattern, which indicates that the cooling areas were similar for the same spray distances achieved by the three cryogens. The axial and radial distributions of spray are displayed in a 2D cloud pattern. As shown in Figure 11, core regions occurred in the R134a, R404A and R32 sprays, in which the droplet temperature was the lowest. The droplet temperature increases when axial and radial distance increase, owing to the enhanced convective heat transfer from warm air.

The spray characteristics (droplet temperature, Sauter mean diameter and axial velocity) obtained with R134a, R404A and R32 are drawn in Figure 12. The droplet temperature and Sauter mean diameter decreased as the axial distance increased, when 15 mm < Z ≤ 70 mm, due to droplet evaporation. Afterward, the droplet temperature and Sauter mean diameter reached their minimum values and remained stable. In the initial evaporation period, the convective heat transferred from the ambient gas cannot supply enough heat induced by the extremely fast evaporation, leading to the rapidly decreasing droplet temperature because of the energy absorption from the remaining droplet liquid. Then, the droplet evaporation rate becomes smaller and smaller as the droplet temperature decreases, due to the decreasing molecular diffusion rate at the cryogen–air interface. Meanwhile, the convective heat transfer is enhanced because of the large temperature difference between the droplet and ambient gas, leading to the gradual reduction in droplet temperature. Finally, the droplet temperature reached its minimum value when the evaporation latent heat was counteracted by the convective heat transfer from the gas. Furthermore, compared with the R134a and R404A, a larger temperature reduction occurred with the R32 spray because of the more significant evaporation effect.

3.2. Optimization of Spray Distance

Figure 13 demonstrates that the maximum heat flux (q_max) on the skin surface obtained using the three cryogens increased due to the droplet temperature reduction (Figure 12a). Afterwards, q_max reached the maximum value, and then decreased slowly because of the reduction in droplet impact kinetic energy caused by the decrease in droplet diameter and velocity. The q_max values of R134a, R404A and R32 achieved their respective peaks of 363.5, 459.9 and 603.6 kW∙m⁻² when the spray distances were 66.0, 43.1 and 22.5 mm, respectively. R32 had the best cooling performance, whereas R134a had the lowest cooling capacity. The maximum surface heat flux for the R32 spray cooling was 1.66 and 1.31 times as large as that for R134 and R404A, respectively. The smallest optimum spray distance increased the cooling concentration, which was beneficial in order to accurately control the cooling for target.

As shown in Figure 14, the spatial and temporal surface heat flux distributions with R134a, R404A and R32 under their optimum spray distances exhibited considerable non-uniformity. The heat flux at the spray center initially increased, and then decreased with an increase in time, a result that was caused by a transient spray with a short spurt duration. The surface heat flux decreased gradually as the radial distance increased. Compared with R134a and R404A, R32 produced a larger surface heat flux and a smaller cooling area, thereby enhancing the cooling selectivity of the skin tissue. However, the large spatial and temporal non-uniformity of surface heat flux was detrimental to the uniform cooling of the skin tissue. Overall, R32 showed great clinical potential in CSC for the laser treatment of Ota’s nevus due to its large cooling capability and low GWP.

3.3. Temperature Distribution Inside Skin Tissue during CSC

The temperature distributions in the central axis section inside the skin tissue during CSC under the optimum spray distance are shown in Figure 15. The surface temperature of the skin tissue decreased to the minimum value with R32 spray cooling, i.e., R32 produced the best cooling effect. However, the cooling depths achieved with R134a, R404A and R32 were almost identical (300 μm) because of the small thermal conductivity of skin tissue.

Further detail on the cooling effect with the three cryogens can be found in Figure 16. The axial temperature increased with an increasing tissue depth. Compared with R134a and R404A, R32 spray cooling achieved the lowest tissue temperature, which allows for greater laser fluence and will improve the therapeutic effect. The tissue temperature further decreased when the cooling duration was extended.

3.4. Optimization of the Interval between CSC and Laser Irradiation

In addition to cooling duration, cooling depth and uniformity are also affected by the time interval (t) between CSC and laser irradiation. The best cooling effect can be obtained by simultaneously considering cooling depth and uniformity with cryogen R32 under the optimum spray distance. However, skin tissue suffers from cold injury with excessive cooling. Thus, t should be determined by considering cooling depth, cooling uniformity, and risk of cold injury, which can be predicted by the multi-scale model.

Given that the diameter of the clinical laser spot is within 4 mm, the radial surface temperature difference (∆T = T_{r = 2 mm} − T_{r = 0 mm}) is defined to represent the cooling uniformity during CSC. As shown in Figure 17a, ∆T increased rapidly to the maximum value at t = 14 and 35 ms with R404A, R32, and R134a. ∆T almost decreased to its minimum value (6.1 °C with R134a, and 7.8 °C with R404A and R32) at t = 90 ms under CSC with the three cryogens. Compared with the cooling non-uniformity under R404A and R32, that under the R134a spray cooling was the smallest. According to the histological anatomy of Ota’s nevus [37], melanin lesions are buried in the deep dermis (500 μm < Z_t ≤ 2000 μm). The temperature distribution in the dermis with a melanin layer (500 μm < Z_t ≤ 2000 μm) should not be affected by CSC. Therefore, the cooling depth (Z_t,d) should be controlled to less than 500 μm.

As shown in Figure 17b, the cooling depth increased as the cooling time increased. The Z_t,d values under CSC with the three cryogens were almost identical. t was determined to be 90–162 ms by considering cooling uniformity and cooling depth. The cooling depth with the function of cooling time was obtained with the least-squares method, as follows:

$$Z_{t,}{}_d=-752.56{\mathrm{e}}^{-t/181.81}+808.69$$

(12)

The dynamic radius of the extra-cellular space (R_ec) and the PIIF at the interface (Z_t = 60 μm) between the epidermis and dermis are drawn in Figure 18 to evaluate the cold injury risk of skin tissue under CSC. As mentioned previously, cold injury occurs when the PIIF is 1 or R_ec reaches its maximum value (R_ec,max). As illustrated in Figure 18a, the R_ec increased exponentially as the cooling duration increased. The R_ec increased to its maximum values when t = 2.3, 1.4, and 1.1 s for R134a, R404A, and R32 spray cooling, respectively. t should not exceed these durations for R134a, R404A, or R32 spray cooling to prevent the cold injury of skin tissue. Thus, a spray duration of 50 ms for clinical use can ensure the safety of CSC application. Under the spray duration of 50 ms and the optimum spray distance of 22.5 mm, the radial surface temperature difference and cooling depth were determined to be 7.8 °C and 354.8 μm, when the interval between CSC and laser irradiation was 40 ms.

The PIIF values under CSC with the three cryogens are depicted in Figure 18b. The variations in the PIIFs were similar to those of R_ec. However, the PIIFs were less than 6% owing to the low thermal conductivity of skin tissue. Although the surface temperature decreased rapidly upon the subcooling droplet’s impact on the skin surface, the temperature inside the skin tissue decreased gradually, suggesting that the cold injury under CSC was dominated by cell dehydration. Moreover, the PIIFs when cold injury occurred were 5.6%, 1.4% and 1.3% for R134a, R404A, and R32 spray cooling, respectively, indicating that CSC entailed a very low possibility of IIF. According to cellular frostbite theory, PIIF has a positive correlation with the effective cell superficial area (A) for water transport and integral time (t). The high cell dehydration rates induced by R404A and R32 spray cooling, with good cooling capability (i.e., lower cell temperature), produced a small A, leading to a small PIIF for skin frostbite with the R404A and R32 sprays. In consideration of the cooling effect, depth and uniformity, and the risk of cold injury, t was optimized as 90–162 ms for the three cryogens.

4. Conclusions

In this work, an accurate, efficient, complete simulation system that considers the entire cryogen spray cooling of gas–liquid flow in the capillary tube, the atomization, evaporation and breakup of spray droplets, and the phase-change heat transfer (cold injury) inside the bio-tissue, were established and validated. The spray characteristics in terms of droplet temperature and diameter, and the heat transfer performance on the skin surface during CSC with R134a, R404A, and R32, were compared. The temperature distributions and cold injury inside the skin tissue were simulated using the multi-scale heat transfer model based on cellular frostbite theory to evaluate the cooling effect during CSC quantitatively.

During spraying, the droplet temperature and Sauter mean diameter decreased exponentially along the axial distance because of droplet evaporation. R32 produced the lowest droplet temperature, the smallest droplet diameter, and the lowest skin tissue temperature, owing to its lowest boiling point and highest latent heat among the three cryogens. This result indicates that R32 had the best cooling capability. The maximum surface heat fluxes of R134a, R404A and R32 were respectively optimized at 363.5, 459.9 and 603.6 kW∙m⁻², at spray distances of 66.0, 43.1 and 22.5 mm, respectively. The maximum surface heat flux under R32 spray cooling was 1.66 and 1.31 times as large as the values under R134a and R404A. Moreover, the smallest optimum spray distance of R32 spray cooling increased the cooling concentration, which was beneficial for the accurate cooling of skin. Given its low GWP, the environmentally friendly R32 has excellent potential as a clinical CSC technique in the laser treatment of Ota’s nevus, in terms of cooling capability and control of the cooling area.

The risk of cold injury in skin tissue was evaluated by considering cell dehydration and IIF simultaneously. The low probability (<6%) of IIF under CSC indicates that the cold injury in the skin tissue was dominated by cell dehydration, because of the low heat conductivity of skin tissue. To prevent cold injury, the intervals between CSC and laser irradiation should not exceed 2.3, 1.4 and 1.1 s for R134a, R404A and R32 spray cooling, respectively. Thus, a spray duration of 50 ms in clinical use can ensure the safety of CSC application. Considering the cooling effect, depth and uniformity, and the risk of cold injury, the interval between CSC with the three cryogens and laser irradiation was determined to be 90–162 ms.