Toward the Simulation of Flashing Cryogenic Liquids by a Fully Compressible Volume of Fluid Solver

: We present a fully compressible single-ﬂuid volume of ﬂuid (VOF) solver with phase change for high-speed ﬂows, where the atomization of the liquid can occur either by the aerodynamics or by the effect of the local pressure. The VOF approximation among a non-miscible phase (non-condensable gas) and a mixture of two ﬂuids (liquid and vapor) represents the liquid core of the jet and its atomization. A barotropic model is used in combination with the equation of state (EoS) to link the mixture density to pressure and temperature. The solver is written with the aim to simulate high-pressure injection in gas–liquid systems, where the pressure of the liquid is great enough to cause signiﬁcant compression of the surrounding gas. Being designed in an C++ object-oriented fashion, the solver is able to support any kind of EoS; the aim is to apply it to the simulation of the injection of liquid propellant in rocket engines. The present work includes the base development; a veriﬁcation assessment of the code is provided by the solution of a set of numerical experiments to prove the boundedness, convergence and accuracy of the method. Experimental measurements of a cavitating microscopic in-nozzle ﬂow, available in the literature, are ﬁnally used for a ﬁrst validation with phase change.


Introduction
Increasing the injection pressure in aerospace propulsion systems is a way to improve the combustion efficiency and to reduce CO 2 emissions. Over the years, this strategy has successfully been applied to jet engines, rocket engines [1][2][3][4], and piston engines [5,6]. On the other hand, the design of an injection system for stable and efficient combustion requires a thorough understanding of fluid behavior at the conditions it operates. The complex physics behind the liquid atomization process makes its study extremely difficult. Phase-change therefore plays an important role in the atomization process [7][8][9][10]: it can help to achieve a finer atomization, but also reduce the spray stability, and it can also damage the injectors and decrease its reliability and durability. Liquid density variation during phase change strongly affects the fuel properties [11]. These conflicting attributes have spurred a renewed interest in understanding the complicated flow physics inside these devices [5,12]. Cavitation is a classical example of a multiphase and multiscale problem in fluid dynamics [13,14]. This problem becomes largely complex if cryogenic fluids are involved, due to the thermodynamics effect on the vaporization process [15][16][17], the detailed mechanics of the vaporization process, the effect of unsteadiness, the energy exchange between gas and liquid phases with phase change, and because the growth or detachment of the cloud cavity in cryogenic fluids remain a challenge for fluid dynamics researchers. Several methods have been developed and tested in the last decade to model phase change and cavitation [18][19][20][21][22][23][24]: examples are scale resolving or direct numerical simulations (DNS) performed on thousand of nodes [25] or to solve atomization [26]. DNS remains an interesting option but with very high-cost demands. Cavitation modeling in high-pressure injection and cryogenic fluids involves the simulation of multiphase, in the context of this work liquid and gas, and multicomponent, (i.e., several instance of the same phase) immersed into a turbulent flow with changes in the thermodynamic properties and density. Different research studies show that unsteady Reynolds-averaged Navier-Stokes (URANS) simulations underestimate the formation and the extent of cavitation because they over-predict the turbulent viscosity in the cavitating zones [27][28][29][30]. A comparison between results from URANS and large eddy simulation (LES) applied to injection and cavitation modeling [31] shows that URANS fails to predict the incipient cavitation at small injection pressures, while LES turbulence modeling better captures the effect at several turbulent scales characterizing the phase change. Phase change and cavitation have been studied by the meshless lattice Boltzmann method (LBM) [32] and by smoothed-particle hydrodynamics (SPH) [33] also. However, these methods suffer from some drawbacks as they are applied to cavitation modeling: they are computationally expensive and the arbitrary distribution of the particle in the domain has a severe impact on the solution. In the framework of the finite volume method (FVM), conventional mesh methods are largely used for the multiphase problem and are grouped in two families, namely the multi-fluid and the single-fluid approach.
In the most general multi-fluid method, there is no mechanical equilibrium (nonzero slip velocity) among the phases, and each phase has its own velocity, pressure and temperature. A generalization of the classic two-fluid approach [34] was successfully used in several applications [35][36][37]; the method requires accurate models for the mass, momentum and energy transfer through the interfaces. The method is available for compressible phases also [38,39]. The single-fluid approach accounts for one mixture velocity and single-fluid pressure and temperature. In the single-fluid approach [20], interfacial models are not required, and the mass transfer among the phases can be included in several ways. Widely used approaches are the level set method [40], the coupled level set-volume of fluid (CLSVoF) [41][42][43] and the VOF interface-capturing method by Hirt et al. [44]. A detailed review of such models is reported in [45], and also represented in Figure 1; they were originally developed for incompressible flows and then extended to compressible cases [46][47][48]. A critical aspect in these formulations is the coupling of the phase-fraction equations with submodels that account for the phase change. In many cases, these models are based on the Rayleigh-Plesset equation [49][50][51][52] and are included into a source term in the phase fraction equations [22]. The Rayleigh-Plesset equation describes the oscillations of a spherical bubble wall under the assumption that the fluid surrounding the bubble is incompressible. An alternative method to handle cavitation with fully compressible VOF solvers consists of using barotropic models. Barotropic models ensure that pressure and density are linked to satisfy the liquid and vapor equations of state (EoS) [20,53], despite the fact that they cannot reproduce the baroclinity term (∇ρ ∧ ∇p)/ρ 2 since the density variation is aligned with the pressure variation [54,55]. An example of the barotropic model is the homogeneous equilibrium model (HEM) that is based on the assumption of the mechanical and thermodynamic equilibrium between the liquid and vapor phase. While this assumption is not general, it has proven to work in many applications [31,[56][57][58][59]. Finally, in the simulation of high-pressure liquid injection, it is important to track the evolution of the non-condensable gases (NCG) in the nozzle, to simulate liquid atomization in the presence of swirl cavitation and, eventually, of hydraulic flip [5,21]. As a result, a three-phase system must be solved [20,22]. Several works have been published in the literature, where solvers, including separate tracking of NCG, were used [10,20,22,23,60], To the authors' knowledge, very few works in the literature combine a barotropic model for phase change with a compressible VOF solver for studying the in-nozzle effects on atomization [15,19,20]. These models differ in the way that the phase change is resolved.  Figure 1. Available methods for the solution of two-phase problems. In this work, an algebraic-type VoF is applied to capture the interface.

Motivation of This Research
In the field of chemical rocket propulsion, oxygen and hydrogen are preferred fuels because of their high specific impulse. In order to minimize the rocket fuel tank structure, oxygen and hydrogen are liquefied at a very low temperature, thereby leading to cryogenic combustion. Two-phase flows resulting from the atomization of liquid jets play a significant role in the proper functioning of cryogenic liquid-propellant rocket engines under subcritical operating conditions. The aim of this research is to work on the development of a framework to simulate liquid atomization in rocket engines, thereby being able to account for the combined effects of both the mechanical breakup and the thermodynamic flashing. The multiscale nature of the physical problem and the typical operating conditions in rocket engines make experimental investigations complex and expensive; for this reason, numerical tools represent a valid alternative to such a study.

Goals and Highlights
The aim of this work is to set up a strategy including the combination a fully compressible single-fluid VOF solver for gas-liquid systems, where the non-condensable gas varies its density for the effect of the liquid pressure, and the heat transfer and a run-time selectable phase change model is used to model either the phase-change induced by the thermodynamics and the mechanical breakup in the primary atomization. The presented solver, developed in the OpenFOAM Technology [61] handles the compressibility effects within the VOF approach by decomposing the phase density equations into an incompressible part and a compressibility correction to distinguish different phase densities and their variation. Different kinds of thermodynamics may be applied; in the present work, phase change in the validation tests is modeled following the HRM (homogeneous relaxation model) that is coupled to the equation of state, to link the density, pressure and temperature of the three-phase mixture. Thanks to its object-oriented structure, the solver supports any kind of equation of state (perfect gas, Peng-Robinson, tabulated) that can be defined at runtime; alternatively, the thermodynamic library CoolProp is coupled to the solver as a dynamic library to be called during the simulation. As a result, the proposed implementation is general and potentially covers a wide range of applications. The extension to a temperature dependence in the EoS can give the chance to formulate a more general and versatile approach. The discretization of the governing equations is based on the finite volume approach. Mass, momentum and energy are solved in a segregated fashion, using the pressure-implicit split-operator (PISO) algorithm [62]. The solution of the momentum, energy and the mass conservation is employed by a PISO algorithm [63,64].

Paper Structure
The paper is organized as follows: the theory of the solver, the discretized solution of the phase fraction equations in the presence of phase change, and its coupling to the segregated solution of the governing equations is presented in Sections 2-5. The solution algorithm is described in Section 6. Code verification is performed on one-dimensional numerical experiments, and it is presented in Section 7, while validation on a threedimensional internal nozzle flow case is reported in Section 8. The main conclusions are drawn in Section 9.

Compressible VOF Solver with Phase Change
An algebraic-type VOF method belonging to the family of interface-capturing methods ( Figure 1) is applied to capture the interface [44] between non-condensable gases and the cavitating mixture of liquid and fuel vapor, which is transported as a single phase. The system is therefore treated as having two components (cavitating liquid-vapor and gas) and three phases (liquid, vapor and non-condensable gases). The phase change between the liquid and fuel-vapor is treated using the barotropic equation of state. In the resulting two-phase system of cavitating fluid and non-condensable gases, each phase has a partial volume V i , that is a fraction of the volume V of the control volume (V i ⊆ V), and it is defined by its local volume fraction α i ∈ [ 0, 1 ]: being the subscript i = lv used for the cavitating mixture and i = nc for the non-condensable gas (nc). The following constraint applies: The single-fluid approximation leads to the so-called mixture variables, namely the following: -Mixture velocity: -Mixture density: ρ = α lv ρ lv + α nc ρ nc (4) -Mixture viscosity: -Mixture heat capacity: -Mixture thermal diffusivity: The phase fraction equations for the multiphase compressible cavitating flow are written as a system of two equations: that has to obey the compatibility condition Equation (2). From the conservation of the phase fractions, it follows that Dα nc Dt = − α nc ρ nc Dρ nc Dt (9) and Equations (9) and (10) can be combined to obtain The LHS of Equation (11) is the total derivative of the compatibility constraint Equation (2): so Equation (11) can be written as Equation (13) is the continuity equation in the compressible form; terms on the RHS of Equation (13) account for flow compressibility and include the density variations of each phase as a function of the temperature and pressure. The flow velocity is calculated from the momentum equation: where p is the fluid pressure, τ is the stress tensor and −g · x∇ρ is the buoyancy term, whose effect is quite negligible in the applications discussed in this work; the last term in Equation (14) represents the source momentum due to the surface tension on the surface. It can be seen as a surface integral on the interface surface: The term σ is the fluid surface tension in [N/m], andn is the unit vector normal to the liquid interface, whose center is located in x s . κ is the interface curvature [m −1 ] defined as S f is the cell faces surface area vector defined as the scalar product between the cell faces normal and the cell-face area. The normal unit vector is defined, taking into account the gradient of the gas phase:n = ∇α nc |∇α nc | This source term is non-zero only at the interface where the Dirac function δ(x − x s ) is active for x = x s . In this approach, the surface tension between the liquid-vapor interface is neglected; this assumption leads to a diffuse interface, which can be acceptable in several applications [20] and is acceptable if vapor and liquid are considered miscible fluids.
The modeling of the surface tension between the air and the cavitating fluid is reached following the same procedure described in De Villiers et al. [65]. The integral defined in Equation (15) cannot be calculated directly; hence, the widely used continuum surface force model (CSF) of Brackbill et al. [66] is applied to overcome this problem: The CSF model may potentially lead to physically unrealistic velocities at the interface [67]; on the other hand, it is a reliable method to account for tangential stresses due to variable surface tension (i.e., Marangoni effect). Although the continuum-surfacestress (CSS) approximation [68] and the ghost fluid method (GFM) [69] can better handle sharp transitions of the surface density and capillary forces, they are computationally more expensive.

Phase-Fraction Equations
From the combination of with Equation (9), an intermediate form of the phase fraction equations can be derived: If Equations (13) and (20) are combined, it follows:

Dρ nc Dt
Adding and subtracting α nc (∇ · u) to Equation (21) and recalling Equation (19), the final form of the phase fraction equation in the compressible form is recovered: The term dg/dt on the RHS of Equation (22) determines the compressibility effects; if dg/dt = 0, Equation (22) is the formulation of a variable density solver.
Being that it follows The compressibility term and α nc (∇ · u) present in Equation (22) are solved in a semi-explicit fashion. The resulting form of the phase fraction equations, is, therefore, No source terms are present on the RHS of Equation (25) because the phase change is handled by the barotropic equation of state. With respect to Equation (22), an artificial interface compression term is added to the LHS of (25); this is similar to what is performed for the scalar-flux second-moment closure in combustion modeling, where a "countergradient" transport is used to model the dynamic of turbulent flames [70] (see also [21,22]). A common closure used for counter-gradient transport has the form where U c is the compression velocity at the interface, and to preserve boundedness [71], it is proportional to α i (1 − α i ). C α is the compression coefficient to ensure interface sharpening that is usually set to unity [71,72].n is the normal direction to the interface and is based on the gradient of the phase fraction. More detail can be found in [21,22]. For a cavitating mixture, two additional phase fractions, namely α l for the liquid and α v for the vapor, are calculated with respect to the volume of the cavitating phase and The proposed formulation is equivalent to the classical formulations adopted in threephase system, where all phase fractions are referred to the volume of the computational cell [22]. Combining Equation (4) with Equation (27), it is possible to recover a mixture density for a three-phase system, consistent with Equation (1). Starting from the void fraction of the cavitating mixture, the vapor phase fraction is calculated [20,53,59,73]: where ρ l sat and ρ v sat are the saturation densities of liquid and vapor. Further details on the phase-change model and how to detect ρ lv are explained in the following pages.

Continuum Barotropic Model for Phase Change
A barotropic model for phase change is used here to simulate systems, including both a free surface and low-pressure vaporization. The choice is motivated by the physics of the validation test, for which experimental data were available. The homogeneous equilibrium model (HEM) [53] is applied under the assumptions of thermodynamic equilibrium (a single temperature is used for the mixture) and mechanical equilibrium (the slip velocity is neglected among the phases). The fluid mixture is treated in an homogeneous form, and the thermodynamics is controlled by a barotropic equation of state (EoS): In Equation (33), ψ lv is the mixture compressibility of the homogeneous cavitating mixture; ψ lv is linked to the phase fractions of liquid and vapor through a correlation, that can be either linear [74] or non-linear [20], or tabulated from the properties derived from CoolProp [75] or modeled [76,77]. Because of its simplicity, a linear correlation for the mixture compressibility ψ lv is used: Liquid and vapor compressibility ψ l and ψ v are calculated by the equation of state. For the liquid, a linear EoS is used: where p sat is the saturation pressure and the liquid compressibility is modeled through a relation coming from the Tait equation [78]: In Equation (36), n T is the material parameter and B is the bulk coefficient for the liquid considered, while ρ ∞ and p ∞ are the fluid density and pressure at a specific reference condition. In this work, n T ≈ 7 and B ≈ 300 MPa are used, while the ambient condition is taken as reference. For the vapor, the perfect gas relation is used: where the compressibility of the vapor ψ v is a function of the fluid temperature: being that R is the universal gas constant in J/(kmol K) and T is the fluid temperature. The density of the cavitating fluid is then calculated by the formulation used by Karrholm [73]: The non-condensable gases are assumed to be a perfect gas: Equations (39) and (40) are combined with Equation (4) and the linear compressiblity approximation (see Equation (34)) to give the complete mixture density for the three-phase system: Equation (41) is the continuous barotropic relation for the compressible three-phase mixture and links pressure and density.

Energy Equation
For a compressible VOF method, the solution of the energy equation allows to include the effect of the flow temperature in the compressibility of the gaseous phase; see (38). According to the HEM [53], the phases are considered in thermal equilibrium at the interface. The conservation of the total energy E can be written as [61]: In Equation (42), potential energy is neglected. The total energy E can be expressed as the sum of specific internal energy and the specific kinetic energy (E = e + K); in Equation (42), the thermal contribution is represented by the heat flux, while the mechanical energy includes the stress forces [79,80]. The energy equation is generally implemented in the form of total energy without the mechanical sources [81].
A heat flux q = −α eff ∇e is assumed [61], where the effective thermal diffusivity α eff is the sum of laminar and turbulent thermal diffusivities.
The energy equation written for a single-fluid mixture can be written for temperature as where α eff is the mixture thermal diffusivity, as defined in Equation (7).

Solution Algorithm
The code resolves the governing equations by the finite volume (FV) solution method; a cell-centered formulation with co-located arrangement is used for the sequential solution of the governing equations on a polyhedral mesh. The segregated solution of the governing equations (mass and momentum) is achieved by a pressure-velocity coupling algorithm. The steps performed by the developed segregated solver over a single time-step n are summarized in Figure 2. (25):

Code Verification
Code verification is performed on the test case of Figure 3, in the following named oscillatingWaterPipe: a one-dimensional pipe filled by water in the middle, and by air at its closed ends (see Figure 3). The pressure difference across the pipe triggers the oscillation of the liquid region. The test is aimed to verify the numerical properties of the implemented solver (in the following, referred to as barotropicFoam) in terms of its ability to (a) capture the air/mixture interface, while maintaining their sharpness, and (b) preserve the conservativeness and the boundedness of the solution of the phase-fraction equations with phase change. Finally, the robustness of the solver and its application to the description of real flow physics is demonstrated on a test case available from the literature [82].
Fluid Properties Oscillating water pipe. The geometry and the case setup (initial and boundary conditions) are shown in Figure 3: the two interfaces between liquid and air are initially positioned at x l = 0.1 m and x r = 0.8 m, and the pipe length is 1 meter. Densities and fluid properties were calculated according the EoS defined in Equations (35), (37) and (40). In the liquid region, a linear profile of pressure was used to initialize the field: A similar setup was presented in [83], while in [84][85][86], an initial velocity field was set instead of the pressure field. No phase-change phenomena are involved in this test. In the simulation, the initial pressure field (Figure 4) promotes an oscillation of the liquid core in the pipe; this oscillation evolves in time, thanks to the fluid compressibility.  The aim of this set of simulations is to verify if the captured interface is sharp and if conservation of mass and of the phase fractions is preserved. An analytical solution does not exist for this problem; in accordance with [83], the incompressible equations for one-dimensional incompressible fluid, which are valid in the early time steps, are assumed as the reference for comparison. This is in accordance with what was proposed in [83].
The momentum conservation in one dimension for an incompressible inviscid flow is Being that u = dx dt , it follows that Since no external forces are applied to the fluid flow ( ∂ u ∂x = 0), an Equation of motion (EoM) is obtained and integrated between the two liquid-air interfaces: Under the assumption of incompressible flow, L = x l − x r is constant. Equation (49) links the liquid acceleration to the pressure drop across the interfaces. With varying density, the liquid acceleration is a function of time of the local density and temperature; then, the approximation is assumed to be reliable only in the early time steps of the simulation, when the compressibility effects are still negligible. In the simulations, the relative mass error for the mixture [86] is monitored: where m i (t) are the different masses evaluated during the simulation through an integration on the volume: and Finally, the velocity U x in the middle section of the computational domain is compared against the incompressible solution.
Simulations were carried out on five different equally spaced grids, ranging from 32 to 512 cells. A second order approximation was applied to the time and space derivatives; adaptive time stepping was used for time marching. Figure 5 shows the evolution of α nc during the first oscillation of the liquid column; in Figure 5a-c, the liquid moves toward the right, while in Figure 5d-f, the liquid is bounced back to recover the original position, Figure 5a. Different grid sizes provide different descriptions of the interface; 128 cells is found to be the minimum number of cell to achieve a decent resolution to capture a sharp interface. With coarser grids (32 and 64 cells), the numerical diffusion smears the interface. In Figure 6a, the evolution of the relative error for the five grids is shown. For each of the grids studied, initial conditions (difference of pressure across the liquid) promote the oscillation of the liquid column. When the simulation starts, the fluid at rest starts oscillating; after some time, the amplitude and the oscillation stabilize. The dissipation caused by the flow viscosity dumps the oscillation amplitude. During the early time steps, a smaller time step favors a more accurate solution of the strong, unsteady nature of the problem. For a given CFL max , a finer grid forces a smaller time step, and a lower relative error (see Figure 6a) is observed. At t > 0.025, unsteadiness is determined mostly by the compressibility, which determines the main oscillating frequency of the system. In this case, the time-step marching of any of the tests presented in Figure 6a ensures a proper temporal resolution to capture the physics of the problem, and the error stabilizes. The error is always very limited; for the finest grid (512 cells), it is lower than 0.01%. In Figure 6b, the exact solution for incompressible flows, Equation (49), is used as reference to verify the computed solutions in the early time steps; it is apparent how the discrepancy between the solutions becomes apparent after t = 0.2 s, when the effects of the interface sharpness start having an impact on the solution. With more than 128 cells, the agreement between the computed and the analytical solution is satisfying. Tests on the same geometry were carried out to check the stability of the solver and, most importantly, to verify its accuracy with different maximum time steps allowed. The maximum Courant numbers were tested, namely CFL = 0.1, CFL = 0.01, and CFL = 0.001. The 512 cell grid was used for this set of tests. In Figure 7a, the computed solution is very similar to the analytical. In Figure 7b, the simulation time is grouped in five intervals, and the distribution of the time-step size (∆t) for each interval is reported. Adaptive time step is applied in the simulations, according to the limitation provided by the CFL. A more precise sinusoidal behavior is reached for smaller time steps (Figure 6a); this is related to the semi-explicit discretization used [71] in the solution of the phase fraction equations ( Figure 8).
As expected, mass conservation is also influenced by the max CFL used, as a consequence of the different error in the solution of the phase fraction equations [22]. From Figure 9, it can be observed an overshoot in the error. Simulations with phase change. The liquidColumn test case presented in [22] is used to assess the conservative properties of the solver developed in this work in the presence of phase change. The test case consists of a one-dimensional domain of height L over the ydirection, open at its top. Half of the volume of the column is filled by the non-condensable gas, while the remaining volume contains the cavitating mixture. In the original test [22], the volume fractions of liquid and fuel vapor in the cavitating mixture are 95% and 5% respectively. The system is initially at rest, and a hydrostatic distribution of pressure is set along the y-axis (see Figure 10). No-slip boundary conditions are applied at the lower boundary, while free slip is set on the side walls. The upper boundary of the domain is an open end. The saturation pressure is set to p sat = 103,000 Pa; as pressure is lower than the threshold, the liquid cavitates, and the non-condensable gases exit through the outlet end. The original test case of [22] is then modified: the open upper boundary is converted into a solid wall, and the domain is transformed into a closed vessel. When the initial pressure distribution makes the liquid cavitate, non-condensable gases are compressed, and their density increases ( Figure 11); this second test case necessarily requires a compressible flow solver. In the following, the two tests are named as liquidColumn and cavitatingTank, respectively. The geometry and the boundary conditions of the two test cases presented are summarized in Figure 10. The initial distributions of pressure, α nc and density in the domain are shown in Figure 12, in Figure 12b, it is possible to see a plot distribution of volume fraction α nc referred to as non-condensable gases. The grids used were made respectively of 1 × 640 cells (this grid will be referred in the following as grid A) and 1 × 1280 cells (grid B, in the following). Variable time stepping is used for the simulations presented in this section to preserve a maximum Courant number CFL max = 0.1. Second-order differencing schemes are applied both for temporal and spatial derivatives. Monitored benchmark quantities are the overall mass conservation and the instantaneous mass balance between the liquid fuel and the fuel vapor.

1.
Evolution of mass for each phase: 2. Time evolution of the relative mass error, to verify if mass is conserved during phase change: 3. Global mass relative error: The hydrostatic distribution of pressure used for initialization (see Figure 12) favors the phase change of liquid from the center of the column, until an equilibrium condition is reached. Fuel vapor is generated by phase change where the liquid is present; buoyancy makes the fuel vapor move over the liquid. For the effect of the increasing volume of the mixture, the non-condensable gas is pushed out of the domain through the open boundary. When the equilibrium condition is reached, the fuel is a uniform mixture of vapor and liquid. Results for grids A and B are similar; minor discrepancies are observed in the late times of the simulation (in Figure 13). In Figure 13, the sum of the mass of the vapor and of the liquid are equal to 1-α nc . After phase change, the propagation of the acoustic waves causes a variable pressure field in space and time; for this reason, condensation is triggered in specific locations of the domain. The temporal evolution of the surface front is shown in Figure 14.
The following quantities were monitored during the simulations: • Volume weighted average void fraction: • Global conservation of the volume weighted void fraction: where in Equations (56) and (57), n c is number of cells, while V is the total volume of the mesh: Figure 15 shows that the error is small and oscillates in time, with very limited peaks; as expected, it decreases for finer grids (see Table 1).  cavitatingTank. The set-up of the cavitatingTank test case is described in Figure 10; the domain is the same as the liquidColumn, but the upper boundary is now closed. As the liquid cavitates, the volume of the mixture increases, and the non-condensable gas is compressed. Flow compressibility therefore plays an important role.
In Figure 16, a pseudo-periodic trend in the temporal evolution of the plotted quantities is apparent; it is a consequence of the acoustic pressure waves traveling over the domain that cause a periodic triggering of vaporization and condensation and a periodic compression and expansion of the non-condensable gas. Oscillations are damped by the fluid viscosity. Again, the mass conservation error (Figure 16c is very limited, and it further decreases as the grid is refined. In Figure 17, the phase change of the cavitating mixture is monitored for a time period of 0.5 s. The amount of generated vapor is quite limited: the non-condensable gas is compressed in the closed domain and limits the expansion of the fuel vapor and the reduction in local pressure. This effect is noticed with any of the grids tested. The error in the conservation of the void fraction is still limited (Figure 18).    The global mass relative error, Equation (55) is very limited for both the grids tested (see Table 2)

Validation: High Pressure Liquid Injection
The experimental setup and measurements available from [82] were used for code validation. The experiment consists of observations on several real-size injectors with different converging nozzles, and it has been largely studied in the literature [53,73,80,87] for studies on high-pressure liquid injection. The geometry simulated in this work is named "U nozzle" and has a rectangular cross-section with 5% of contraction. The width is w = 301 µm, and it has an inlet diameter D in = 301 µm and an outlet D out = 284 µm, with an inlet radius of r = 20 µm, while the nozzle length is L = 1 mm. The geometry simulated is reported in Figure 19, while fluid properties are summarized in Table 3.  Different operating conditions were tested; the upstream pressure in the experiment was set at 10 bar, while the downstream pressure varied between 1.5 bar and 8 bar. In this way, different ∆p across the nozzle were tested.  [82] in a single position along the nozzle, at a distance of 57 µm in the nozzle hole; in Figure 19, it is labeled as v 1 .
A polyhedral mesh featuring 5 M cells was used for the simulations. Refinement boxes were applied near the nozzle, where phase change could occur, Figure 20. The wall-adapting local eddy-viscosity model (WALE) [88] was used in the simulations. Timeaveraging on the monitored quantities was performed after 55 µs of simulation for five flow-through times. for studies on high-pressure liquid injection. The geometry simulated in this work is named "U nozzle" and has a rectangular cross-section with 5% of contraction. The width is w = 301 µm, and it has an inlet diameter D in = 301 µm and an outlet D out = 284 µm, with an inlet radius of r = 20 µm, while the nozzle length is L = 1 mm. The geometry simulated is reported in Figure 19, while fluid properties are summarized in Table 3.  Different operating conditions were tested; the upstream pressure in the experiment was set at 10 bar, while the downstream pressure varied between 1.5 bar and 8 bar. In this way, different ∆p across the nozzle were tested.  [82] in a single position along the nozzle, at a distance of 57 µm in the nozzle hole; in Figure 19, it is labeled as v 1 .
A polyhedral mesh featuring 5 M cells was used for the simulations. Refinement boxes were applied near the nozzle, where phase change could occur, Figure 20. The wall-adapting local eddy-viscosity model (WALE) [88] was used in the simulations. Timeaveraging on the monitored quantities was performed after 55 µs of simulation for five flow-through times. In the first operating condition simulated (∆p = 55 bar), no cavitation is present, and the set-up used can be tested without phase change. As documented in [82], phase change is triggered at ∆p = 60 bar. For these conditions, simulation results were compared in Figure 21 against LDV measurements of the flow velocity at a distance of 57 µm in the nozzle hole. The discrepancy observed in the plots can be the combined effect of multiple sources, namely (a) the resolution of the grid, that was limited to cope with the computational resources available for this work; (b) the linear approximation applied to In the first operating condition simulated (∆p = 55 bar), no cavitation is present, and the set-up used can be tested without phase change. As documented in [82], phase change is triggered at ∆p = 60 bar. For these conditions, simulation results were compared in Figure 21 against LDV measurements of the flow velocity at a distance of 57 µm in the nozzle hole. The discrepancy observed in the plots can be the combined effect of multiple sources, namely (a) the resolution of the grid, that was limited to cope with the computational resources available for this work; (b) the linear approximation applied to model the mixture compressibility; and (c) the uncertainties on the inflow conditions and on the LDV measurements. The mass flow rate (MFR) for each operating condition is calculated at the nozzle outlet aṡ m = n ∑ i ρ i (x, t)|u i (x, t) · S f i | (59) where ρ i and U i are the mixture density and velocity at the i-th cell of the outlet and The calculated MFR value is compared against the experiments in Figure 21. A comparison with the experimental results is also provided at incipient cavitation in Figure 22. Finally, the temporal evolution of the vapor fraction ( Figure 23) shows that the vapor clouds form at the nozzle entrance near the detachment region. Early cavitation pockets are symmetric in the spanwise direction; as the cavitation intensity becomes stronger and the vapor clouds evolve, the interaction with the turbulent flow contributes to form non-symmetrical three-dimensional flow structures.

Conclusions
The fully compressible three-phase VOF solver combines the VOF compressible formulation and the HRM model, and it is able to account for thermal effects and phase compressibilities. Different numerical tests were presented for code verification. A first test case, oscillatingWater, was employed to test the compressible interface-capturing method and the conservative properties of the solver. A second test case, cavitatingTank, was used to verify the conservative properties of the solver in a closed domain, where flow compressibility becomes relevant. It is shown that the interface is captured from the phase fraction equation of the non-condensable gas, which does not include any source term for phase change; this aspect favors strong conservation properties of the solver. Finally, code validation is proposed on an experiment proposed in the literature [82]. In the solver, the mixture compressibility of the homogeneous cavitating mixture is linked to the phase fractions of the liquid and vapor through a correlation, which can be linear [74], non-linear [20], tabulated or directly calculated by CoolProp [75], or modeled [76,77]. Although it has been thought to be reserved for aerospace applications, the strategy is applicable to a wide range of problems.