Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a

Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a

Applied Thermal Engineering xxx (2014) 1e12 Contents lists available at ScienceDirect Applied Thermal Engineering journal homepage: www.elsevier.com...

2MB Sizes 0 Downloads 50 Views

Applied Thermal Engineering xxx (2014) 1e12

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a D. Li a, B. Chen a, *, W.J. Wu a, G.X. Wang a, b, Y.L. He a a b

State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Department of Mechanical Engineering, University of Akron, Akron, OH 44325-3903, USA

g r a p h i c a l a b s t r a c t

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 January 2014 Accepted 16 March 2014 Available online xxx

In laser dermatologic surgery, cryogen spray cooling (CSC) has been proved to be an efficient cooling technique to avoid thermal damage from skin burning due to the light energy absorption by melanin in epidermis. R134a is now the only cooling agent in commercial laser and has been proved to be effective for light pigmented skin. R407c and R404a could provide better cooling effect for darkly pigmented skin than R134a because both of them have much lower boiling point. In order to investigate the potential cold injury mechanism prior to the further clinic use, this paper presents a multi-scale model to simulate the cooling process of the skin and estimate the potential cold injury. In the model, the skin tissue is treated as multi-layered geometry and the heat transfer within this multi-layer skin is described by a macro-scale bio-heat transfer model. A general dynamic relation is introduced on the surface of skin to quantify the convective cooling of CSC with various cooling agents. Meanwhile, the micro-scale mass transfer and the ice formation in cell during the cooling are evaluated in a Krogh unit. The cold injury is recognized once the cell is dehydrated or the ice formed intracellularly. The results show that the surface cooling effect of spray cooling is well related with the boiling point of cryogen. Much lower surface and inner skin temperature will be achieved by using cryogens with lower boiling poring, e.g. R404a and R407c, which is benefit to thermal damage protection for darkly pigmented skin. Recognized as cell dehydration, the spray durations to cause cold injury are 3.3 s, 2.2 s and 1.9 s for R134a, R407c and R404a, which proved that three cooling agents are all safe for epidermis protection in CSC with spurt duration of tens of milliseconds in clinic. Ó 2014 Elsevier Ltd. All rights reserved.

Keywords: Cryogen spray cooling Laser dermatology Cold injury Micro-scale heat and mass transfer

1. Introduction * Corresponding author. E-mail address: [email protected] (B. Chen).

Port Wine Stain (PWS) birthmark is a congenital and progressive vascular malformation of the capillaries in dermis, which occurs in

http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034 1359-4311/Ó 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

2

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

Nomenclature A a C c h L p PIIF R r S s T t V z

cross-sectional surface area for cellular water transport (m2) control parameter in energy equation molar concentration (mole/m3) specific heat (J/kg K) heat transfer coefficient (w/m2 K) latent heat (kg/kJ) cellular permeability constant probability of intracellular ice formation radial coordinate (m) radius of the cylindrical extra-cellular space in Krogh unit (m) surface area of the Krogh cylinder (m2) phase-change interface (m) temperature ( C) time (s) volume of Krogh unit(m3) axial coordinate (m)

Greek symbols osmotic coefficient density (kg/m3) thermal conductivity (w/m K) cell type dependent constants contributing to the thermodynamic of crystal nucleation on cell membrane surface

s r l k

approximately 0.3% of children. Most PWSs are pale pink patches at birth that progressively darken and thicken with age. This often leads to increased cosmetic disfigurement and psychological distress, prompting many patients and their families to seek treatment as early as possible. Nowadays, laser has been the best choice for the treatment of PWS. Based on the principle of selective photothermolysis [1], light (laser) energy will be preferably absorbed by hemoglobin in PWS blood vessels causing irreversible thermal damage, thrombosis, and, eventually, permanent photocoagulation in the PWS vessels. Unfortunately, a significant amount of energy in laser beam will be absorbed by melanin in the basal layer of the epidermis, leading to unwanted thermal damage and successive skin dyspigmentation or hypertrophic scarring [1e3]. Cryogen spray cooling (CSC) has been developed to protect the epidermis from non-specific heating by pre-cooling the skin prior to laser irradiation [4e8]. R134a, R407c and R404a are non-toxic and environment friendly candidates for the spray cooling, among which R134a (boiling point is 26.1  C at 1 atm) has been approved by U.S. Food and Drug Administration for the use in clinical laser dermatology. Dai et al. [9] clarified that environmentfriendly cryogens R404a and R407c (boiling points are 43.6  C and 46.5  C at 1 atm) can improve the efficiency in protecting dark human skin from thermal injury since their boiling points are much lower. During CSC, the skin is cooled by continuous impact of atomized cryogen droplets at fairly low temperature (e.g. below 70  C) [10] due to evaporation in fly. Such low temperature is potential to cause cold injury in skin tissue. Several studies have been conducted to investigate the potential cold injury during skin cooling [11e14]. In their work, the potential cold injury is determined by the lethal temperature. However, the cellular-level biophysical events which occur in tissue as it freezes (i.e., cell dehydration and intracellular ice formation) have long been recognized as a

U h s n Dx

cell type dependent constants contributing to the kinetic components of crystal nucleation on cell membrane surface intracellular solution viscosity normalized time ionic dissociation coefficient dimension of Krogh unit

Subscripts and superscripts a air c cellular csc cryogen spray cooling d dermis e epidermis f frozen i ice iso isotonic conditions max maximum value ns non-solvent o initial condition ph phase change pf partial frozen u unfrozen w water * normalized

mechanistic way to explain cell injury as a result of the freezing process. Once the macro-scale skin-tissue experiences cooling and the temperature decreases below the phase change point, the micro-scale freezing will initiate within the cells. The states of the cell can be clarified as unfrozen, partial frozen and complete frozen stages. As freezing began, the cell would experience a partial frozen state before getting complete frozen: ice firstly formed extracellularly and the water transported form intracellular to extracellular space as a result of concentration difference. If the cooling rate is high enough, the water in cell will not be able to escape to extracellular space, and ice will form intracellularly. As a result, the cell will die from intracellular ice formation (IIF). Otherwise if the cooling rate is slow, the water will transport from intracellular to extracellular space till the cell dehydration. Consequently, the cell will die from concentration injury. Therefore, the dynamics of both the thermal conditions inside macro-scale frozen tissue and the micro-scale cellular-level biophysical events must be evaluated simultaneously to understand the mechanism of tissue freezing causing by cryogen spray cooling. The mathematical description of above micro-scale freezing process originated from Rubinsky [15]. They represented the liver by a large number of identical tissue/vascular structure known as Krogh unit, which consists of a sinusoid (vascular space) surrounded by hepatocytes (cellular space). In their work, the governing equations for mass transport and intracellular ice formation during freezing process were derived for an individual Krogh unit. The results show that their micro-scale freezing model is capable of qualitatively predicting the amount and extent of intracellular ice formation and dehydration of the tissue under arbitrary cooling conditions with the known thermal history of the tissue. The micro-scale freezing model in Krogh unit, which has not been used in the skin tissue before, can also be adopted to quantify the potential cold injury during CSC.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

To provide thermal history within skin tissue, the Pennes bioheat transfer model can be employed to describe the macro-scale heat transfer within skin tissue during cryogen spray cooling. A convective thermal boundary condition on the skin surface is usually adapted to depict the surface cooling effect in the model. One challenging issue is how to determine the surface heat transfer coefficient induced by spray cooling. To date, most studies implemented a constant convective heat transfer coefficient throughout the entire spray cooling process over the entire spray area [12,16e 18]. Nevertheless, experimental results have found that heat transfer coefficient shows a strong dynamic variation during and after the spray [19] with large fluctuation over the edge of the spray area [20]. In our previous work [13,21], we pursued a way to simulate the cooling process with realistic boundary condition by presenting a simplified correlation of convection heat transfer coefficient from experimental data. The simulation results match the experimental data qualitatively, but difference still exists, which may attribute to the difference in thermal properties between experimental sample and human tissue. For this reason, the corrections on the correlation of surface heat transfer coefficient are needed to match the experiment data well. In this paper, a theoretical multi-scale model will be established to investigate the cooling effect of cryogen spray cooling with various cryogens and potential cold injury by coupling macro-scale bio-heat transfer model and micro-scale freezing model in a skin tissue. The main objectives of this study are: 1) to find a general model for the cooling effect of spray cooling; 2) to describe the freezing process and examine the possible cold injury in tissue during the spray cooling with various candidate cryogens; 3) to propose a criterion of spray duration to protect skin from potential cold injury. For the first purpose, Newton’s cooling law was used as boundary condition in heat conduction model, in which surface heat transfer coefficient and cryogen temperature need to be specified for different cryogens. A general correlation of surface heat transfer coefficient obtained from experimental data will be employed as we did before [13,21]. Different from our previous work, the maximum value of heat transfer coefficient is a self-adapted parameter to make skin surface temperature stabile at the corresponding boiling point of the given cryogen after liquid cryogen pool is formed on the skin surface. Individual cryogen film temperature on the skin surface will be determined from recently experimental results by Zhou et al. [10]. For the second purpose, a set of micro-scale mass transfer equations are adopted in each Krogh unit to predict the response of tissue cell to cooling and intracellular ice formation model is used to check the potential of cold injury during spray. Finally, some suggestions on spray duration will be presented to avoid the cold injury within skin tissue.

and the micro-scale biophysical response at various locations within the tissue domain during a freezing process, simultaneously. Cryogen spray cooling is a complex process to be simulated. As a preliminary study, we assume a steady jet flow of cryogen spray impingement on the skin surface during the spurt duration. The influences of instant conditions of spray impingement, spatial droplet distribution, film formation are not considered. Geometry symmetric distributions of droplet and cryogen residual layer are assumed on the skin surface. As a result, the computation can be conducted in a two-dimensional axisymmetric domain. Also, properties of the fluid and skin in the cooling process are assumed constant. As illustrated in Fig. 1, the skin tissue was divided into two parallel layers: epidermis and dermis. It is assumed that phase change (freezing) and cold injury only occurs in dermis because the water content in the epidermis is quite low. For this macro-scale process, the freezing of tissue were clarified as unfrozen, partial frozen and frozen stages, and the governing equations of heat transfer are specified for each skin layer and each freezing stage. Robin boundary condition is used in the heat conduction equation to simulate the cooling effect of cryogen spray cooling. Heat convection coefficient hcsc is calculated from general correlation proposed by our previous work [13,21], and the cryogen film temperature is obtained from recently experiment result from Zhou et al. [10]. In Fig. 1 we show a typical Krogh unit composed of a cylindrical extra-cellular compartment and surrounding intra-cellular space, which was used to model the micro-scale physical events during freezing. The simulation starts from unfrozen stage when temperature is higher than that of phase change. Water transportation between intracellular and extracellular space only happens in partial frozen cell. Complete frozen stage (cold injury) would be recognized in cells where mass transport terminates, which could be either intracellular ice formed (probability of IIF exceeds 1) or cell dehydrate (cylinder radius exceeds maximum value). 2.1. Macro-scale heat transfer model and corresponding boundary condition of cryogen spray cooling The macro-scale heat transfer during cooling can be described by the following bio-heat transfer equations for epidermis and dermis, respectively. In this study, blood perfusion and the metabolic heat generation are not considered since the spurt duration is

2. Multi-scale model The multi-scale model consists of two components: 1) macroscale heat transfer model to calculate the heat conduction in twolayer human skin during cryogen spray cooling; 2) micro-scale model carried out in Krogh unit to predict the cellular volume shrinking as water transport from intracellular space to extracellular space and predict intracellular ice formation (IIF) during cooling. The coupling of the micro-scale biophysical phenomena into the calculation of phase change is accomplished by creating a latent heat function in the macro-scale heat transfer model. The latent heat released from cellular water freezing for each control volume is determined by the cellular water transport and intracellular ice formation processes in a micro-scale Krogh unit. The major advantage of the coupled thermal/biophysical model is its unique ability to predict both the macro-scale thermal response

3

Fig. 1. Schematic of two-layer skin model and Krogh unit.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

4

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

short (less than 5 s). Thus, the Pennes bio heat transfer equation can be reduced to heat conduction equations as below:

re ce

h*



vT ¼ le V2 T vt

R* ; s



(1)

8 < h* ðsÞ hcsc ðR; tÞ o ¼ ¼ : 51  R* =½3h* ðsÞ ho;max 8 > h* ðsÞ > > < o

o

correlation between normalized convection heat transfer coefficient (h* ¼ hcsc/ho,max) for the spray distance ranging from 30 to 60 mm and the non-dimensional time (s ¼ t/tmax) and nondimensional space (R* ¼ R/Rspray) has been proposed as the following [13,21]:

0  R*  0:4

s  1:0

0:4 < R*  1:0

  hcsc ðR; tÞ   * ¼ h* R* ; s ¼ 5R  s  1 > ho;max > * >  ½0:09ðs  1Þ  h*o ðsÞ : ho ðsÞ þ ð4  sÞ h* ðR* ; sÞ ¼

rd;j cd;j

hcsc ðR; tÞ ¼ h*o ðsÞ ho;max

0  R*  0:2ðs þ 1Þ 



0:2 s þ 1 <

R*

 1:0

s  4:0

0  R*  1:0

 vT 2pr rL dr ¼ ld;j V2 T þ a ð1  PIIFðT; tÞÞ vt Dx2 dt  ðVc  Vns ÞrL dPIIFðT; tÞ þ dt Dx3

with h*o ðsÞ as the normalized convection heat transfer coefficient at the center of spray:

(2)

where r is the density, c is the specific heat, l is the thermal conductivity, T is the temperature, and t is the time. Subscripts e in Eq. (1) and d in Eq. (2) represent epidermal layer and dermal layer, respectively. In epidermis, only heat conduction is considered. However, all the three stages of the tissue freezing should be considered in dermis. In Eq. (2), subscripts j ¼ u, f and pf, which denotes the unfrozen, complete frozen and partial frozen stages, respectively. The last term on the right hand of Eq. (2) is used to describe the latent heat released during tissue freezing within cryogen spray cooling: (1  PIIF(T,t))2prrL/Dx2dr/dt represents the contribution of the intracellular water which freezes after being transported from intracellular to extracellular space, while (Vc  Vns)rL/Dx3dPIIF(T,t)/dt represents the contribution of water freezes intracellularly. Vc and Vns are the volumes of cell and nonsolvent water in cells, respectively. PIIF is the probability of intracellular ice formation and L the latent heat of water. Dx is the dimension of the Krogh unit and r is the radius of the cylindrical extra-cellular space in Krogh unit. a is a control parameter, which is 1 if the tissue is under partial frozen stage and 0 if the it is under unfrozen or complete frozen stage. On the skin surface (z ¼ 0), the standard convection condition was used to describe the cooling effect of CSC:

 8 vT  > ¼ hcsc ½Tð0; tÞ  Tc  during CSC > < l vz  z¼0  >  > : l vT  ¼ ha ½Tð0; tÞ  Ta  after CSC vz z¼0

(4)

1:0 < s < 4:0

(3)

where Tc is cryogen film temperature, Ta is the air temperature which is set to be 25  C. ha and hcsc are heat convection coefficient during and after CSC, respectively. In most of existing models of CSC for laser dermatology, Newton’s law of cooling was employed and the value of a heat convection coefficient hcsc was usually assumed constant during the entire duration of spray [12,16e18]. Recent experiments [19,20] found that hcsc is a strong function of both time and space. Based on the experimental data, a generalized

h*o ðsÞ ¼

8 > > <

s

hcsc ð0; sÞ 1:0  0:35ðs  1Þ ¼ 0:3  0:02ðs  3Þ > ho;max > : 0:2  0:0125ðs  8Þ

s  1:0 1:0 < s  3:0 3:0 < s  8:0 8:0 < s (5)

where hcsc is surface heat convection coefficient and ho,max is the maximum value of heat transfer coefficient during spray, t is the spray duration and tmax is the time when ho,max occurs. It is noticed that above correlation includes two free parameters, i.e., tmax and ho,max, which have to be experimentally determined and are available only for R134a [19,20]. In this paper, we assume that this correlation is applicable to other cryogens like R407c and R404a, as suggested from our recent experiments [10]. In the model, ho,max is treated as a free parameter. For a given CSC condition, it will selfadapt until the skin surface temperature is stabilized at the corresponding boiling point of the given cryogen after liquid cryogen pool is formed on the skin surface [19]. After freezing began in the tissue, boundary conditions between the unfrozen, partially frozen and completely frozen regions in the tissue must be specified. Between the unfrozen and partially frozen

Table 1 Geometric, thermal, and mass transfer properties of the tissue [16,22]. Properties

State

Epidermis

Dermis

D (mm) l (W/(m K))

/ Unfrozen Frozen Unfrozen Frozen Unfrozen Frozen / /

50 0.209 0.209 1000 1000 3530 3530 / /

550 0.552 2.25 998 921 4230 1230 335 22 4.5 0.3(pr2oDx) 1 2 0.154 5  107

r (kg/m3) c (J/(kg K) L (kJ/kg) Dx (mm) ro (mm) Vns (mm3)

s n

Cco (moles) p (m/mole s))

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

regions, the boundary condition is the initial freezing of water in extracellular part of the Krogh unit:

8 > pro2 vsðtÞ vT vT > < lph  lu ¼ rL vz vz ðDxÞ2 vt > > : TðsðtÞÞ ¼ T ph

(6)

where s is the position of the phase change interface between the unfrozen and partially frozen regions as a function of time t. ro represents the initial radius of the Krogh cylinder, and subscript ph represent phase change. The boundary condition between the partially and completely frozen regions can be either geometric or thermodynamic. The geometric condition is when all unbound cellular water is transported into the extra-cellular space and thus the Krogh cylinder reaches a maximum radius, when the tissue is all frozen. The thermodynamic condition is met when the probability of IIF gets 1, at which all unfrozen water in the cells changes phase. The frozen of the tissue will complete if either condition is satisfied. Geometric, thermal and mass transfer properties of each skin layer in the unfrozen and frozen regions are given in Table 1. The skin initial temperature is assumed to be 37  C. In partial frozen regions, the thermal properties are the function of the amount of ice and water, which can be expressed by the radius of Krogh cylinder,

Yðrðx; tÞÞ ¼

ððDxÞ3  prðx; tÞ2 DxÞYw þ prðx; tÞ2 Yi Dx D x3

(7)

5

  Tphw  TðR; z; tÞ drðR; z; tÞ Vco  Vns ¼ p þ Cco dt Vc ðR; z; tÞ  Vnx 1:86sy

(12)

where Tph is the freezing point of water (0  C); s and n are osmotic coefficient and ionic dissociation coefficient, respectively; Vns is the volume of non-solvent water in cells. p is the cellular permeability constant. Sv is the surface area of the Krogh cylinder. Cc is the cellular space molar concentration. The probability of intracellular ice formation is determined using the model proposed by Toner [23]:

2 6 PIIFðT; tÞ ¼ 1  exp4 

ZT

AUo

Tseed

T Tfo

!1=2

2 . 1=4 3 3 k 6 iso Tf Tfo 7 7  exp4 5dt 5

2 T  Tf T 3

ho A h Ao

(13)

where Tseed is the temperature when ice is seeded in extracellular space; Tf is the equilibrium phase change temperature of the intracellular water, A is the surface area of the cell available for water transport. k and U are cell type dependent constants contributing to the thermodynamic and kinetic components of crystal nucleation on cell membrane surface, respectively. The subscript iso denotes the value of that parameter under isotonic conditions. h is the intracellular solution viscosity, which is approximated using the following expression [22]:



1:64

where Y denotes to the thermal properties of thermal conductivity, density or heat capacity, and the subscript w and i denotes the water and ice, respectively.

h ¼ 0:139  103

2.2. Micro-scale model for cell freezing

2.3. Numerical strategy

The Krogh model is defined by cylinder (extracellular space) and its surrounding (cellular space) tissue as illustrated in Fig. 1. The geometric of Krogh unit defined below: The volumes of Krogh unit and cellular space are:

The above mathematical equations constitute a multi-region two-dimensional transient heat conduction problem for cryogen spray cooling. It is a multi-scale model, where micro-scale freezing model is used to quantify the mass transport and intracellular ice formation during cell freezing and macro-scale model is used to describe the heat transfer in tissue skin during CSC. As mentioned in the beginning of Section 2, the potential cold injury during CSC can be predicted mathematically by coupling the macro- and micro-scale model. The problem is numerically solved in a rectangular computational domain with R  z ¼ 5 mm  600 mm. Grid-independent tests are conducted on four different uniform grids: 125  15, 250  30, 500  60 and 1000  120. The solution by 500  60 grid is adopted in the following calculations because the calculated surface temperature is only about 5% deviated from that of 1000  120 grid. Similar examinations were also conducted for the time discretization and 0.1 ms is the independent time step for the first-order implicit discretization used in this study. The convergence criteria for governing energy equation is defined as that the relative error between two successive iterations is less than 106 in each time step. The flow char of simulation was indicated in Fig. 2. The solution procedure includes the following steps:

V ¼ ðDxÞ3

(8)

Vc ¼ V  pr 2 Dx

(9)

The initial and the maximum value of the radius of the Krogh unit can be represented as below,

 rðx; 0Þ ¼ ro and rmax ¼

Dx3  Vns pDx

1=2 (10)

where subscripts o and ns represent initial value and non-solvent components within the cell space, respectively. Cellular compartment volume (Vc) change due to flow of water from the intracellular space to extra-cellular space during freezing can be written as [22]:

  Tph  TðR; z; tÞ dVc ðR; z; tÞ Vco  Vns ¼ pSv þ Cco dt 1:86sy Vc ðR; z; tÞ  Vnx

(11)

This leads to an equation describing the change of the radius of the extra-cellular compartment, or Krogh cylinder, as a function of time and position [22],

T 1 225

(14)

1) To solve the macro-scale heat transfer Eqs. (1) and (2) together with the convective boundary condition Eqs. (4) and (5) with an implicit time scheme by the tri-diagonal matrix algorithm (TDMA) with block correction.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

6

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

Fig. 2. Flow chart of our calculation.

2) To track the phase change front by solving a Stephen problem (Eq. (6)) if the dermal temperature decreases below the phase change point. 3) To quantify the mass transport from intra-cellular space to extra-cellular space and the probability of intracellular ice formation to each control volume by solving micro-scale freezing model Eqs. (12) and (13) based on the known node temperature if the dermal temperature decreases below the phase change point. 4) To couple the macro- and micro-scale model by incorporating the latent heat releases from mass transport and intracellular ice

formation into the energy Eqs. (1) and (2) and update the nodal temperature with the updated thermal properties (Eq. (7)).

2.4. Model validation To validate our model, the movement of phase-change interface (interface between unfrozen and partial frozen regions) is checked by solving a 1-D water-ice phase change process in half-space. The water at initial uniform temperature (To ¼ 0  C) is confined to a semi-infinite space z > 0. At time t ¼ 0, the temperature Tz ¼ o at the

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

7

Fig. 5. Dynamic variation of skin surface temperature as a function of spray duration with different cryogen: R134a, R404a and R407c, respectively. Fig. 3. Comparison of numerical simulation and analytical solutions for ice front propagation.

boundary surface (z ¼ 0) is lowered to 1  C. As such, solidification starts at the surface z ¼ 0 and the ice front moves in the positive z direction. Fig. 3 displays the comparison between numerical results for ice front location z versus time t by our model with the analytical solution given by Ozisik [24]. As we can see from the figure, the ice front propagation obtained by numerical solution agrees well with the analytical solution, which verified the validation of our model.

3. Results and discussion 3.1. Dynamic temperature variation at skin surface and epidermis/ dermis interface during cryogen spray cooling with various cryogens Fig. 4 compared the dynamic temperature variation in spray center at cooling surface using variable heat transfer coefficient model with traditional constant heat transfer coefficient (4 kW/ m2 K) [16] and the experiment measurement of R134a spray cooling under the clinic condition [19,20]. In the experiment, polymethyl methacrylate resin was selected as the sample surface

Fig. 4. Comparison of the dynamic temperature variation in spray center at skin surface with constant and variable heat transfer coefficient.

because of its thermal conductivity is similar to that of human skin. Surface temperature was measured by a fast response type-K thermocouple with 50 mm bead diameter and heat transfer coefficient was calculated by an inverse heat conduction problem algorithm. With constant heat transfer coefficient, a continuous decreasing of the skin surface temperature is predicted, which overestimated the cooling capacity of the cryogen spray cooling. With the variable heat transfer coefficient, however, the skin surface temperature decreased very fast to a minimum value as the cryogen droplets reach skin surface. Then the surface temperature slowly increases to an almost constant value until the end of the spray as a result of the formation of liquid cryogen pool. This interesting dynamic behavior is in consistent with the experimental observation. In the initial stage of the spray, the simulated temperature decreases faster, which we believe is due to the indirect measurement of the surface temperature under a thin layer of aluminum foil in the experiment [19]. As the spray continues, the predicted temperature by using variable hcsc agrees very well with the experiment measurement. Therefore, the dynamic transient cooling process during CSC can be better predicted by the variable

Fig. 6. Dynamic variation of epidermis/dermis interface temperature as a function of spray duration with different cryogen: R134a, R404a and R407c, respectively.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

8

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

heat transfer coefficient than traditional constant heat transfer coefficient. Fig. 5 illustrates the dynamic variation of skin surface temperature as a function of spray duration (totally 100 ms) with different cryogens: R134a, R404a and R407c. As the spray began, cryogen droplets reach skin surface with a very low temperature which is even lower than its boiling point due to the evaporation during flight. As a result, the temperature of skin surface (z ¼ 0 mm) decreased very fast to a minimum value. Then after about 20 ms, the temperature of the skin surface will slowly increase to a constant value corresponding to the boiling points of each cryogen (R134a: 26.1  C, R404a: 46.5  C and R407c: 43.6  C), due to the formation of the cryogen pool on the skin surface. However, their difference is not enlarged as the spray duration prolong. It implies that the surface cooling effect of spray cooling is well related with the boiling point of cryogen. The cryogen with lower boiling point (R404a, R407c) would perform better cooling effect than that with higher boiling point (R134a). Fig. 6 plots the dynamic variation of epidermis/dermis interface temperature as a function of spray duration with different cryogens: R134a, R404a and R407c, respectively. When the spray began, the interface temperature decreased as the spray duration prolonged. The epidermis/dermis interface temperature was lower by cooling with lower boiling point cryogen than that with higher boiling point cryogen. Other than the surface temperature variation, the interface temperature decreased as the spray duration prolonged. The interface temperature difference between different cryogen was enlarged as the spray duration further prolonged. As shown in Figs. 5 and 6, the cryogens with different boiling point will perform distinctly different cooling effect on the epidermis. By using cryogen with lower boiling point, the temperature of tissue surface and epidermis/dermis interface is much lower than that with higher boiling point cryogen. It can be expected that the spray cooling effect will be stronger by using cryogens with lower boiling point. Secondly, the skin surface temperature will be stable on boiling point of cryogens other than further decreased as spray duration prolongs. However, the temperature of epidermis/dermis interface would keep decreasing. It suggested that the cooling effect difference between cryogens with different boiling point will increase as the spray duration prolongs. 3.2. Spatial temperature distribution inside skin tissue immediately after cryogen spray cooling with various cryogens The spatial temperature distributions in the skin after cryogen spray with R134a are displayed in Fig. 7 at three time instants during spray. It is interesting to notice that the skin surface temperature drops to about 30  C at 10 and 30 ms after the spray starts (see Fig. 7(a) and (b)). At this instant, however, the temperature within the skin including majority of the epidermal layer still remains the same initial value. After 60 and 100 ms of CSC (Fig. 7(c) and (d)), the skin surface temperature remains low (it recovers a bit due to a dynamic variation in hcsc), but the entire epidermal layer is cooled down (below 10  C). The entire skin temperature keeps decreasing if the spurt duration extended to 200 ms (see Fig. 7(e)). At this time, part of the dermis is also cooled down. Available models for cryogen spray cooling of laser dermatological surgery assume a spatially uniform value of hcsc. However, experimental measurements have clearly shown a spatial variation in both the heat flux and temperature [10,20]. We have introduced such spatial variation in present model, and the spatial variations of temperature in the skin during cryogen spray with R134a can be observed in Fig. 7. The temperature inhomogeneity over the spray area is clearly shown in above figures. Since the cooling is not uniform, a much higher minimum temperature is obtained for the

Fig. 7. Calculated temperature distributions within the skin during cryogen spray at times after 10(a), 30 ms (b), 60 ms (c), 100 (d) and 200 ms (e).

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

9

Fig. 8. Temperature distribution along the tissue depth with different spray duration by R134a.

Fig. 10. Temperature distribution along the tissue depth with different spray duration by R407c.

location farther away from the center of the spray. At the beginning of the spray, the center area drops to the subzero, while the temperature near the boundary is still high. Such inhomogeneity could give the risk to insufficient protection, which is unwanted in the clinical practice. After about 100 ms spray (see Fig. 7(d)), the cooling effect covered the entire area. After 200 ms spurt duration, the uniform cooling can be achieved over the entire spray area (see Fig. 7(e)). The variation of the temperatures distribution along the tissue depth at R ¼ 0 in Fig. 7 is summarized and presented in Fig. 8. Due to the limitation of the low thermal conductivity of human tissue and short spray duration, the temperature decrease in deep epidermis is much slower than that at the skin surface. For example, 30 ms spray duration is sufficient for the surface temperature to reach its minimum value. However, the temperature in the deep epidermis is not low enough to avoid thermal injury and longer spray duration is needed. On the other hand, for the precooling of laser treatment of PWS in dermis, the spray duration should be controlled to cool the epidermis selectively. Excessive spurt would cool the dermis instead of further cooling of the epidermis. For example, there is not much temperature difference in epidermis between 100 ms and 200 ms spray. However, the temperature in dermis is much lower with 200 ms spray than that with 100 ms spray. Take skin depth of 200 mm for example, the skin

tissue temperature is 30  C and 20  C after 100 and 200 ms spurt duration, which is not beneficial to the treatment of PWS because the final coagulation of the PWS blood vessels will be weakened by the excessive cooling in dermis. It can be concluded that epidermis cooling will get less efficient as the spray duration prolongs. As the spray duration further increase, the temperature of dermis may fall

Fig. 9. Temperature distribution along the tissue depth with different spray duration by R404a.

Fig. 11. Predicted extra-vascular radius (a) and IIF probability (b) as a function of R134a spray duration at different depth: 61 mm, 83 mm, 105 mm, 127 mm, 149 mm and 171 mm, respectively.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

10

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

below the freezing point of the water and potential cold injury could occur. Temperature distributions at the spray center with different spray duration by R407c and R404a are respectively shown in Figs. 9 and 10 as a function of tissue depth. As we can see from the two figures, the tissue temperature in both epidermis and dermis is much lower than that spray with R134a (Fig. 8) due to much stronger cooling effect of R407c and R404a. The minimum skin surface temperatures are 50  C and 55  C after R407c and R404a spray, which are much lower than that of R134a (36  C). Lower epidermis temperature will improve the laser threshold fluence and benefit the final treatment outcomes. Nevertheless, as we can see from the figures, much lower dermal temperatures are also predicted and the subzero area of the tissue extended to a depth of 115 mm and 136 mm after R407c and R404a cooling, respectively. The results indicate that cryogens with lower boiling points will cool down the epidermis more effective than R134a. Nevertheless, their potential cold injury needs to be carefully investigated. 3.3. Analysis of potential cold injury Fig. 11 exemplifies predicted extra-vascular radius of Krogh unit (Fig. 11(a)) and IIF probability (Fig. 11(b)) as a function of R134a spray duration at different tissue depths: 61 mm, 83 mm, 105 mm, 127 mm, 149 mm and 171 mm, respectively. As we can see from Fig.

Fig. 13. Predicted extra-vascular radius (a) and IIF probability (b) as a function of R404a spray duration at different depth: 61 mm, 83 mm, 105 mm, 127 mm, 149 mm and 171 mm, respectively.

Fig. 12. Predicted extra-vascular radius (a) and IIF probability (b) as a function of R407c spray duration at different depth: 61 mm, 83 mm, 105 mm, 127 mm, 149 mm and 171 mm, respectively.

11(a), the temperature in the control volume is below the 0  C after about 175 ms spray duration at the depth of 61 mm, which is the depth of first Krogh unit adjacent to the epidermis/dermis interface. As a result, the water in extra-cellular compartment began to freeze and the concentration difference between vascular and cellular compartment gets higher. The water in the cellular compartment transported from cell to the extra-cellular compartment and the radius of the extra-cellular compartment becomes bigger. The mass transport at deeper tissue depth occurred subsequently as the cooling conducted inward. As we can see from the figure, the radius of Krogh unit reaches the maximum firstly at skin depth of 61 mm with 3.3 s spurt duration. On the other hand, the probability of IIF is low (Fig. 11(b)). This indicted that cold injury occurred in dermis as cell dehydration with the spurt duration as long as 3.3 s. Fig. 12 plots predicted extra-vascular radius of Krogh unit (Fig. 12(a)) and IIF probability (Fig. 12(b)) as a function of R407c spray duration at different tissue depths: 61 mm, 83 mm, 105 mm, 127 mm, 149 mm and 171 mm, respectively. As can be seen from Fig. 12(a), the mass transport of the intracellular water to extra cellular space occurred after 92 ms spurt duration. Again, no IIF is predicted (Fig. 12(b)) and cell dehydration is predicted with 2.2 s spurt duration (Fig. 12(b)). Fig. 13 illustrates predicted extra-vascular radius of Krogh unit (Fig. 13(a)) and IIF probability (Fig. 13(b)) as a function of R404a spray duration at different tissue depths: 61 mm, 83 mm, 105 mm,

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

127 mm, 149 mm and 171 mm, respectively. The mass transport of the intracellular water to extra cellular space occurred after 80 ms spurt duration. For R404a, no IIF is predicted (Fig. 13(b)) and cell dehydration is predicted with 1.9 s spurt duration (Fig. 13(b)). The cooling effect of cryogen with lower boiling point is stronger that with higher boiling point. However, excessive spray duration may lead to the cooling of dermis and cause the phase change as the spurt duration prolonged. In summary, the phase change in dermis would initiate with 175 ms R134a spray, 92 ms R407c spray and 80 ms R404a spray based on the simulation result. Recognized as cell dehydration, the spray durations to cause cold injury are 3.3 s, 2.2 s and 1.9 s for R134a, R407c and R404a. Here the cold injury of cell dehydration is predicted prior to the IIF due to the low thermal conductivity of human skin and resulting in the low cooling rate within the dermal tissue. Cryogen spray cooling used in laser treatment of PWS is typical of tens of milliseconds. Therefore, the cold injury is not likely to occur during the spray even R404a is used. However, as an efficient surface cooling technology in laser dermatology, spurt duration should be carefully controlled to avoid the potential cold injury in further application of the cryogen spray cooling in laser clinic.

4. Conclusion In laser dermatologic surgery, cryogen spray cooling (CSC) is used to avoid unwanted damage such as scars from skin burning due to the melanin absorption of the laser beam. To describe the cooling process of cryogen spray cooling with various cooling agents, we presented a multi-scale model by coupling a micro-scale mass transport model to predict intracellular ice formation during cooling based on Krogh model with a macro-scale heat transfer model to predict the heat transfer in tissue skin during CSC. In the macro-scale model, a general dynamic relation of surface heat transfer coefficient from experimental data was used to quantify the surface cooling during cryogen spray cooling. The temperature dependent thermal properties of water and the biophysical properties of the tissue are also considered in the model. By using the above model, R134a, R407c and R404a which have different boiling point were compared. Important conclusions and suggestions to clinicians are as follows: 1) The surface cooling effect of spray cooling is well related with the boiling point of cryogen. The skin surface temperature would achieve the minimum value with 30 ms spurt duration and then increase to the boiling point temperature of the cryogens. The skin surface temperature will be stable on boiling point of cryogens after the formation of the cryogen pool on the skin surface other than further decreased as spray duration prolongs. However, the temperature of epidermis/dermis interface would keep decreasing. 2) Strong spatial inhomogeneity of the skin temperature is observed during cryogen spray cooling. Large temperature gradient in space (tissue) is caused by the low thermal conductivity of skin and the short spray duration. A lower surface and inner skin temperature will be achieved by using cryogens with lower boiling poring, e.g. R404a and R407c, which benefits the thermal damage protection. Epidermis cooling will get less efficient as the spray duration prolongs and the dermal temperature will decrease as the spurt duration prolonged. As the spray duration further increase, the temperature of dermis will fall below the freezing point of the water and potential cold injury could occur. Furthermore, the final coagulation of the PWS blood vessels will be weakened by the excessive cooling in dermis, which is not beneficial to the laser treatment of PWS.

11

3) Based on the numerical simulation for cold injury, the phase change in dermis would initiate with 175 ms R134a spray, 92 ms R407c spray and 80 ms R404a spray. By quantifying the microscale mass transfer and IIF, we found that the spray durations to cause cold injury (recognized as cell dehydration) are 3.3 s, 2.2 s and 1.9 s for R134a, R407c and R404a, respectively. 4) Cryogens with lower boiling point such as R407c and R404a are recommended to be used in clinically treatment of PWS because they can provide much stronger cooling protection for epidermis than traditionally used R134a and cold injury in dermal tissue will not take place with the spurt duration clinically used.

Acknowledgements The work is supported by Joint Research Fund for National Natural Science Foundation of China (51336006, 51228602), the International Science &Technology Cooperation Plan of Shaanxi Province (No. 2013KW30-05) and the Fundamental Research Funds for the Central University.

References [1] R.R. Anderson, J.A. Parrish, Microvasculature can be selectively damaged using dye laser: a basic theory and experimental evidence in human skin, Lasers Surg. Med. 1 (1981) 263e276. [2] J.W. Tunnell, J.S. Nelson, J.H. Torres, B. Anvari, Epidermal protection with cryogen spray cooling during high fluence pulsed dye laser irradiation: an ex vivo study, Lasers Surg. Med. 27 (2000) 373e383. [3] G. Aguilar, B. Majaron, W. Verkruysse, Y. Zhou, J.S. Nelson, E.J. Lavernia, Theoretical and experimental analysis of droplet diameter, temperature, and evaporation rate evolution in cryogenic sprays, Int. J. Heat. Mass Transf. 44 (2001) 3201e3211. [4] G. Aguilar, W. Verkruysse, B. Majaron, L.O. Svaasand, E.J. Lavemia, J.S. Nelson, Measurement of heat flux and heat transfer coefficient during continuous cryogen spray cooling for laser dermatologic surgery, IEEE J. Sel. Top. Quant. Electron. 7 (2001) 1013e1021. [5] B.M. Pikkula, J.H. Torres, J.W. Tunnell, B. Anvari, Cryogen spray cooling: effects of droplet size and spray density on heat removal, Lasers Surg. Med. 28 (2001) 103e112. [6] J.S. Nelson, T.E. Milner, B. Anvari, B.S. Tanenbaum, S. Kimel, L.O. Svaasand, S.L. Jacques, Dynamic epidermal cooling during pulsed laser treatment of portwine stain: a new methodology with preliminary clinical evaluation, Arch. Dermatol. 131 (1995) 695e700. [7] J.S. Nelson, B. Majaron, K.M. Kelly, Active skin cooling in conjunction with laser dermatologic surgery: methodology and clinical results, Semin. Cutan. Med. Surg. 19 (2000) 253e266. [8] J.S. Nelson, T.E. Milner, B. Anvari, B.S. Tanenbaum, L.,O. Svaasand, S. Kimel, Dynamic epidermal cooling in conjunction with laser-induced photothermolysis of port wine stain blood vessels, Lasers Surg. Med. 19 (1996) 224e 229. [9] T. Dai, M.A. Yaseen, P. Diagaradjane, Comparative study of cryogen spray cooling with R-134a and R-404a: implications for laser treatment of dark human skin, J. Biomed. Opt. 11 (2006) 04116 1e04116 11. [10] Z.F. Zhou, W.T. Wu, B. Chen, G.X. Wang, L.J. Guo, An experimental study on the spray and thermal characteristics of R134a two-phase flashing spray, Int. J. Heat. Mass Transf. 55 (2012) 4460e4468. [11] H.H. Zenzie, G.B. Altshuler, M.Z. Smirnov, et al., Evaluation of cooling methods for laser dermatology, Lasers Surg. Med. 26 (2000) 130e144. [12] W. Verkruysse, B. Majaron, B.S. Tanenbaum, J.S. Nelson, Optimal cryogen spray cooling parameters for pulsed laser treatment of port wine stains, Lasers Surg. Med. 27 (2000) 165e170. [13] D. Li, Y.L. He, G.X. Wang, J. Xiao, Y.W. Liu, Numerical analysis of cold injury of skin in cryogen spray cooling for dermatologic laser surgery, in: 2007 ASME International Mechanical Engineering Congress, 2007. Seattle, Washington. [14] F. Sun, G. Aguilar, K.M. Kelly, G.X. Wang, Thermal analysis for Cryosurgery of Nodular basal cell Carcinoma, 2006 ASME International Mechanical Engineering Congress and Exposition, Chicago, Illinois, USA. [15] B. Rubinsky, D.E. Pegg, A mathematical model for the freezing process in biological tissue, Proc. Phys. Soc. Lond. Sect. B 234 (1988) 343e358. [16] G. Aguilar, S.H. Diaz, E.J. Lavernia, J.S. Nelson, Cryogen spray cooling efficiency: improvement of port wine stain laser therapy through multiple-intermittent cryogen spurts and laser pulses, Lasers Surg. Med. 31 (2002) 27e35. [17] T.J. Pfefer, D.J. Smithies, T.E. Milner, M.J.C. van Gemert, J.S. Nelson, A.J. Welch, Bioheat transfer analysis of cryogen spray cooling during laser treatment of port wine stains, Lasers Surg. Med. 26 (2000) 145e157.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034

12

D. Li et al. / Applied Thermal Engineering xxx (2014) 1e12

[18] B. Majaron, W. Verkruysse, K.M. Kelly, J.S. Nelson, Er:YAG laser skin resurfacing using repetitive long-pulse exposure and cryogen spray cooling: II theoretical analysis, Lasers Surg. Med. 28 (2001) 131e137. [19] G. Aguilar, G.X. Wang, J.S. Nelson, Effect of spurt duration on the heat transfer dynamics during cryogen spray cooling, Phys. Med. Biol. 48 (2003) 2169e 2181. [20] W. Franco, J. Liu, G.X. Wang, J.S. Nelson, G. Aguilar, Radial and temporal variations in surface heat transfer during cryogen spray cooling, Phys. Med. Biol. 50 (2005) 387e397.

[21] D. Li, Y.L. He, Y.W. Liu, G.X. Wang, Numerical analysis of cryogen spray cooling of skin in dermatologic laser surgery using realistic boundary conditions, in: 22nd International Congress Refrigeration Conference, 2007 (BeiJing, China). [22] J.C. Bischof, B. Rubinsky, Microscale heat and mass transfer of vascular and intracellular freezing in the liver, ASME J. Heat. Transf. 115 (1993) 1029e1035. [23] M. Toner, E.G. Cravalho, M. Karel, Thermodynamics and kinetics of intracellular ice formation during freezing of biological cells, J. Appl. Phys. 67 (1990) 1582e1593. [24] M.N. Ozisik, Heat Conduction, second ed., John Wiley & Sons, New York, 1993.

Please cite this article in press as: D. Li, et al., Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Applied Thermal Engineering (2014), http://dx.doi.org/10.1016/j.applthermaleng.2014.03.034