Numerical investigation of the effect of the shape of the nozzle on the flash boiling phenomenon

. In this paper analysis of influence of the nozzle shape on the flash boiling phenomena is presented. The mixture model was applied to account for two-phase flow through the nozzle. The Zwart-Gerber-Belamri model was used to describe the dynamics of the water-vapor phase change process, with the polynomial relationship linking the saturation pressure and temperature of the fluid. The effect of the shape of the nozzle on the mass flow rate of the two-phase mixture was numerically investigated by considering the stepwise and conical geometry of its divergent part. At first analysis of the influence of grid size as well as the ways of turbulence and near wall region modelling on the mass flow rate of two-phase mixture were studied for the diameter of nozzle neck of 0.42 mm. Predictions were found independent on these factors. Then simulations were carried out for three nozzle neck diameters (i.e., 0.62, 0.72 and 0.82 mm) and for several pressures (i.e., from 5 to 7 bar) and undercoolings (i.e., from 1 to 50 K) of water at the inlet. The shape of divergent part of the nozzle was found to not have effect on the mass flow rate of flashing water flowing through nozzle.


Introduction
When the subcooled or saturated fluid is discharged through the nozzle an abrupt pressure drop and a decrease of the saturation temperature below actual local liquid temperature occur resulting in the intensive boiling. This process is called the flash boiling flow. The phase change in flashing flows is thermally driven, since the rate of vapor generation depends on the heat transfer rate at the bubble interface. Flash boiling flows are found in many systems, e.g., in geothermal steam generators, in distillation devices, in two-phase ejectors, in fuel injection and atomization, etc. They are also observed during accidents in power plants or industry when hot pressurized fluids flow away through cracks in pipes or tanks. Therefore, accurate heat and mass transfer models of flashing flows are needed for safety and operability of many systems.
The flash boiling flow is characterized by intricate physics, which results in the complexity of mathematical and numerical models. In this phenomenon the heat transfer is strongly coupled with the mass transfer. Depending on the pressure difference and size of the opening in the nozzle heat transfer rates between phases in the flowing mixture may range from slow to very fast which significantly affects the way of modelling.
Numerous of approaches were proposed for modelling of different devices, which utilize flash boiling phenomenon. In the nineties of the last century Bilicki and Kestin [1] proposed meaningful modification of the classical homogeneous equilibrium model (HEM) by adding a single relaxation equation to describe the non-equilibrium interphase mass transfer between liquid and vapor phase. The new model called the homogeneous relaxation model (HRM) took into account the non-equilibrium evaporation, which is leading to the metastable conditions of the liquid phase, while the HEM assumed equilibrium vapor generation. In the next paper, Downar-Zapolski et al. [2] presented a correlation for the relaxation time which was applied to close the HRM. The correlation provided the characteristic time for the non-equilibrium evaporation process as a function of void fraction and non-dimensional pressure difference. The model proposed in [2] was implemented in a one-dimensional space (1D) and was successfully validated for the flashing water flow. In [3] the HRM was applied to study thermal non-equilibrium two-phase flows with flash boiling and condensation in short channel. Simulations were carried out for a twodimensional (2D) axisymmetric geometry. Obtained results were validated against experimental measurements. In the next work, Colarossi et al. [4] developed a 2D axisymmetric model of heat and fluid flow in the two-phase condensing ejector. The carbon dioxide (CO2) was the working medium in the considered device. The Eulerian pseudo-fluid approach was applied, i.e., the fluid was treated as being in a nonequilibrium thermodynamic state. The modified HRM was employed to account for phase transition. The model was implemented in the open-source OpenFOAM software. Simulations revealed that a pressure rise in the ejector was comparable to the experimental data. In turn, Smolka et al. [5] presented a model of a compressible transonic single-and two-phase flow of a real fluid through a heat pump ejector. In the energy equation they used the enthalpy as independent variable. Moreover, a thermodynamic and mechanical equilibrium between gaseous and liquid phases were assumed for the twophase flow, which meant that both phases had the same pressure, temperature, velocity, turbulence kinetic energy and turbulence dissipation rate. The proposed HEM was implemented in the commercial software ANSYS Fluent. The model was tested for a single-phase flow of the R141b and for a two-phase flow of the R744 fluid (CO2) in a three-dimensional (3D) geometry of the ejector motive nozzle. The 3D approach allowed to consider two configurations of the secondary flow, e.g., the tangent and radial directions. Successful validation procedure was also carried out by comparing obtained results with the experimental data for the primary and secondary mass flow rates in the ejector device as well as the pressure at the selected points in the computational domain. Next, Palacz et al. [6] proposed a more complex two-phase flow model of the CO2 ejector, i.e., the HRM which included metastable effects. Then they evaluated the accuracy of the new model against the previously proposed HEM [5] by comparing calculated results with the experimental data for the typical operating regimes of CO2-based refrigeration systems. The errors for the HRM predictions were approx. 5% lower than those of the HEM. Subsequently, the HRM from [6] was modified by applying the genetic algorithm optimization of the relaxation time correlation by Haida et al. [7]. Then the studies on the impact of the relaxation time definition on the CO2 two-phase flow behaviour as well as on the model accuracy were carried out. In the next paper, Battistoni et al. [8] compared a homogeneous mixture model in conjunction with the HRM for phase change with a multifluid model with the Rayleigh bubble dynamics model for phase change. These approaches were implemented in commercial packages CONVERGE CFD and AVL FIRE, respectively, and were applied for modelling of the cavitation in small nozzles, like those used in high-pressure diesel or gasoline fuel injectors. The compressible flow of three components, i.e., liquid, vapor and air were modelled. The obtained data were assessed against quantitative high-resolution experimental data. The amount of void calculated by the multifluid model was in a good agreement with measurements, while the mixture model overpredicted the amount of void. In turn, Janet et al. [9] tested three models for wall nucleation during flashing flow through converging-diverging nozzle. The thermoflow model was developed in the ANSYS CFX applying the Eulerian two-fluid model with phase change. They compared obtained results with the experimental data and found a good agreement with the critical flow rate and axial profiles, but the agreement with the radial void fraction data was not satisfactory. A 3D two-fluid model accounting for drag and non-drag forces as well as interphase heat transfer induced phase change during the flow of initially subcooled water in a vertical circular converging-diverging nozzle was developed by Liao and Lucas [10]. Interphase source terms in mass, momentum and energy conservation equations were modelled by constitutive relations. The model was implemented in the ANSYS CFX. Its reliability was confirmed by validation against the experimental data and results from a 1D driftflux model. Satisfactory prediction of the mass flow rate and cross-section averaged parameters was achieved. However, calculated radial distribution of void fraction was much flatter than the measured one. In succession, Karathanassis et al. [11] examined models of phase change mechanisms during the flashing flow through converging-diverging nozzle. Different approaches, i.e., the kinetic theory of gases (Hertz-Knudsen equation), HEM, bubble-dynamics considerations using the Zwart-Gerber-Belamri model (ZGBM), as well as HRM with semi-empirical correlations were tested. It was found that phase change models based on the kinetic theory of gases produced the most accurate predictions for all the cases investigated, while the validity of the HRM and ZGBM models was found to be situational. The multiphase model based on unsteady compressible Reynolds-Averaged Navier-Stokes (RANS) equations and on homogeneous mixture theory was proposed by Jin et al. [12] to study cavitating and flushing flows. Two different phase change mechanisms, i.e., one which accounted for pressure and the other for thermal effects were adopted. The results obtained showed an overall good agreement with the experimental data for cavitating and flushing flows. Next, Le et al. [13] developed a model of the flashing flow, which considered the nonequilibrium thermal effect. These effects were accounted for by using ad hoc modelling of the boiling delay. The two-phase mixture approach was assumed with the phase change which was depended on the difference between the vaporization pressure and the vapor partial pressure. The model was implemented with the aid of the ANSYS Fluent and was applied to simulate 2D axisymmetric flow of subcooled water through convergent-divergent nozzle. The numerical results showed a good agreement with the available experimental data in terms of the mass flow rate, average vapor fraction and static pressure along the nozzle. Moreover, sensitivity analysis was carried out concerning the grid independency, turbulence modelling approaches, near-wall treatment approaches, turbulence inlet parameters and semi-empirical coefficients. Recently, Lyras et al. [14] presented the compressible RANS solver FlashFOAM for calculating the phase change within various nozzle geometries undergoing rapid pressure drops. The solver was developed in the open source code OpenFOAM. The proposed model accounted for the interphase heat transfer by applying the HRM. The surface forces due to liquid-gas interfacial instabilities were modelled in a framework of a novel coupling of the HRM with the volume of fluid method. The novel formulation of the pressure equation was proposed. Model predictions were successfully compered with published experimental measurements, also for long nozzles. Moreover, new tests for flow characteristics and vapor generation in cryogenic liquid cases were carried out.
The main objective of this paper is analysis of the influence of the shape of the divergent part of the nozzle on the mass flow rate of water-vapor mixture flowing through it. In order to perform analysis, the model of the flashing flow of subcooled water was developed applying two-phase mixture model supplemented by the ZGBM to account for the phase change. Moreover, the sensitivity analysis of the solution on the grid size, turbulence modelling and near-wall treatment was carried out.
The article is organized as follows. First, the description of the considered nozzle geometries is presented. Second, mathematical and numerical models are described. Next, the sensitivity analysis is carried out. Lastly, mass flow rates calculated for two nozzle geometries and for three nozzle neck diameters (i.e., 0.62, 0.72 and 0.82 mm) as well as for several inlet pressures (i.e., from 5 to 7 bar) and undercoolings (i.e., from 1 to 50 K) of water at inlet are presented and discussed.

Geometries of nozzles
Two nozzle shapes, schematically presented in Fig. 1, were considered. Assumed shapes of nozzles corresponded to two possible ways of their manufacturing by applying the milling technique. The first nozzle (Fig. 1A) had the stepwise shape of the decompression zone, while the second one (Fig. 1B)  There is a question how these shapes of nozzles, particularly shapes of their divergent parts affect the mass flow rate of flashing water and what is the optimal shape of the nozzle.

Mixture model
For modelling of the unsteady turbulent flash boiling flow through the nozzle the two-phase mixture model was applied. Water was set as the primary (continuous, q) phase, while vapor as the secondary (dispersed, p) one with spherical shape of bubbles. The influence of gravitational acceleration was neglected. Considering above assumptions continuity, momentum and energy equations for the mixture model were the following: and the total energy, mass-averaged velocity, effective thermal conductivity, effective dynamic viscosity and density were calculated by applying the following formulae: where: subscript k denotes k-th phase, h -specific enthalpy, I -identity matrix, k -turbulence kinetic energy, p -pressure, t -time, a -volume fraction, Fviscous-dissipation function, l and lt -molecular and turbulent thermal conductivity, µ and µt -molecular and turbulent dynamic viscosity. Specific enthalpies of water and vapor were given by the following relations: where: cw is the specific heat of water, cp,v -specific heat at constant pressure of vapor, l -latent heat of vaporization. The reference temperature was equal to Tref = 273.15 K. The drift velocity which appeared in momentum equation, eq. (2), was calculated from the following formula:  (11) where ukq was the velocity of the k-th phase relative to the q-th phase. The mass fraction of the k-th phase was calculated from: The slip velocity, i.e., the relative velocity between secondary (p) and primary (q) phase, for turbulent flow was given by the following algebraic equation: (13) where: st is the Prandtl number (st = 0.75) and htturbulent diffusivity. The bubble relaxation time was given by: where dp is the bubble diameter (dp = 10 -6 m). The drag function was defined by applying the following Schiller and Naumann model: The acceleration of bubbles of the secondary phase was calculated from: The continuum surface force model was applied to account for surface tension effect. In this model the following source term for two phases flow was added to momentum equation, eq.

Numerical implementation
The unsteady thermo-fluid model of the flashing boiling flow through the nozzle, which was based on the mixture model, was implemented with the aid of the ANSYS Workbench environment. Geometries of nozzles were prepared in the ANSYS DesignModeler, grids were generated in the ANSYS Meshing, while simulations were carried out in the ANSYS Fluent. In sensitivity analysis stepwise and conical nozzles with the necking diameter of D2 = 0.42 mm were used. Then six meshes were generated for each nozzle (see meshes no. S1 to S6 and C1 to C6 in Table 1). For the parametric analysis of influence of the inlet pressure and undercooling, neckings of diameter of D2 = 0.62, 0.72 and 0.82 mm were applied (see meshes no. S7 to S9 and C7 to C9 in Table 1). Meshes parameters, i.e., total number of elements, number of elements along the necking radius, maximal value of the skewness and aspect ratio are given in Table 1. Geometries were specially prepared for meshing by dividing them into rectangular/trapezoid and triangular subdomains. Then rectangular and triangular elements were generated, respectively, in these subdomains. Rectangular elements dominated, while the triangular ones were generated only where it was required due to mesh quality reasons. Meshes were refined close to walls, especially close to the wall of the necking (i.e., bias factor along the radius of the necking was equal to 8). Exemplary meshes for necking regions are presented in Fig. 2  Simulations were carried out applying the segregated (pressure based) solver which decouples the solution of governing equations and is usually used for incompressible flows. However, in considered flashing flows compressibility effects were accounted for by using compressible water and vapor models. Solutions were obtained using the SIMPLE pressure-velocity coupling scheme, the second order implicit temporal discretization, second order upwind spatial scheme for density, momentum, energy and turbulence equations, the QUICK spatial scheme for void fraction equation, the PRESTO! spatial scheme for the pressure equation and the least squares cell based gradient reconstruction method. Moreover, underrelaxation factors were adjusted referring to the stability of the solution. Residuals were set to 10 -5 for all equations. The time step size was equal to 10 -6 s, which allowed for stable calculations. Simulations were conducted for 0.02 s, which ensured obtaining of steady mass flow rates through nozzles, i.e., asymptotic values which did not varied were attained.

Turbulence modelling
The basic turbulence model was the realizable k-e model with the scalable wall function (SWF) to account for presence of boundary layer. Moreover, in the sensitivity analysis two additional approaches were tested, i.e., the realizable k-e with the enhanced wall treatment (EWT) and the SST k-w model. Compressibility effects were accounted for in the all turbulence models. For other parameters and options standard settings of the ANSYS Fluent were applied. The turbulent dynamic viscosity and turbulent thermal conductivity, which appeared in eq. (6) and (7), were calculated by applying formulae relevant for used turbulence model.

Boundary and initial conditions
At the inlet to the computational domain three value of the inlet pressure of the subcooled water were applied, i.e., pin = 5, 6 and 7 bar. Moreover, eleven values of water undercooling were used, i.e., DT = 1, 2, 4, 6, 8, 10, 15, 20, 30, 40, 50 K. The turbulence intensity was set to 2%. At the outlet from the nozzle, the fluid discharged to the surroundings at psur = 1 bar. Walls of the nozzle were adiabatic. In the initial state the nozzle was filled with motionless water kept in inlet conditions.

Fluid properties
Fluids properties required in the simulations were either taken from the database or calculated using models of properties available in the ANSYS Fluent. For water the compressible liquid density model was used, while the other properties were constant. For vapor the Soave-Redlich-Kwong real gas model and specific heat defined by piecewise polynomial function were used, while the other properties were constant.

Sensitivity analysis
The sensitivity analysis was carried out to check the solution independence on the grid size, turbulence model and the way of modelling of turbulence in the boundary layer. In this analysis the necking diameter was set equal to D2 = 0.42 mm and inlet pressure and undercooling of water were equal to 6 bar and 10 K. Simulations were preformed for six meshes which were generated for each nozzle shape, i.e., meshes no. S1 to S6 and C1 to C6 for stepwise and conical nozzle, respectively. Three approaches to turbulence modelling were applied, i.e., the realizable k-e with SWF, realizable k-e with EWT and SST k-w. The first turbulence model was applied for all meshes, while the second and third one for three most dense meshes only due to required mesh refinement close to walls. Results, i.e., steady values of the mass flow rate flowing through the nozzle m & flash, are presented in Table 2 for the all analysed cases. Relative differences between obtained results for the smallest and largest mesh were small and did not exceed 3.5, 1.8 and 1.7% for the realizable k-e SWF, realizable k-e EWT and SST k-w, respectively. The lowest relative differences were found between results obtained by applying the realizable k-e SWF and SST k-w, while the highest between realizable k-e EWT and SST k-w. They did not exceed 0.7 and 1.9%, respectively. Low values of relative differences between predictions suggests that obtained results are independent on the grid size and turbulence modelling. The lack of influence of the turbulence modelling on the m & flash suggests that flashing flow behaviour and m & flash are not affected by the turbulence phenomena in the boundary layer and bulk stream. Moreover, for corresponding cases for stepwise and conical nozzle the same values of the steady mass flow rate of flashing water were calculated, but this is analysed in detail in the next section.

Analysis of influence of the nozzle shape
In order to find relation between the mass flow rate of two-phase mixture flowing through the nozzle and the shape of the divergent part of the nozzle two different nozzle geometries, i.e., stepwise and conical, were studied. Simulations were preformed for three diameters of the nozzle neck, i.e., D2 = 0.62, 0.72 and 0.82 mm, for three inlet pressures of water, i.e., 5, 6 and 7 bar, and for eleven inlet undercoolings of water, i.e., 1, 2, 4, 6, 8, 10, 15, 20, 30, 40, 50 K. Meshes used for simulations were following: S7, S8 and S9 as well as C7, C8 and C9 for the stepwise as well as conical nozzle, respectively -see Table 1. As the sensitivity analysis revealed grid size independent solution, applied grids were not very dense and correspond to grids S2 and C2. Results of simulations, i.e., steady values of the mass flow rate of flashing water flowing through the nozzle m & flash, are summarized in the Table 3. The shape of the expansion part of the nozzle had no influence on the m & flash. Most of obtained values differed on the third decimal place. However, increase of the m & flash with the rise in the necking diameter as well as with the increase of inlet pressure and undercooling of water were observed. This suggest that behaviour of the flash boiling flow depends on the pressure difference between inlet and outlet to the nozzle as well as initial undercooling of the water. Moreover, more important is the necking shape (e.g., length and diameter) than the shape of convergent and divergent part of the nozzle. and was used for analysis of two shapes of the nozzle, particularly their divergent parts. The stepwise and conical nozzle shapes were considered. In the first step the sensitivity analysis was carried out. Influences of the grid size, turbulence modelling and near-wall treatment on obtained results were investigated. It was found that the model predicted independent values of the mass flow rate of flashing water flowing through the nozzle. This suggested that that flashing flow behaviour was not affected by the turbulence phenomena in the boundary layer and bulk stream. In the second step influence of the nozzle shape was examined. Three nozzle necks diameters (i.e., from 0.62 to 0.82 mm) as well as three pressures (i.e., from 5 to 7 bar) and eleven undercoolings (i.e., from 1 to 50 K) of water at inlet were considered. The shape of the expansion part of the nozzle had no influence on the mass flow rate of flashing water. For both geometries mass flow rates were the same.