Numerical Simulation of Flash-Boiling Through Sharp-Edged Orifices

In cases where a liquid, stored in a high-pressure vessel, flows towards a low-pressure environment, some interesting phenomena might be triggered that can change the morphology of the flow. If the liquid’s local pressure drops below the local saturated vapour pressure, boiling process, generally termed flashing, occurs. This phenomenon has been exploited in clean combustion technologies to improve mixing and combustion processes, and multiple aerosol industrial applications where small droplets in the vicinity of the nozzle exit are desired. Flashing is known to produce fine sprays and like cavitation has an impact on the atomisation and spray dynamics [1]. The phase change is manifested by bubbles forming within the liquid, changing progressively the regime of the flow [2]. The exact time and position of the bubble generation is not a trivial task especially in the presence of turbulence. The flow is strongly dependent on the initial pressure and superheat or subcooling degree. The flow characteristics dominantly affect the subsequent atomisation upstream of the nozzle exit. The regime of the jet varies depending on the local topology including the length of the tubes, size and shape of the orifice.

In flashing, extensive bubble nucleation occurs in the liquid. Thus, multiphase CFD codes dealing with flash boiling need to be able to account for the impact of the phase change throughout the liquid. The growth of a single bubble was firstly described by Rayleigh [3] and Plesset [4]. Examples of incorporating the Rayleigh–Plesset model in CFD methods for calculating bubble growth and collapse can be found in Kawano et al. [5] and Giannadakis et al. [6]. Some studies assumed thermodynamic equilibrium using the homogeneous equilibrium model (HEM). The HEM is a simple model capable of simulating specific flows. It is popular among one-fluid Eulerian approaches and in long pipes where there is sufficient time to reach thermodynamic equilibrium [2, 7]. For sharp-edged orifices, the flow is not expected to have time to reach equilibrium. The liquid’s state delays the initiation of the boiling process, exhibiting a metastable condition. The modelling of this delay was studied in Ref. [8]. A similar formulation is used in this present study, employing the homogeneous relaxation model (HRM) developed by Downar-Zapolski et al. [9]. HRM has previously been success- fully used for laminar simulations for flashing water in Ref. [10–12]. Usually, the simulations considered large nozzles, with relatively large ratios of nozzle length to nozzle diameter greater than four. Turbulence which might be triggered across the flow could either be sustained or damped [13].

In the present study, the Reynolds-Averaged-Navier-Stokes (RANS) approach is used. Two equations for the vapour mass fraction and the presence of the ambient air in the mixture have been added to the generic compressible solver within the open source computational fluid dynamics (CFD) code OpenFOAM [14] to facilitate the simulation of multiphase flows. The source term of the vapour mass fraction is developed within the solver according to Schmidt et al. [11]. Additionally, the pressure–velocity coupling follows the same implementation as developed and described in [10, 11] including the contribution of phase change into the PIMPLE algorithm. RANS simulations were performed to simulate the two-phase flow within the orifices, calculating how the primitive variables change with respect to the storage conditions. Even along the small length of an orifice, interfacial heat transfer is observed, and the discharging liquid mass transforms gradually from a continuous liquid phase to bubbly with separate vapour bubbles. This same mechanism will transform the bubbly flow into a continuous gas phase with separate liquid droplets prevailing the whole flash-boiling atomisation process.

The classical Navier–Stokes equations are solved in a fully Eulerian one-fluid approach, which solves one set of equations for the whole flow field. The liquid and vapour phases are considered to have the same velocity. The presence of ambient air in the atomisation region and its influence in the calculation of density are considered. The equations of continuity, momentum, energy and vapour mass fraction are written in a compressible formulation as:

Here, $\rho, u_i, p$ and $h$ are, respectively, the mixture density, velocity vector, pressure and enthalpy. The variable $x$ is the vapour mass fraction (quality) and is used to calculate the mixture composition. The mixture properties such as density $\rho$, enthalpy $h$, viscosity $\mu$ are, in case of pure liquid-vapour mixture, the weighted average value of the liquid and vapour properties, using the quality as a weight function, and $\alpha_{\text {eff }}$ is the effective thermal diffusivity of the mixture. In Newtonian fluids the deviatoric viscous stress tensor is $\tau_{i j}=2 \mu S_{i j}-\frac{2}{3} \mu \delta_{i j} S_{k k}$, where the strain rate tensor is calculated as $S_{i j}=1 / 2\left(\partial u_i / \partial x_j+\partial u_j / \partial x_i\right)$.

The source term Γ in eqn (4), stands for the interfacial mass transfer. It is closely related to vapour generation and flashing. An additional relationship is needed to close the above set of transport equations. The HRM implemented here, assumes that the instantaneous vapour mass fraction x relaxes at the local equilibrium value, $\bar{x}$ over a timescale Q following [9, 11]. This first-order approximation of the term G is given as follows:

Based on the experimental findings in Ref. [15] for superheated liquids, the isentropic assumption for the equilibrium quality $\bar{x}$ may lead to unrealistic predictions. Thus, the equilibrium quality is calculated assuming isenthalpic process,

Here, the liquid and vapour enthalpies denoted as $h_{l, \text { sat }}$ and $h_{v, \text { sat }}$ are calculated at the saturation conditions. Knowing the void fraction $\alpha=\left(\rho_l-\rho\right) /\left(\rho_l-\rho_v\right)$, the quality can be directly computed as:

The timescale $\Theta$ is calculated using the following relationship [9]:

Where $\Theta_0$ is a constant with time dimensions and is equal to $3.84 \cdot 10^{-7}[\mathrm{~s}]$, and $\varphi$ is the dimensionless pressure given by:

Equations (8) and (9) were derived from experiments for flashing water and are valid for pressures higher than 10 bar.

The standard $k-\varepsilon$ model is used. Preliminary computations assuming laminar conditions were also performed but showed in general less accuracy and demanded more time to reach steady state for the current simulations (typically, steady state was achieved after 4ms for RANS).

The segregated approach is used as a general framework for the present work. The equations are solved in a finite volume approach. The presented method was developed as a dedicated new solver HRMSonicELSAFOAM within the framework of OpenFOAM, which gives it the ability to handle arbitrary shape unstructured polyhedral meshes in two and three dimensions with all field variables stored in the centre of the control volumes and the interpolated fluxes are stored at the cell faces. Since the developed code is a pressure based solver, the PIMPLE algorithm is used for the pressure–velocity coupling. Up to three PIMPLE iterations are used for every simulation. The pressure equation is solved using a biconjugate gradient stabilised method that offers smoother convergence than the traditional biconjugate gradient methods (PBiCG). For the rest of the variables, PBiCG is applied. The solver uses a run-time selectable preconditioner. In the cases presented here, a diagonal-based incomplete LU smoother suitable for asymmetric matrices is used. The fluxes are evaluated using a gamma TVD scheme. Also, a second-order scheme that uses the least square method is selected for the gradient of pressure and the Laplacian is evaluated using a Gauss linear scheme.

Numerical simulations concerned a series of NASA experiments [16] involving liquid nitrogen flowing through sharp-edged orifices. The schematics of the experiments are shown in Fig. 1. Liquid nitrogen was released from a high-pressure storage vessel with pressures above nitrogen saturation point. The flows are simulated as axisymmetric. This experimental work is one of the very few flashing experiments where two-phase liquid nitrogen is the working fluid and provided the dataset to test the model for different cryogenic scenarios rather than the trivial water channel flows usually considered in other experiments. The nozzle diameter was D = 0.358 cm and L/D = 0.19. The small length to nozzle diameter ratio was kept the same for all the test cases considered in the present study. The tests covered a wide range of temperatures, from 0.75 to 1.035 times the critical temperature (126.3 K), and the pressures were slightly above and up to two times the critical pressure (33.958 bar). In Table 1 the inlet and outlet conditions are shown. Three series of experiments corresponding to different initial temperatures are considered here and the measured properties are denoted by in and out for inlet and outlet respectively. Due to the pressure difference on both sides of the orifice, the liquid flows, and at some point, the liquid pressure drops below the saturation pressure initiating the flashing process.

The maximum flow rates were measured at both upstream and downstream orifices with the difference between the two found to be less than 1.75% for all the experiments. Here, the upstream data are used for comparison.

Figure 2 shows the mass flow rate per unit area for constant stagnation temperature equal to approximately 95.1 K and various inlet pressures. The temperature at the outlet varied from 83.2 K for the lowest stagnation pressure at 8 bar up to 91.7 K for the maximum stagnation pressure of 67.2 bar. The mass flux was calculated from the mixture density averaged by the vapour mass fraction and velocity. The mass flux is minimum at the lowest pressure and increases in a rather parabolic manner up to the critical pressure. This dependency should not be attributed to choking which is less likely to occur in small tubes and L/D without any strict limits. The experimental observations were in line with this argument at least for subcooled stagnation conditions. The predicted mass fluxes linearly increase for inlet pressures greater than the critical pressure and achieve the highest value at the highest inlet pressure 67.2 bar. Similar trends occur when increasing the inlet temperature to 110.1 K (Fig. 3) and 119.4 K (Fig. 4). The mass flow rates show an increasing pattern by raising the inlet pressure with a linear behaviour as before for pressures above the critical one. The inlet temperature is shown to have a significant impact on the dynamics of the jet and consequently the mass flow rate. The difference between the CFD results and the experimental data remained in most of the cases investigated less than 5% and always below 8.5%. For the same pressure and geometry, the flow rate decreases as the inlet temperature increases. Higher stagnation temperature manifests the density to be more sensitive to pressure change and smaller average density, thus smaller mass flux. As mentioned above, eqns (8) and (9) for the HRM were initially proposed for water. The good agreement achieved here suggests that the formulations are also suitable for jets of cryogenic nitrogen. The correlation for the timescale Θ influences the vapour mass production rate as shown in eqn (5).

The better predictions of the mass flux for lower inlet temperatures compared to the ones for $T_{\text {in }}=119.4 \mathrm{~K}$ might be an indication of the need for alternative correlations for liquid nitrogen. Additionally, simulations using a different correlation for $\Theta$ (valid for pressures below 10 bar for water, proposed in [9]) resulted in better mass flux predictions in some cases for low pressures at $T_{\text {in }}=119.4 \mathrm{~K}$. The flow is highly turbulent with $O(R e) \sim 10^6$. In all cases, the Reynolds number increases with the increase of the inlet pressure. Simulations' results presented in Figs. 3, 4 and 5 show that the $k-\varepsilon$ can accurately simulate flows of cryogenic jets through sharp-edged orifices.

As mentioned before, flashing jets of cryogenic liquids are usually assigned to metastability and are subject to nucleation at different stages and regimes. Here, the metastability refers to the phase change that drastically occurs to the liquid jet. Studies in nozzles indicate that the initiation of this phase change is possible inside the nozzle [17] and is usually pronounced in long nozzles (large L/D). For short nozzles or orifices uncertainties may exist implying that geometry might have less influence. As reported in the NASA experiment [16] considered in the present study, the flashing jet vaporises within the orifice leading to an alteration of the flow regime and the measured flow coefficient. In the present numerical simulations, the metastable jet exists inside the orifice and phase change starts at the tip of the inlet corners. This phase change, which results in vapour generation, is also reported in the literature for different geometries [e.g. 18] and seems that does not dependent on the working (Newtonian) fluid. Although, for outlet pressures greater than the saturation pressure, a single-phase jet is expected according to the experiment and for pressures lower than the saturation pressure the jet vaporises. It is interesting to track the velocity inside the orifice and examine how it changes until it emerges in the low-pressure region. In Fig. 5 the axial velocity is shown at the upstream and downstream orifices for different inlet temperatures at approximately the same pressure. The results shown correspond to half of the domain (symmetric to the flow stream axis). The flow along the orifice radius R (R = 0.179 cm) follows a similar trend in all cases. The upstream velocities are very close to each other and have the maximum value at the centreline of the jet (r/R = 0). As soon as the jet flows towards the exit, nucleation takes place decreasing the mixture density, and higher velocity occurs with a maximum value at approximately 70% of the orifice radius. The velocity finally minimises to zero due to the no-slip condition. In Fig. 6 the full domain is considered offering a general overview of the jet at the downstream orifice. Keeping the inlet temperature constant and varying the pressure, the axial velocity changes correspondingly giving a parabolic-like profile. The peaks of the velocity indicate that vapour is present in the mixture along with the radial direction.

HRMSonicELSAFoam is developed from a generic compressible solver in OpenFOAM tool- box for simulating flashing jets. The model takes into consideration the rapid pressure drop and delay in the phase change which characterises flash boiling and interfacial heat transfer. The current solver is based on the RANS approach. Validations have been conducted using the NASA experiments [16] involving liquid nitrogen jets flowing within short-edged orifices for a large range of pressures and different inlet temperatures. The flash-boiling model employs a semi-empirical formulation, derived for water but proved to be applicable for cryogenic nitrogen as well. The results for the mass fluxes are found to be in reasonably good agreement with the experiments. The manifested boiling that is observed in the experiments is also captured by the numerical simulations. The accurate predictions for mass flux, imply that the velocity, mixture density and other fundamental properties needed for simulating the atomisation pro- cess can also be accurately computed despite the metastable jet inside the orifice that vaporises changing the flow regime. Subject to further validations studies, the approach developed here should also be applicable to both subcooled and superheated liquids that are subject to sudden depressurization through nozzles and offer a good understanding of the thermodynamic changes that happen, providing insight to the primary atomisation.