Initially sub-cooled liquid flashing flow CFD modeling, closure parameters for a non-symmetric interfacial area density formulation

A bubble to droplet transition non-symmetric interfacial area density formulation was developed for flash boiling nozzles CFD modeling. This formulation is used in a two fluid Euler-Euler modeling with a thermal phase change model. The regime transition formulation is based on the number of bubbles density and number of droplet density values. The first affects mainly the flow regime close to the nozzle throat and the second parameter affects mainly the flow close to the outlet. The physical factors affecting the calibration values of these parameters are analyzed in this paper with the aim of defining appropriate closure models. The analysis is done using experimental data available in literature. The comparison is done on the mass flow rate and the nozzle outlet velocity. It appears that the number of bubbles density is related to the inlet sub-cooling conditions and that the number of droplets density is affected by the outlet pressure and Reynolds number.


INTRODUCTION
Flash boiling nozzles are the motive elements of various energy conversion components such as twophase impulse turbines [1] or two-phase ejectors [2]. Both have been identified as key components for energy efficiency improvement in mechanical power and cold power production thermodynamic systems. The design and optimization of such components is very challenging due to the complexity of their internal two-phase flow. The behavior of flash nozzles is strongly related to the evaporation regime during the expansion observed by the flow. This affects the critical flow rate and the pressure, void fraction and velocity profiles of the flow [3]. Computational fluid dynamics tools have been extensively used for flash flows modeling and research on associated models are object of intense research as can be read in [2] or [4]. Besides, models for the interfacial mass transfer term are specially investigated either for nuclear safety applications [4] or for CO 2 two-phase ejectors [2].
Whatever the modeling approach is, a majority of models in literature need the use of one or various adaptation coefficients to adjust the model's predictions to reality. However when dealing with motive flash nozzles the type and/or number of adaptation coefficients may not be adequate to reproduce by simulation the observed nozzle performance. For example in CO 2 two-phase ejector modeling, the calibration of the motive nozzle flow rate can be in general attained with relatively low discrepancies (depending on the inlet enthalpy), but the entrainment capacity, mainly related to the motive nozzle outlet velocity, presents high discrepancies and is still subject to extensive research work. This can be read in a recent review paper of Ringstad et al. [5]. In critical flashing expansion conditions, the outlet conditions don't affect the upstream flow regime. So if adaptation of the evaporative model is needed close to the outlet it should avoid modifying the upstream evaporative situation. However this is not possible with models based on only one accommodation coefficient and hard to attain with others. A proof to that complexity is the development by some authors of specific calibration procedures based on genetic algorithms [6]. The purpose of this paper is to present a model allowing the adjustment of the evaporative term close to the nozzle throat and close to the nozzle outlet independently. The formulation describes the liquid-vapor interfacial area density as a function of 2 the void fraction and two adjustment parameters: the number of bubbles density and the number of droplets density. Using physically significant adjustment parameters aims to develop a generalist way of modeling such flows. The model is tested on literature water flashing nozzle flow rate and outlet velocity data obtained for various inlet sub-cooling degrees and outlet pressures. Beside, the paper analyses the adjustment parameters relations to the flow conditions in order to investigate closure relationships for the number of bubbles and the number of droplets densities definition. The objective is to propose a simpler-to-adjust and robust method for flashing nozzles performances prediction. The paper is organized as follows: first in the modeling section, the conservation equations used for the simulations are presented followed by the main literature interfacial mass transfer models presentation before introducing the proposed model. Then a simple sensitivity analysis of the models to their respective adaptation coefficients is conducted in order to discuss and compare the models. After that, in the simulation section, the experimental case study, the simulation configuration and the adjustment procedure are presented before discussing the results obtained by simulation. Finally the adjustment parameters analysis is proposed before the conclusion.

SET OF EQUATIONS
Conservation equations: The general formulation of conservation equations for a two phase flow requires a phase per phase description. In the commercial CFD code Ansys CFX 16.0 [7] these equations, describing the so called Euler-Euler formulation, are described by equations 1 to 4 for liquid. The same kind of equations is used for liquid and gas; index 1 is used for the liquid and 2 for the gas. is the interfacial heat transfer rate. The interfacial transfer terms and take into account the flux direction as follows:

4
The CFX solver solves also the volume continuity equation, or volume fraction transport equation, expression 5, usually introduced in multiphase problems. Thus, in the general case 6 (+1) equations are solved.

5
The index p refers to the phase (1 or 2) and the index p' to the remaining phase (i.e. 2 or 1). This two fluid approach has already been implemented for flash nozzles by Yixiang for example [8]. Reagrding the mass transfer term, a lot of formulations exist in literature. A summary of them is done in the next section.
Interfacial mass transfer models, review and discussion: This section presents the most popular mass transfer formulations for flashing flows modeling. These models are then discussed and a sensitivity test based on simple case is done. This analysis aims to highlight the effect of the adjustment parameters usually used in these models especially regarding the ability of the models to separately act on the early evaporation regime and the latter evaporation regime of the expanding flow. For flash boiling flows, the interfacial mass flux is in general formulated as a function of a pressure, temperature or enthalpy difference between the actual and the saturation states. The state difference, which is also illustrating the metastable situation of the liquid, is then associated to a relaxation term. Three types of formulations are discussed here. The list of models doesn't pretend to be exhaustive but represent, following the readings of the author, the most representative model studied today in literature 3 concerned by water flash nozzles and two-phase ejectors.
Pressure difference based models: A reference model for the pressure difference driven formulations is the Singhal model which is based on Rayleigh-Plesset bubble growth analysis. It was also used for depressurized flows [9]. However for the last the Hertz-Knudsen (HK) formulation is more frequently used. That is the case in the work of Karathanassis [10] and Le [11] on water nozzles or Yazdani [12] for CO 2 two phase ejectors. The general expression of the mass flux is based on the kinetic theory of gases and it is named Hertz-Knudsen relation: 6 The relaxation term described here by the first division can take a multitude of forms [10] [11] [12] but the one presented by equation 6 is relatively common. The mass flux is thus function of an accommodation coefficient , the interfacial area density ( ) and interface temperature assumed to be equal to the local cell temperature [10]. The interfacial area density is expressed as follows: 7 The bubbles number density and the bubble diameter are usually constants. The calibration of such models is done by the accommodation coefficient. Classical values for these parameters can be read in [10].
Vapor quality difference based model: This approach was proposed by Downar [13] and reused by Karathanassis [10] for flash flow nozzles, by Lorenzo [14] for transient flash flows, and by Palacz [15] for a CO2 ejector. It is called Homogeneous Relaxation Model (HRM) and its general form is: is the equilibrium vapor quality and a relaxation time. Then the final form of the mass flux is given by Downar :

9
The authors [13] called metastable liquid enthalpy; since the mass transfer does not affect the energy equation they used (same for [15]), it seems to correspond to inlet liquid enthalpy. However, one can read in [5] that liquid enthalpy may change during expansion for HRM model. So the liquid actual enthalpy is the more likely to be used the last being estimated by energy conservation. The relaxation time is defined by semi-empirical correlations presented in the original paper: 10 Here is the vapor volume fraction, is a non dimensional pressure difference; for pressures lower than 10bar it is defined as:

11
The constants in equation 10 are fitting coefficients.
The recommended values can be read in [13]. This formulation was named "homogeneous" because the two phases were supposed to have the same velocity but not the same temperature.
Temperature difference based model : The most notable model based on this difference is the so called thermal phase change model. The basic idea is to compute the heat flux associated to the phase change when mass transfer occurs. The mass flux is: 12 where is the interfacial heat transfer coefficient, is the interfacial area density (in m²/m3) and is the total enthalpy difference between the gas and liquid. The interfacial heat flux ( ) is applied also to the energy balance. This expression shows the relation between the capacity of the fluid to internally transfer heat from liquid to vapor and its ability to change phase. The liquid to saturation temperature difference is usually named superheat degree. This formulation was used by Wu [3]and Yixiang [8] for loose of coolant situations in nuclear power plant safety studies and by Bussac [16] who worked of two phase impulse turbines for geothermal power energy conversion. Regarding the heat transfer coefficient, Yixiang [8] proposes to use the Aleksandrov heat transfer coefficient that takes into account the gas/liquid velocity difference through the Péclet number. 13 where is the liquid's thermal conductivity and is the interfacial characteristic length. The Jakob number is: 4 where and are respectively the isobaric heat capacity, and the latent heat of the liquid. The Péclet number is: 15 where is the thermal diffusivity of liquid. In the cited Yixiang model, the heat transfer coefficient is multiplied by a correction factor in order to add a calibration parameter to the formulation. We could also cite in this section the work of Giacomelli [17].and Geng [18] who use a dimensionless temperature difference for a CO 2 and R134a two-phase ejector primary nozzle modeling. They used the same relaxation factor than in equation 6 multiplied by . This leads to similar behavior regarding the adjustment of the model. The interfacial area density can modeled on several ways. Based on the work of work of Wu [3], Yixiang [8] assumed in her work a flow containing equal sized spherical bubbles leading to the following expression:

16
The value of the bubble number density was used as an adjustment parameter by Yixiang [8].
Proposed model: The proposed model is based on the thermal mass transfer model but has a specific interfacial area formulation. The formulation is inspired by the one used by Wu [3] who proposed to model the interfacial area density at high vapor volume fractions based on a mist flow situation assumption using a number of droplets density ( ). The equation 16 was used for low values of vapor volume fraction ( ). For high values of ( ), the expression of for a homogenous sized droplet mist flow regime was [3]: 17 Nevertheless, Wu proposed a symmetric interfacial area density model stating that bubble and droplet numbers densities were identical and that maximal interfacial area density was attained at 0.3 vapor volume fraction. The model adjustment was done modifying .
Here and will not be assumed identical and a simple linear interpolation between and is proposed. Both and will be used for the model adjustment. This proposed model is called in the following lines TABD model. In Ansys-CFX, for Euler-Euler problems composed of two continuous phases, the interfacial area density is expressed as follows: 18 The interfacial length scale has to be set by the user. That can be a constant or a variable. It was here defined by a CEL expression (in CFX, a user defined function). In order to avoid getting non defined values of at values of 0 and 1, the volume fraction is capped to a certain minimum value; a value of 1e-7 was taken for this work. Note that liquid volume fraction is and and are null at 0 and 1 void fractions respectively.
Discussion on the models: The aim here is to study the calibration ability of the models regarding the nozzle flow rate and the outlet mixture velocity. Assuming that the first is mainly affected by the early evaporation close to the throat (low ) and the last is affected by the evaporation close to the outlet (high ), the flexibility of the models to adjust the evaporative profile at different stages of the expansion would be very useful for model calibration. Thus this analysis is focused on the sensitivity of the interfacial mass flux at both low and high vapor volume fractions. The limit between low and high is set to 0.5. So let's assume a linear expansion between inlet pressure of 5bar and an outlet at 1bar assuming that flashing inception begins at 3.5bar, a linear superheat reduction from 5K to 1K and a linear vapor volume fraction evolution from 0 to 1. These conditions are close to those of the test campaign conducted by [1]. This type of analysis omits 2D effects on the flow.  Figure 1. The results are shown in function of ; thus the beginning of the expansion is at and the end at For this example the parameters of the mass transfer models were: for HK model, , a=-0.25 and b=-2.24 for HRM model and =1e 9 and =1e 9 for TABD. The analysis was conducted as follows: first the nominal parameters were used to compute the mass flux for each pressure point of the assumed linear expansion. Then each parameter was changed by +/-10%; for models with various parameters, each one was modified at time. The positive and negative variations serve in showing an eventual non-linearity. Then the average mass fluxes relative variations ( ) were observed for and . The results are presented in Table 1, Table 2 and  Table 3 for HK, HRM and TABD models respectively. The tables give the variations applied to the parameters in bold text and the average variations observed on the mass flux in normal text. In the case of HK model the adjustment parameter variation affect linearly the both low and high void fraction fluxes. There is no difference between the two regions. So with such a formulation it is not possible adjusting the evaporation intensity at nozzle outlet without varying early evaporation intensity.
In the case of HRM, increasing affects inversely the flux; the same for 'a' and 'b' parameters; please note that 'a' and 'b' are applied to variables ranging between 0 and 1. The sensitivity has not the same magnitude for low and high values of except for . So in order to separately adjust the flux at low and high values, the parameters need to be modified in inverse magnitudes and try, thanks to the sensitivity magnitude differences, to modify the flux in the desired region of . This explains the HRM model calibration complexity mentioned before. Adj.
As a summary, it can be observed that models with a unique accommodation coefficient have limited adjustment flexibility. Besides, the HRM model contains three adaptation coefficients so it is flexible; however the sensitivities to the adaptations coefficients are not linear and are strongly coupled since it is not possible changing a parameter without acting on the effect (and the sensitivity) of the others. Finally, the proposed TABD model presents almost non coupled linear sensitivities making it easy to calibrate. Case study: In the early 90's a Japanese team [1] worked on waste heat recovery by impulse turbines using flash nozzles. Various nozzle geometries were tested. They were widely studied for a wide outlet pressure range. The researchers measured mass flow rate, thrust and pressure profiles. Here, only one nozzle geometry is explored. The dimensions and the operating conditions are shown in Figure 2 and Table  4 respectively.  Note that the mass flow in Table 4 does not depend on the outlet pressure; that means that the nozzle operates in critical conditions. The uncertainties of the measured data reported in [1] are within 2% for mass flow rate, 2% for thrust, 1% for pressure and 0.5% for temperature. This leads to an uncertainty of 2.8% for velocity by uncertainty propagation (knowing that thrust is equal to ). The outlet velocity is defined as the outlet mass fraction averaged mixture velocity.

Solver
configuration and adjustment procedure: The CFX solver is a coupled solver using a pseudo-transient formulation; the coupled option was selected for volume fraction as well. A bounded second-order upwind scheme was selected for advection. Please refer to Ansys CFX documentation for details on numerical resolution. A steady state simulation was done. The physical time step was set in a range between 1.10 -4 s and 1.10 -8 s depending on the boundary conditions. This parameter acts like a relaxation coefficient. The simulation was supposed converged when the mass and energy imbalance was lower than 0.5%, the inlet mass flow rate and the outlet velocity were steady; all residuals were in this situation lower than 1.10 -5 . The total energy formulation of the energy conservation equation was selected. The flow field was initialized at 0 vapor volume fraction, at inlet temperature and pressure and at 0m/s velocity. The kw-SST closure was used as turbulence model. The liquid and vapor properties of water were computed from standard IAPWS IF97 tables available in Ansys CFX. The liquid properties ( ) were computed as a function of temperature computed from the enthalpy resulting from the energy balance. Vapor was considered to be saturated. The drag force between phases was computed using a constant 0.44 coefficient after comparison to more elaborated models results. A mesh sensitivity analysis was done on the nozzle mass flow rate and outlet velocity for 1bar and 0.45bar outlet pressure and 1.5K degree of sub-cooling leading to the choice of a 12000 nodes mesh resulting from 15 radial subdivisions and a 2.10 -4 m of axial mesh size in the central part of the control volume.

Results:
Calibration or adjustment was done by first obtaining the measured mass flow rate adapting and then obtaining the outlet velocity changing for the highest inlet temperature and outlet pressure case (1.5 SC, and 1bar). Using the calibration parameters for this first case, all the cases were simulated. Then, for each case the calibration was done changing when changed and when changed. The calibration was stopped when discrepancy of the target variable was lower than 1%. CFD simulation results are presented as follows: the final calibration parameters are presented, and then the discrepancies before and after individual calibration are presented.
The final calibration parameters are presented in Table 5. Globally, when increasing sub-cooling, increases and decreases. When reducing the outlet pressure decreases. Mass flow rate changed slightly when changing however the average absolute discrepancy over the four outlet pressures was always lower than 0.7%. The effect of calibration on outlet velocity prediction is shown in Figure 3 where the red lines were obtained for calibration parameters of 1.5k SC and 6°2 8°3 . 5 14.7 114.5 50.8 7 1bar outlet pressure case and the blue lines represent data after calibration.

Figure 3 : outlet velocity discrepancies
The discrepancies where increasing with the inlet sub-cooling and where pressure dependent. To obtain the discrepancies in term of nozzle isentropic efficiency one has to multiply the values of Figure 3 by two ( ). After calibration all discrepancies were lower than 1% Adjustment parameters analysis: This section provides observations on the relationships between and parameters and some meaningful flow parameters. Figure 4 shows (in a log scale) as a function of inlet sub-cooling degree; an exponential relation cloud be used here as a first closure model; however further investigation needs to be operated in order to check why such a sudden increase of appears between 6K and 12K of sub-cooling. Various regions could exist needing a more complex closure model. In the original work of Wu [3], the calibrated values appeared to be related to the liquid superheat degree at inception point. However this is relatively complex to predict. Also, when plotting the calibration results obtained by Yixiang [8], a clear relation of with the inlet sub-cooling degree is observed. The values found here are of the same other of magnitude than the ones presented by Wu or Yixiang; however is these works, for a similar subcooling range the calibrated values ranged from 1.10 9 to 8.10 10 . Here ranges between 3.10 10 and 1.2.10 12 .This differences could be attributed to the difference in the interfacial area formulation and nozzle geometry difference. In order go further in this investigation supplementary geometries and operating conditions should be explored with the present model. Regarding the droplet number density, , it should be affected by parameters responsible of liquid ligaments fragmentation in the last part of expansion such as Weber or Reynolds numbers. After various parameters exploration, a relation appears between and the outlet phase averaged Reynolds number defined as follows:

19
Where is the volume fraction averaged mixture density, is the mass fraction averaged mixture velocity, is the volume fraction average mixture dynamic viscosity and is the outlet diameter. Figure 5 shows the calibrated values of as a function of outlet Reynolds number for the three inlet sub-cooling degrees.

CONCLUSION
This article presents a mass transfer model based on a non-symmetric interfacial area density formulation. A simplified sensitivity analysis showed a better flexibility and calibration ability of the presented model compared to literature models frequently used for CO 2 two phase ejectors modeling for example. This makes the simulation easily adjustable in terms of nozzle mass flow rate and outlet velocity. Adjustment parameters and appear to be correlated to inlet liquid sub-cooling degree and outlet two phase Reynolds number respectively. A closure model was proposed for ; for further investigation needs to be done using supplementary experimental data. Besides, the adjustment parameters need to be tested on other nozzle geometries in order to identify eventual geometry effects on interfacial area density formulation in the perspective of geometry optimization processes.