The performance of droplet evaporation model in predicting droplet dynamics and thermal characteristics for R134a single isolated droplet and two-phase flashing spray

Aerospace Science and Technology 93 (2019) 105363

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

The performance of droplet evaporation model in predicting droplet dynamics and thermal characteristics for R134a single isolated droplet and two-phase flashing spray Zhi-Fu Zhou, Guan-Yu Lu, Dong-Qing Zhu, Lu Zhang, Jia-Feng Wang, Bin Chen ∗ State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an, 710049, China

a r t i c l e

i n f o

Article history: Received 27 February 2019 Received in revised form 11 June 2019 Accepted 23 August 2019 Available online 29 August 2019 Keywords: Evaporation Single isolated droplet Flashing spray R134a OpenFOAM

a b s t r a c t Flashing spray with low saturation point and high volatility mediums is of great importance in aerospace field. It involves complex droplet dynamics and heat and mass transfer processes in a turbulent, two phase flow. This paper comparatively evaluates the predictive performance of a selected number of droplet evaporation models that focus on convective and blowing effects. The studies span from a single, isolated R134a droplet that evaporates in a convective environment, to a fully turbulent, flashing spray formed through an accidental release of high pressure R134a liquid. An in-house developed code for single isolated droplet evaporation and a modified sprayFoam solver in OpenFOAM for flashing spray are used to calculate droplet and spray behaviors. The results show droplet evaporation model greatly affects the evolutions of droplet diameter, velocity and temperature for single isolated R134a droplet, that the C-R-S model predicts the lowest droplet diameter and velocity, and highest droplet temperature; the HN-R model predicts the largest droplet velocity and lowest temperature; the A-S and N-G-R-M models predict almost identical results. In contrast to the great impact on droplet evolution for single isolated droplet modeling, droplet evaporation model has little influence on spray and thermal characteristics for R134a two phase flashing spray simulation. However, the A-S model predicts quite different radial profile of droplet temperature at spray periphery compared with other models, which is much lower than the experimental value. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction Flashing spray occurs when a high-pressure liquid or liquefied gas is injected into the low-pressure environment to make the liquid superheated, characterized by explosive atomization of superheated liquid to generate fine droplets [1]. Meanwhile, strong evaporation leading to rapid changes in droplet temperature and diameter is an important feature, which distinguishes this phenomenon from other traditional pressure or hydraulic spray. Flashing spray is a common phenomenon and has important applications in many industrial fields. In the aerospace field, flashing spray is a key process in the thrust chamber of microthrusters or liquid rocket engines [2,3]. The refrigerant R134a is usually used as propellant for vaporizing liquid microthrusters, a micro-propulsion system applied to small satellites [4]. R134a is also one type of operating gases in air separators or vortex tubes that are used to provide an appropriate air stream to helicopter engines [5]. In the thermal

*

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

https://doi.org/10.1016/j.ast.2019.105363 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

management field for ultra-high power devices such as electronic and LED chips and laser weapons, flashing spray with R134a has demonstrated to be an effective way to solve the extreme heat dissipation problems [6,7]. In safety area, the behavior and characteristics of flashing two-phase spray significantly influence the hazard zone due to the accidental release of high-pressure flammable and poisonous liquid, and R134a is commonly used to generate the flashing spray to simulate the accidental release [8,9]. In cutaneous laser surgeries, especially for the treatment of port wine stain (PWS), cryogen spray cooling with R134a has been widely used for its excellence in protecting epidermis and improving the outcomes of laser treatment [10]. The importance of numerous applications has motivated researches on the evaporation and spray characteristics of R134a flashing spray, most of which was carried out through experimental means. Aguilar et al. [11,12] studied spray morphology and the variations of droplet diameter, velocity and temperature along the axial distance for R134a flashing spray with different nozzles by non-intrusive and intrusive methods, e.g. phase Doppler particle analyzer (PDPA) and fine thermocouple. Yildiz et al. [13–16] also conducted similar study to investigate droplet size, velocity and

2

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

Nomenclature BM BT c CD D g hs L Nu m P Pr R Re Sc Sh t T U

Spalding mass transfer number Spalding heat transfer number specific heat drag force coefficient droplet diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m gravitational acceleration . . . . . . . . . . . . . . . . . . . . . . . . . m/s2 specific enthalpy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J/kg latent heat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J/kg Nusselt number mass pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pa Prandtl number radius or radial distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Reynolds number Schmidt number Sherwood number time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . s temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s

temperature distributions of R134a flashing spray using different diameter nozzles. Vu et al. [17] further investigated the effect of nozzle length on spray behavior of R134a flashing spray, but no obvious difference of external spray characteristics was observed. Zhou et al. [1] conducted a comprehensive experimental study on the spray and thermal characteristics of R134a flashing spray and found a warm core in spray central region near the nozzle exit. Further, Zhou et al. [18] did a comparative study on the spray behavior and heat transfer dynamics using different volatile cryogens (R134a, R407C, and R404A), and found that droplets in spray with large velocity, small diameter, and low temperature led to a strong cooling effect on the surface. Experimental studies have provided basic data for R134a flashing spray, however, current optical methods, e.g. PDPA, phase Doppler Anemometry and particle image velocimetry, fail to provide reliable data of spray information in the near nozzle field regions due to the non-equilibrium nature and high spray density. Moreover, it is extremely difficult to track single droplet in a dense flashing spray using current measurement methods. Numerical simulation offers an alternative to explore the behavior and mechanism of R134a flashing spray. Wang et al. [19] developed a three-dimensional vortex method to investigate R134a flashing spray near nozzle field, while flash boiling phenomenon of superheated liquid was not taken into consideration in this study. Zhou et al. [20] conducted a numerical study on R404A flashing spray using OpenFOAM and considered the main phenomena of droplet transport, breakup, flashing boiling, evaporation, and heat transfer as well as the interaction between liquid and gas phases involved the complex flashing spray. Recently, Zhou et al. [21] proposed a coupled droplet evaporation model to predict the droplet temperature variation of R407C flashing spray, rather than the complicated CFD simulation, considering the coupling of heat and mass transfer between a droplet surface and its surrounding region of influence. Their study demonstrated that the coupled evaporation model could greatly improve the prediction of droplet temperature compared with the one-way droplet evaporation model. In flashing spray, droplet evaporation plays a decisive role in determining momentum, heat and mass transfer between the liquid phase and its ambient gas phase, especially for low saturation and high volatile mediums, e.g. R134a with saturation temperature of ∼ −26.2 ◦ C at atmospheric pressure. The explosive atomization and intense evaporation will make R134a droplet tempera-

S Y

source term mass fraction

Greek symbols

α Γv λ

μ ρ

thermal diffusivity or overall heat transfer coefficient binary diffusion coefficient thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . W/(m·K) dynamic viscosity density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3

Subscripts atm b eff g l s v

∞

ambient gas boiling effective ambient gas liquid surface vapor far field

ture decrease continuously due to that the convective heat transfer from ambient gas can not provide the large amount of latent heat but absorbing heat from its remaining liquid for evaporation, which distinguish it from other water and gasoline sprays [22]. Thus, droplet evaporation is crucial for the accurate prediction of spray and thermal characteristics for flashing spray. There are many evaporation models proposed continuously trying to predict the droplet evolution, especially for single droplet evaporation and heating. Since the first steady state model of single droplet evaporation in a stationary gas environment by Spalding [23], lots of efforts have been made to investigate the effects of convection, Stefan flow, blowing and turbulence on the process of heat and mass transfer between droplet and its ambient gases. Heat and mass transfer of an evaporating droplet is traditionally expressed by the correlations in forms of Nusselt (Nu) and Sherwood (Sh) numbers. Progress involved in various respects of droplet evaporation has been reviewed by Faeth [24], Law [25], Chiu [26], and Sazhin [27,28]. Miller et al. [29], Sazhin et al. [30] and Zhou et al. [31] conducted comparative analysis on the predictive performance of typical evaporation models for single droplet modeling with high saturation temperature, e.g. water and fuels. They found that all models predicted nearly same results and agreed well with experimental data in the cases of low evaporation rate, while great deviation among model predictions appeared at high evaporation rates. Although great progress of droplet evaporation model has been made in the past, rare research focuses on the model choice on the evaporation for low saturation temperature and high volatile cryogen droplet. Therefore, there is still a lack of consistency on how to select an appropriate evaporation model for cryogen droplet, and how the evaporation model influences the dynamics and heat and mass transfer for single isolated droplet and droplets in a flashing spray. For this reason, this study evaluates a wide range of typical droplet evaporation models focusing on the convective and blowing effects on heat and mass transfer for R134a droplet. The comparison and evaluation of the considered models is conducted following a two-stage procedure. Firstly, the performance of evaporation models in predicting temporal diameter, velocity and temperature evolutions for single isolated droplet is assessed through a in-house developed C++ code. Then, a modified sprayFoam solver in the open source tool OpenFOAM is used to evaluate the performance of evaporation models in predicting droplets diameter,

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

3

Table 1 Droplet evaporation models used in single droplet modeling and flashing spray simulation. Model Classical

Nu ( f T Nu0 )

Sh ( f M Sh0 )

1 ln(1+ B T ) (2 + 0.552Re 2 BT

1

A-S

1

1

2 Pr 3 ln(1+ B T ) (2 + 0.552Re ) BT F (B T ) 0.7 ln(1+ B T ) F ( B T ) = (1 + B T ) BT

H-N-R

(1 +

C-R-S

(1 + B T )−0.7 (2 + 0.454Re0.615 Pr0.98 )

N-G-R-M

(1 + B T )− 3 (2 + 0.552Re 2 Pr 3 )

Asano

1

Sc 3 )

0 < Re < 200 Elwardany et al. [33]

1

2 Sc 3 ln(1+ B M ) (2 + 0.552Re ) BM F (B M ) 0.7 ln(1+ B M ) F ( B M ) = (1 + B M ) BM

1 1 B T )−0.7 (2 + 0.57Re 2 Pr 3

2

Comments

1 ln(1+ B M ) (2 + 0.552Re 2 BM

1

Pr 3 )

)

(1 +

1 1 B M )−0.7 (2 + 0.87Re 2 Sc 3

0 < Re < 2000 Abramzon et al. [34]

)

(1 + B M )−0.7 (2 + 0.39Re0.54 Sc0.76 ) 1

1

2

30 < Re < 200 Chiang et al. [36] 1

1

(1 + B M )− 3 (2 + 0.552Re 2 Sc 3 )

1 2+0.552Re Pr 3 0.3+0.7(1+ B T )0.88

1 2+0.552Re Sc 3 0.3+0.7(1+ B M )0.88

1 2

1 2

velocity and temperature distributions along axial and radial directions in a R134a flashing spray, through Reynolds averaged Navier Stokes simulation with Eulerian Lagrangian framework for continuous gas and discrete droplets phases respectively.

0 < Re < 100 Haywood et al. [35]

Zhou et al. [31] 0 < Re < 1000 Asano [37]

where c pg and L are the specific heat of the mixture of gases near the droplet surface and the latent heat of evaporation. The momentum equation of the droplet, taking into account the drag force due to the relative motion between the droplet and its surrounding gas while gravity is negligible, is expressed as follows,

2. Single isolated droplet evaporation modeling

dU Several assumptions are first made to simplify the analysis for single droplet evaporation: 1) Droplet is spherical with uniform temperature distribution throughout droplet; 2) Ambient gases are not soluble in droplet; 3) Radiation heat transfer is neglected; 4) Droplet interface is always at local equilibrium state. According to Fick’s laws of mass diffusion and Stefan’s flow effect, the droplet vaporization rate dD /dt can be written as [32]:

dD Γv ρ g = −2 Sh ln(1 + B M ) dt D ρl

(1)

where Sh, t, D and Γ v are the Sherwood number, time, droplet diameter, and binary diffusion coefficient of liquid vapor in the ambient gas, respectively. ρ g and ρl are the average density of the mixture of gas and vapor near droplet surface and droplet density, respectively. B M is the Spalding mass transfer number, which is calculated as:

BM =

Y v ,s − Y v ,∞

(2)

1 − Y v ,s

Y v ,s is the vapor mass fraction on the droplet surface which can be obtained from the Clausius-Clapeyron equation [29], and Y v ,∞ is the vapor mass fraction at the far field away from the droplet surface, which is equal to zero. Since the droplet temperature is uniform within droplets, heat transfer only occurs at the droplet interface through convective heat transfer between the droplet and the surrounding ambient gases. Therefore, the energy equation for droplet evaporation is given as follows:

dT dt

  Nuλ g 1 1 dD = −( T − T g ) 6 2 +3 D

ρl c pl

B T D dt

c pg ( T g − T ) L

=−

3 ρg CD 4

ρl D

(U − U g )2

(5)

where U , U g and C D are the droplet velocity, ambient gas velocity (U g = 0) and drag force coefficient. And the equation of C D is shown below.

 CD =

24 (1 + 0.15Re0.687 ) Re

Re ≤ 1000

0.44

Re > 1000

(6)

The correlations of Nu and Sh numbers determine the heat and mass transfer on the surface of an evaporating droplet according to Equations (1) and (3), which refers to the evaporation model or evaporation gas phase model. A representative selection of such evaporation models examined in this study are shown in Table 1, which are deduced from the analytic and numerical analysis. All droplet evaporation gas phase models mainly consider the effects of convection and blowing on heat and mass transfer during evaporation process. The convective effect is due to the relative motion between droplet without evaporation or solid particle and the ambient gas or liquid, and is commonly expressed by Nu0 and Sh0 depended only on Reynolds and Prandtl numbers (Re and Pr). The blowing effect is caused by the vapor moving normal to droplet surface due to the intense evaporation, i.e. mass transfer, which changes the boundary layer on droplet surface and thus influences heat and mass transfer. As a result, the correction factors f T and f M as a function of B T and B M are introduced into the equations of Nu0 and Sh0 , accounting for the blowing effect on heat and mass transfer. Therefore, all droplet evaporation models can be expressed in the form of f T Nu 0 and f M Sh0 for vaporizing droplet. The details of these models can be found in our previous work [31].

(3) 3. Numerical modeling of R134a flashing spray

where Nu, T and T g are the Nusselt number, droplet and ambient gas temperatures; λ g and c pl are the thermal conductivity of the mixture gases near the droplet surface and liquid specific heat; and B T is the Spalding heat transfer number.

BT =

dt

(4)

3.1. Gas phase equations The continuity conservation is given by

∂ρ + ∇ · (ρ U ) = S evap ∂t

(7)

4

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

where ρ is density, U is the velocity vector, ∂∂t is the partial derivative with respect to time, and S evap is the source term that is from droplets evaporation. The momentum conservation equation uses the Reynolds Average Navier Stokes assumption, and S mom is the momentum source term. The turbulence model used in the numerical study is the standard k-ε model.

∂ρU + ∇ · (ρ U U ) − ∇ · μeff ∇ U ∂t    2  − ∇ · μeff (∇ U )T − tr (∇ U )T I = ρ g − ∇ p + S mom

Table 2 Initial parameters for R134a single isolated droplet modeling.

where μeff is the effective dynamic viscosity, ∇· is the divergence operator, ∇ is the gradient operator, T stands for the transpose operation, tr denotes the trace operator, I is the identity matrix, p and g are pressure and gravitational acceleration. The species conservation equation is given by

∂ρYi + ∇ · (ρ U Y i ) − ∇ · μeff ∇ Y i = S evap,i ∂t

(9)

where S evap,i is the source term that is from spray evaporation. The energy conservation equation is given by

∂ ρ hs Dp + ∇ · (ρ U h s ) − ∇ · (ρ U h s ) − ∇ · αeff ∇ h s = + S heat ∂t Dt (10) where h s is specific enthalpy,

αeff is effective thermal diffusivity,

Dp Dt

is the total derivative of p, and S heat is the source term that is from heat transfer model. 3.2. Liquid phase models The major liquid models considered in the numerical simulations are droplet motion, droplet breakup, droplet evaporation/boiling, droplet heat transfer and injection models [20]. The model concerning droplets flashing evaporation is described as below. High-pressure R134a is injected into the low-pressure environment to make it superheated quickly. Therefore, apart from the equilibrium evaporation, the non-equilibrium boiling process should also be considered. When the saturated vapor pressure of the liquid is lower than the ambient pressure, the evaporation model should be used in the simulation; on the contrary, when it is higher than the ambient pressure, the non-equilibrium boiling model should be employed [38,39]. For the boiling condition, the overall evaporation rate derives ˙ flash ) and the evaporation from superheating evaporation rate (m ˙ t) rate due to external heat transfer (m

dm dt

˙ flash + m ˙t =m

(11)

According Adachi’s experimental correlations [40], the flash evapo˙ flash ) is given by ration rate (m

˙ flash = π D 2 α S m

(T − T b ) L

(12)

where α S is the overall heat transfer coefficient, T b is boiling temperature. The overall heat transfer coefficient is also given by Adachi correction [40]

⎧ 0.26 ⎪ ⎨ 760 × ( T − T b ) , 0 < T − T b < 5 2.33 α S = 27 × ( T − T b ) , 5 ≤ T − T b < 25 ⎪ ⎩ 13800 × ( T − T b )0.39 , 25 ≤ T − T b

(13)

Value

Ambient temperature, T ambi (K) Ambient pressure, P ambi (MPa) R134a droplet initial temperature, T 0 (K) R134a droplet initial velocity, V 0 (m/s) R134a droplet initial diameter (μm)

298.15 0.1 246.15 60 100

Table 3 Initial parameters for R134a two phase flashing spray simulation.

(8)

3

Parameter

Parameter

Value

Ambient temperature, T ambi (K) Ambient pressure, P ambi (MPa) R134a initial temperature, T 0 (K) R134a injection pressure, P inj (MPa) R134a injection angle (◦ ) Diameter of nozzle, dnozzle (mm) Rosin Rammler distribution parameter d (μm) Rosin Rammler distribution parameter n Duration of injection, t dur (ms) Parcels per second

298.15 0.1 298.15 0.7 90 0.81 12 1.7 50 20,000,000

˙ t ) is given The evaporation rate due to external heat transfer (m by

˙t =πD m

λg

Sh





˙ flash m

ln 1 + 1 + ˙t c pg 1 + m˙ flash m ˙t m



 BT

(14)

The evaporation rate under the normal equilibrium evaporation is given by equation (1). Energy equation of heat transfer between droplets and continuous gas is same with that for single isolated droplet, i.e. equation (3). The injection, breakup and motion models are Rosin-Rammler, Reitz Diwakar and O’Rourke models, details of which can be found in [20]. 4. Numerical methods The modeling of R134a single isolated droplet evaporation has been performed using in-house developed C++ code to calculate the mass, energy and momentum equations through a complete implicit iteration algorithm. The initial parameters for modeling are shown in Table 2. As for the numerical simulation of R134a two phase flashing spray, it has been performed using modified sprayFoam solver in OpenFOAM-2.4.0 open source code, based on the Eulerian Lagrangian framework for the continuous and discrete phases respectively. The gas phase is computed through solving the Reynoldsaveraged Navier Stokes equations. For the liquid phase, Lagrangian approach is adopted to trace the evolution of parcels, each of which represents a group of real droplets with the same properties, i.e. diameter, velocity and temperature. As the interaction between the two phases, a two-way coupling is introduced in the numerical study, where source terms calculated in the sub-models using Lagrangian approach, are introduced into gas phase equations in Eulerian approach. The initial parameters for the numerical simulation of R134a flashing spray are shown in Table 3, based on our previous experiment conditions [1]. The computation domain is a rectangular prism with the size of 200 × 200 × 500 mm3 with uniform grid, at the top center of which the single-hole nozzle is set. To verify the grid independence, three different size grids are tested and spray tip penetration is used as the basis of comparison. As shown in Fig. 1, no apparent distinction can be seen with grid number larger than 101 × 101 × 250, therefore, the grid used in this study is 101 × 101 × 250, about 2,550,250 cells totally.

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

5

Fig. 1. Spray tip penetration for various grid sizes.

5. Results and discussions 5.1. Single isolated R134a droplet evaporation modeling results In this section, we studied the effect of evaporation model on the evolutions of diameter, velocity and temperature for single isolated R134a droplet. Fig. 2 presents single isolated droplet modeling results predicted by six different evaporation models. As shown in Fig. 2a, although all diameter curves show similar variation trend, droplet diameter variation along distance is greatly influenced by the evaporation models except that the A-S and NG-R-M models predict nearly same results. At the initial stage of droplet evaporation within ∼ 80 mm, the slope of droplet diameter presents a gradual decrease with distance increasing, while it seems to have the opposite trend at the later evaporation stage with distance larger than 80 mm. The H-N-R and classical models predict the largest and second largest evaporation rates at the initial evaporation stage, while the Asano model predicts the smallest evaporation rate. The results predicted by the C-R-S, A-S and N-GR-Y models lie in the middle of them. As the droplet goes into the latter evaporation stage, the six evaporation models contribute to three groups of results. The C-R-S evaporation model predicts noticeably smallest droplet diameter (about 76.9 μm) at the distance of 200 mm; the Asano evaporation model predicts the second smallest droplet diameter (80.6 μm) at the distance of 200 mm; and droplet diameter obtained by other four evaporation models gradually becomes close along the distance and reach about 81.7 μm at the distance of 200 mm. In summary, the effect of the choice of evaporation model on the diameter variation of single isolated R134a droplet seems to be relatively great at the initial evaporation stage and become relatively weak at the last 100 mm distance. The variation of droplet velocity along the distance is shown in Fig. 2b. All six evaporation models also predict a similar trend of droplet velocity variation that both droplet velocity and its decreasing rate gradually decrease with distance increasing due to the drag force acting on droplet. However, discrepancy among model predictions still can be observed. The C-R-S evaporation model results in the lowest droplet velocity along the distance; the H-N-R evaporation models results in the highest droplet velocity along the distance, and the other four evaporation models predict droplet velocity between the results of C-R-S and H-N-R evaporation models. It is also noticed that the results predicted by the A-S and N-G-R-M evaporation models appear almost identical curves once again. Combining Fig. 2a and Fig. 2b, larger evaporation rate leads to smaller droplet velocity, probably attributed to the greater mass loss during droplet evaporation. As for the variation of droplet temperature along the distance, which is shown in Fig. 2c, all curves predicted by the six evap-

Fig. 2. Evolution of (a) droplet diameter, (b) droplet velocity, and (c) droplet temperature as the function of distance for single isolated R134a droplet predicted by six evaporation models.

oration models also show a similar pattern that both the value of droplet temperature and the slope of the curves gradually decrease along the distance. However, the choice of evaporation model has a more significant influence on the droplet temperature variation compared with its effect on droplet diameter and velocity. The predictive results by six evaporation models begin to differ with each other largely at the very beginning of droplet evaporation. The C-R-S evaporation model results in almost highest droplet temperature along the distance and H-N-R evaporation model results in the lowest droplet temperature along the distance. This can be explained by the result shown in Fig. 2b, which means lower droplet velocity generates higher droplet temperature and vice versa, due to the lower evaporation rate. Once again, the A-S and H-G-R-M

6

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

Fig. 4. Spray tip penetration as the function of time predicted by six evaporation models.

that of N-G-R-M evaporation model. As for other four evaporation models, there is no noticeable difference among them.

Fig. 3. Comparison of the spray profiles: (a) liquid phase of R134a flashing spray from simulation; (b) gas phase of R134a flashing spray from simulation; (c) R134a flashing spray image from experiment.

evaporation models predict almost identical temperature curves with the final steady value of ∼ −58.2 ◦ C. The Asano evaporation model results in the same value of droplet steady temperature with that predicted by A-S and H-G-R-M evaporation models, although it predicts higher temperature at the beginning. 5.2. R134a two-phase flashing spray simulation results In this section, we studied the effect of evaporation model on the simulation results of spray and thermal characteristics of R134a flashing spray using OpenFOAM. All the results shown in the following figures correspond to the steady state of R134a flashing spray at 30 ms injection duration except the time-dependent results of spray penetration distance. 5.2.1. Spray morphology All the droplet evaporation models give almost same numerical predictions of spray morphology, thus it is unnecessary to show all the predictive spray morphology. Fig. 3 shows the morphologies of R134a flashing spray predicted by the classical model reaching the steady state in a 30 ms spray, and the experimental result of spray image recorded by a high speed camera is also presented for comparison [1]. The numerical results (both the liquid phase and gas phase) show quick expansion once the liquid is released from the nozzle exit. Afterwards, the liquid phase keeps relative stable, while gas phase undergoes a slow but continuous expansion towards the lateral direction with axial distance increasing. Moreover, the spray central region has much higher concentrations of droplet and vapor near nozzle field compared with the periphery regions, while vapor concentration becomes more uniform in the downstream area. This kind of spray morphology agrees reasonably with the experimental result. Thus, it provides the validation of the two-phase flashing spray numerical simulations. Fig. 4 presents the effect of different evaporation models on the predictive spray tip penetration. As shown in Fig. 4, the evaporation model has very limited influence on spray penetration in R134a flashing spray simulation. The N-G-R-M evaporation model results in relative shorter spray penetration. Moreover, the spray penetration of A-S evaporation model gradually becomes close to

5.2.2. Axial variations of droplet diameter, velocity and temperature Fig. 5 presents the evolutions of droplet dynamics and thermal characteristics, i.e. droplet Sauter diameter (D 32 ), axial velocity and temperature along the spray centerline in R134a flashing spray simulation using six different evaporation models given in Table 1. As shown in Fig. 5a, all six evaporation models result in nearly same variations of D 32 , which decrease fast at the near nozzle field and then have small change in the downstream area (> 50 mm). The predictive results of D 32 agree reasonably with the experimental data obtained by PDPA in our previous study [1], which further validates the flashing spray simulation. The effect of evaporation model on D 32 in R134a flashing spray simulation can be negligible, although it has noticeable influence on droplet diameter variation of single isolated droplet as shown in Fig. 2a. There might be at least two sound reasons contributing to this phenomenon. Firstly, the breakup process that is taken into consideration in the numerical simulation of R134a flashing spray is neglected in single droplet modeling, but breakup phenomenon is vitally important in real flashing spray [41]. Secondly, D 32 is a statistical quantity with the statistical significance of multiple droplets in spray, thus the accumulation of individual differences may lead to the subtle difference as a whole. As shown in Fig. 5b, the variation curves of droplet axial velocity along the centerline are also quite similar. However, there are still some slight differences among the results predicted by the six evaporation models. During the distance of first 50 mm away from the nozzle exit, the N-G-R-M evaporation model results in the highest droplet velocity along the axial distance, while other five evaporation models produce almost same results. However, during the distance from 100 mm to 200 mm, the A-S evaporation model results in the highest droplet velocity along the axial distance, while other five evaporation models lead to almost same results. Once again, the effect of the choice of evaporation model on the variation of droplet velocity along the axial distance in numerical simulation of R134a flashing spray is apparently weaker than that in single isolated droplet modeling. Apart from the two reasons mentioned above, this could also be explained by the fact that the R134a vapor velocity is not calculated in the single droplet modeling while there is a two-way coupling between the vapor and liquid phases in spray numerical simulation. Fig. 5c shows the variation of droplet temperature along the axial distance in numerical simulation of R134a flashing spray using six evaporation models. The experimental data of droplets average temperature measured through a invasive thermocouple is also plotted for comparison [1]. There are some differences among the

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

7

Fig. 6. Variation of droplet Sauter mean diameter as the function of radial distance at four axial distances in R134a flashing spray simulation using six evaporation models.

the invasive measurement method of droplet temperature which has great impact on flashing spray, or probably attributed to the uncertainties of various models involved in the full numerical simulation. Nevertheless, the minimum value of droplet temperature predicted by N-G-R-M evaporation model is closest to the experimental data.

Fig. 5. Variations of (a) droplet Sauter diameter, (b) droplet axial velocity, and (c) droplet temperature as the function of axial distance in R134a flashing spray simulation using six evaporation models.

temperature results predicted by different evaporation models but the effect of evaporation model on the variation of droplet temperature in spray numerical simulation is still much weaker than that in single isolated droplet modeling. Although the N-G-R-M and A-S evaporation models generate almost identical droplet temperature curves in single isolated droplet modeling, the N-G-R-M evaporation model results in the highest droplet temperature and the A-S evaporation model results in the lowest droplet temperature along the axial distance in the spray numerical simulation. Moreover, the other four evaporation models seem to generate very similar droplet temperature curves. As for the comparison with the experimental data, the droplet minimum temperatures predicted by the six models are always a bit lower than the experimental data. This slight discrepancy is probably attributed to

5.2.3. Radial variations of droplet diameter, velocity and temperature Fig. 6 presents the variation of D 32 along the radial distance at four axial distances in R134a flashing spray simulation using six evaporation models. All the curves have the same variation trends at the three upstream sections, i.e. z = 50, 50 and 130 mm, that D 32 first decreases slightly and then increases with radial distance increasing. At the last section of 170 mm, D 32 increases monotonously from spray center. It is noticeable that A-S evaporation model predicts the largest D 32 in the region with radial distance larger than 5 mm, while N-G-R-M evaporation model results in the smallest D 32 . There is no obvious difference in D 32 predicted by other four models. In contrast to the significant impact on droplet diameter for single isolated droplet evaporation, the choice of droplet evaporation model has much less effect on the distributions of D 32 along both axial and radial distances. Fig. 7 depicts the variation of droplet axial velocity along the radial distance at four axial distances. The highest velocity always locates at the spray center for all spray sections, and droplet velocity decreases monotonously towards spray edge. The decreasing rate of velocity becomes small as the axial distance increases. As for the choice of evaporation models, it seems that their impact on the radial distribution of droplet velocity depends on spray distance. The radial distributions of droplet axial velocity at z = 50 and 130 mm are almost identical for six evaporation models, while the radial distributions differ a little among them at z = 90 and 170 mm. A-S evaporation model results in the lowest droplet velocity along the radial distance at z = 90 mm, while it results in the highest value along the radial distance at z = 170 mm. The choice of droplet evaporation model also has much less effect on droplet velocity in flashing spray simulation than in single isolated droplet evaporation modeling. Fig. 8 illustrates droplet temperature variation along the radial distance in R134a flashing spray simulation using six evaporation models. All evaporation models produce similar radial temperature profile at four axial distance that there is a relatively hightemperature region near the spray centerline and a relative low-

8

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

6. Conclusion In this study, the performance of droplet evaporation models in predicting droplet dynamics and heat and mass transfer for R134a single isolated droplet and two-phase flashing spray is examined systematically, based on an in-house developed code and a modified sprayFoam solver in OpenFOAM respectively. Six typical evaporation models are selected, focusing on the convective and blowing effect due to high evaporation rate. Numerical simulations of R134a flashing spray are validated by previous experimental results in terms of spray morphology and droplet diameter. The main conclusions can be summarized as follows:

Fig. 7. Variation of droplet axial velocity as the function of radial distance at four axial distances in R134a flashing spray simulation using six evaporation models.

1) The choice of evaporation model greatly affects the evolutions of droplet diameter, velocity and temperature in the modeling of R134a single isolated droplet, while the A-S and N-G-R-M evaporation models predict almost identical results among the selected typical evaporation models. The C-R-S model predicts the lowest droplet diameter and velocity, and highest droplet temperature. The H-N-R model predicts the largest droplet velocity and lowest temperature. 2) The predictive results of spray morphology and droplet Sauter mean diameter along spray centerline agree reasonably with the previous experimental results. In contrast to the great impact on droplet evolution for single isolated droplet, droplet evaporation model has much less effect on the numerical simulation results of spray penetration distance, the distributions of droplet diameter, velocity and temperature along both axial and radial directions for R134a two-phase flashing spray. 3) The A-S evaporation model predicts quite different radial profile of droplet temperature at spray periphery in R134a flashing spray simulation compared with other models, which is much lower than the experimental data. Thus, the A-S evaporation model is not suitable in the prediction of droplet temperature distribution, and other five evaporation models can be recommended to predict droplet dynamics and thermal characteristics of R134a flashing spray. Declaration of competing interest

Fig. 8. Variation of droplet temperature as the function of radial distance at four axial distances in R134a flashing spray simulation using six evaporation models.

No competing interests. Acknowledgements

temperature region near the spray edge. The temperature profile along the radial distance becomes uniform as the axial distance increases. Moreover, the effect of the choice of evaporation model on droplet temperature variation gradually decreases as the axial distance increases. At the region near the nozzle exit (z = 50 mm), the temperature profile along the radial distance differs a lot among the six evaporation models. The A-S evaporation model results in the lowest droplet temperature along the radial distance; the H-N-R evaporation model results in the second lowest droplet temperature along the radial distance; the N-G-R-M evaporation model results in the highest droplet temperature along the radial distance; and other evaporation models result in similar values. This pattern of difference in predictive droplet temperature by six evaporation models still exits as the axial distance increases, although the magnitude of difference becomes small. It must point out that the A-S evaporation model leads to much lower droplet temperature than those by other models, which predicts the lowest temperature below −70 ◦ C. It is also much lower than the experimental data of ∼ −60 ◦ C [1]. This large gap between predictive and experimental results indicates that the A-S evaporation model might not be suitable in predicting droplet radial temperature distribution for R134a flashing spray.

This study was supported by the Key Program for International Science & Technology Cooperation Projects from the Ministry of Science and Technology of China (Grant No. 2017YFE0134200) and the International Science & Technology Cooperation Plan of Shaanxi Province (Grant No. 2019KW-021). The authors also thank the fund from the State Key Laboratory of Engines of Tianjin University (K2019-04), China. References [1] Z. Zhifu, W. Weitao, C. Bin, W. Guoxiang, G. Liejin, An experimental study on the spray and thermal characteristics of R134a two-phase flashing spray, Int. J. Heat Mass Transf. 55 (15–16) (2012) 4460–4468. [2] J. Kang, J. Heo, H.G. Sung, Y. Yoon, Dynamic characteristics of a cryogenic nitrogen swirl injector under supercritical conditions, Aerosp. Sci. Technol. 67 (2017) 398–411. [3] A. Kebriaee, G. Olyaei, Semi-analytical prediction of macroscopic characteristics of open-end pressure-swirl injector, Aerosp. Sci. Technol. 82–83 (2018) 32–37. [4] M.G. De Giorgi, D. Fontanarosa, A novel quasi-one-dimensional model for performance estimation of a vaporizing liquid microthruster, Aerosp. Sci. Technol. 84 (2019) 1020–1034. [5] S.E. Rafiee, M.M. Sadeghiazad, Experimental and 3D CFD analysis on optimization of geometrical parameters of parallel vortex tube cyclone separator, Aerosp. Sci. Technol. 63 (2017) 110–122.

Z.-F. Zhou et al. / Aerospace Science and Technology 93 (2019) 105363

[6] Z.-F. Zhou, T.-Y. Xu, B. Chen, Algorithms for the estimation of transient surface heat flux during ultra-fast surface cooling, Int. J. Heat Mass Transf. 100 (2016) 1–10. [7] Z.-F. Zhou, Y.-K. Lin, H.-L. Tang, Y. Fang, B. Chen, Y.-C. Wang, Heat transfer enhancement due to surface modification in the close-loop R410A flash evaporation spray cooling, Int. J. Heat Mass Transf. 139 (2019) 1047–1055. [8] G. Polanco, A.E. Holdo, G. Munday, General review of flashing jet studies, J. Hazard. Mater. 173 (1–3) (2010) 2–18. [9] R.K. Calay, A.E. Holdo, Modelling the dispersion of flashing jets using CFD, J. Hazard. Mater. 154 (1–3) (2008) 1198–1209. [10] Z.-f. Zhou, B. Chen, R. Wang, F.-l. Bai, G.-x. Wang, Coupling effect of hypobaric pressure and spray distance on heat transfer dynamics of R134a pulsed flashing spray cooling, Exp. Therm. Fluid Sci. 70 (2016) 96–104. [11] 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 (17) (2001) 3201–3211. [12] G. Aguilar, B. Majaron, E. Karapetian, E.J. Lavernia, J.S. Nelson, Experimental study of cryogen spray properties for application in dermatologic laser surgery, IEEE Trans. Biomed. Eng. 50 (7) (2003) 863–869. [13] D. Yildiz, P. Rambaud, J.v. Beeck, Break-up, droplet size and velocity characterizations of a two-phase flashing R134a jet, in: 5th International Conference on Multiphase Flow, Yokohama, Japan, 2004. [14] D. Yildiz, P. Rambaud, J. Van Beeck, J. Buchlin, Evolution of the spray characteristics in superheated liquid jet atomization in function of initial flow conditions, in: 10th International Conference on Liquid Atomization and Spray Systems, Kyoto, Japan, 2006. [15] D. Yildiz, J.P.A.J. van Beeck, M.L. Riethmuller, Feasibility exploration of laserbased techniques for characterization of a flashing jet, Part. Part. Syst. Charact. 21 (5) (2004) 390–402. [16] D. Yildiz, P. Rambaud, J. Van Beeck, J. Buchlin, Thermal characterization of a R134A two-phase flashing jet, in: 9th International Conference on Liquid Atomization and Spray Systems ICLASS, 2003. [17] H. Vu, O. García-Valladares, G. Aguilar, Vapor/liquid phase interaction in flare flashing sprays used in dermatologic cooling, Int. J. Heat Mass Transf. 51 (23–24) (2008) 5721–5731. [18] Z.-F. Zhou, B. Chen, R. Wang, G.-X. Wang, Comparative investigation on the spray characteristics and heat transfer dynamics of pulsed spray cooling with volatile cryogens, Exp. Therm. Fluid Sci. 82 (2017) 189–197. [19] R. Wang, B. Chen, X.-S. Wang, Numerical simulation of cryogen spray cooling by a three-dimensional hybrid vortex method, Appl. Therm. Eng. 119 (2017) 319–330. [20] Z.-F. Zhou, G.-Y. Lu, B. Chen, Numerical study on the spray and thermal characteristics of R404A flashing spray using OpenFOAM, Int. J. Heat Mass Transf. 117 (2018) 1312–1321. [21] Z.-F. Zhou, M.-Y. Hu, H. Xin, B. Chen, G.-X. Wang, Experimental and theoretical studies on the droplet temperature behavior of R407C two-phase flashing spray, Int. J. Heat Mass Transf. 136 (2019) 664–673. [22] Y. Li, H. Guo, Z. Zhou, Z. Zhang, X. Ma, L. Chen, Spray morphology transformation of propane, n-hexane and iso-octane under flash-boiling conditions, Fuel 236 (2019) 677–685.

9

[23] D.B. Spalding, Experiments on the burning and extinction of liquid fuel spheres, Fuel 32 (1953) 169–185. [24] G.M. Faeth, Current status of droplet and liquid combustion, Prog. Energy Combust. Sci. 3 (4) (1977) 191–224. [25] C.K. Law, Recent advances in droplet vaporization and combustion, Prog. Energy Combust. Sci. 8 (3) (1982) 171–201. [26] H.H. Chiu, Advances and challenges in droplet and spray combustion. I. Toward a unified theory of droplet aerothermochemistry, Prog. Energy Combust. Sci. 26 (4–6) (2000) 381–416. [27] S.S. Sazhin, Advanced models of fuel droplet heating and evaporation, Prog. Energy Combust. Sci. 32 (2) (2006) 162–214. [28] S.S. Sazhin, Modelling of fuel droplet heating and evaporation: recent results and unsolved problems, Fuel 196 (Complete) (2011) 69–101. [29] R.S. Miller, K. Harstad, J. Bellan, Evaluation of equilibrium and non-equilibrium evaporation models for many-droplet gas-liquid flow simulations, Int. J. Multiph. Flow 24 (6) (1998) 1025–1055. [30] S.S. Sazhin, T. Kristyadi, W.A. Abdelghaffar, M.R. Heikal, Models for fuel droplet heating and evaporation: comparative analysis, Fuel 85 (12–13) (2006) 1613–1630. [31] Z. Zhifu, W. Guoxiang, C. Bin, G. Liejin, W. Yueshe, Evaluation of evaporation models for single moving droplet with a high evaporation rate, Powder Technol. 240 (2013) 95–102. [32] Z.-F. Zhou, W.-Y. Li, B. Chen, G.-X. Wang, A 3rd-order polynomial temperature profile model for the heating and evaporation of moving droplets, Appl. Therm. Eng. 110 (2017) 162–170. [33] A.E. Elwardany, S.S. Sazhin, A. Farooq, Modelling of heating and evaporation of gasoline fuel droplets: a comparative analysis of approximations, Fuel 111 (Supplement C) (2013) 643–647. [34] B. Abramzon, W.A. Sirignano, Droplet vaporization model for spray combustion calculations, Int. J. Heat Mass Transf. 32 (9) (1989) 1605–1618. [35] R.J. Haywood, R. Nafziger, M. Renksizbulut, A detailed examination of gas and liquid phase transient processes in convective droplet evaporation, J. Heat Transf. 111 (2) (1989) 495–502. [36] C.H. Chiang, M.S. Raju, W.A. Sirignano, Numerical analysis of convecting, vaporizing fuel droplet with variable properties, Int. J. Heat Mass Transf. 35 (5) (1992) 1307–1324. [37] K. Asano, Mass Transfer: From Fundamentals to Modern Industrial Applications, Wiley-VCH, Weinheim, 2006. [38] D. Kawano, H. Ishii, H. Suzuki, Y. Goto, M. Odaka, J. Senda, Numerical study on flash-boiling spray of multicomponent fuel, Heat Transf. Asian Res. 35 (5) (2006) 369–385. [39] B. Zuo, A.M. Gomes, C.J. Rutland, Modelling superheated fuel sprays and vaporization, Int. J. Eng. Res. 1 (4) (2000) 321–336. [40] M. Adachi, V.G. McDonell, D. Tanaka, J. Senda, H. Fujimoto, Characterization of Fuel Vapor Concentration Inside a Flash Boiling Spray, SAE International, 1997. [41] S.S. Sazhin, W.A. Abdelghaffar, E.M. Sazhina, M.R. Heikal, Models for droplet transient heating: effects on droplet evaporation, ignition, and break-up, Int. J. Therm. Sci. 44 (7) (2005) 610–622.