Next Article in Journal
Velocity Distribution in Channels with Submerged Vegetation
Previous Article in Journal
Unsteady Fluid Flows in the Slab Mold Using Anticlogging Nozzles
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Daniel Angel Palomino Solis
1,* and
Federico Piscaglia
2
1
Brembo SGL Carbon Ceramic Brakes S.p.A, 24040 Stezzano, Italy
2
Department of Aerospace Science and Technology (DAER), Politecnico di Milano, 20156 Milano, Italy
*
Author to whom correspondence should be addressed.
Fluids 2022, 7(9), 289; https://doi.org/10.3390/fluids7090289
Submission received: 22 July 2022 / Revised: 21 August 2022 / Accepted: 24 August 2022 / Published: 30 August 2022

Abstract

:
We present a fully compressible single-fluid volume of fluid (VOF) solver with phase change for high-speed flows, 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 fluids (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 significant 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 verification 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 flow, available in the literature, are finally used for a first validation with phase change.

1. Introduction

Increasing the injection pressure in aerospace propulsion systems is a way to improve the combustion efficiency and to reduce CO2 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 (non-zero 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.

1.1. 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.

1.2. 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].

1.3. 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 Section 2, Section 3, Section 4 and Section 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 three-dimensional internal nozzle flow case is reported in Section 8. The main conclusions are drawn in Section 9.

2. 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 ] :
α i = V i V
being the subscript i = lv used for the cavitating mixture and i = nc for the non-condensable gas (nc). The following constraint applies:
α l v + α n c = 1
The single-fluid approximation leads to the so-called mixture variables, namely the following:
-
Mixture velocity:
u = α n c u n c + α l v u l v
-
Mixture density:
ρ = α l v ρ l v + α n c ρ n c
-
Mixture viscosity:
μ = α n c μ n c + α l v μ l v
-
Mixture heat capacity:
c v = α n c c v n c + α l v c v l v
-
Mixture thermal diffusivity:
α eff = α n c α eff , n c + α l v α eff l v
The phase fraction equations for the multiphase compressible cavitating flow are written as a system of two equations:
( ρ n c α n c ) t + · ( u ρ n c α n c ) = 0 ( ρ l v α l v ) t + · ( u ρ l v α l v ) = 0
that has to obey the compatibility condition Equation (2). From the conservation of the phase fractions, it follows that
D α n c D t = α n c ρ n c D ρ n c D t
and
D α l v D t = α l v ρ l v D ρ l v D t
Equations (9) and (10) can be combined to obtain
D α n c D t + D α l v D t = α n c ρ n c D ρ n c D t + α l v ρ l v D ρ l v D t
The LHS of Equation (11) is the total derivative of the compatibility constraint Equation (2):
D α n c D t + D α l v D t = t ( α n c + α l v = 1 ) + · u ( α n c + α l v = 1 ) = · u
so Equation (11) can be written as
· u = α n c ρ n c D ρ n c D t + α l v ρ l v D ρ l v D t
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:
( u ρ ) t + · ( ρ u u ) = p + · τ g · x ρ + f σ
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:
f σ = S ( t ) σ κ n ^ · δ ( x x s ) d S
The term σ is the fluid surface tension in [ N / m ] , and n ^ is the unit vector normal to the liquid interface, whose center is located in x s . κ is the interface curvature [ m 1 ] defined as
κ = · ( n ^ · S f )
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 ^ = α n c | α n c |
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:
f σ = S ( t ) σ κ n ^ · δ ( x x s ) d S σ κ α
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-surface-stress (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.

3. Phase-Fraction Equations

From the combination of
· ( u α i ) = u ( α i ) + α i ( · u )
with Equation (9), an intermediate form of the phase fraction equations can be derived:
α n c t + u ( α n c ) + α n c ( · u ) = · ( u α n c ) = α n c ρ n c D ρ n c D t
If Equations (13) and (20) are combined, it follows:
α n c t + u ( α n c ) = α n c 1 ρ n c D ρ n c D t α n c ρ n c D ρ n c D t α l v ρ l v D ρ l v D t = α n c [ 1 ρ n c D ρ n c D t ( 1 α n c ) = α l v α l v ρ l v D ρ l v D t ] = α n c α l v 1 ρ l v D ρ l v D t 1 ρ n c D ρ n c D t
Adding and subtracting α n c ( · u ) to Equation (21) and recalling Equation (19), the final form of the phase fraction equation in the compressible form is recovered:
α n c t + · ( u α n c ) = α n c α l v 1 ρ l v D ρ l v D t 1 ρ n c D ρ n c D t d g / d t + α n c ( · u )
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
ψ i = D ρ D p
it follows
d g / d t = α n c α l v 1 ρ l v D ρ l v D t 1 ρ n c D ρ n c D t = α n c α l v ψ l v ρ l v ψ n c ρ n c D p D t
The compressibility term and α n c ( · u ) present in Equation (22) are solved in a semi-explicit fashion. The resulting form of the phase fraction equations, is, therefore,
α n c t + · ( u α n c ) + · ( α n c α n c U c ) = α l v α n c 1 ρ l v D ρ l v D t 1 ρ n c D ρ n c D t + α n c ( · u ) α l v = 1 α n c
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 “counter-gradient” 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
U c = min [ C α | u | , max [ | u | ] n ^
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 V l v = V l + V v :
α l = V l V l v α v = V v V l v
and
α l + α v = 1
with
ρ = ( ρ v α v + ρ l α l ) α l v + α n c ρ n c = α l v α v ρ v + α l v α l ρ l + α n c ρ n c
so
α l v α v = V l v V V v V l v = V v V
α l v α l = V l v V V l V l v = V l V
The proposed formulation is equivalent to the classical formulations adopted in three-phase 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]:
α v = ρ l v ρ l s a t ρ l s a t ρ v s a t
where ρ l s a t and ρ v s a t are the saturation densities of liquid and vapor. Further details on the phase-change model and how to detect ρ l v are explained in the following pages.

4. 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):
D ρ l v D t = ψ l v D p D t
In Equation (33), ψ l v is the mixture compressibility of the homogeneous cavitating mixture; ψ l v 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 ψ l v is used:
ψ l v = α l ψ l + α v ψ v
Liquid and vapor compressibility ψ l and ψ v are calculated by the equation of state. For the liquid, a linear EoS is used:
ρ l = ρ l s a t + ( p p s a t ) ψ l = ρ l s a t + ψ l p
where p s a t is the saturation pressure and the liquid compressibility is modeled through a relation coming from the Tait equation [78]:
ψ l = ρ l ρ l / ρ n T n T p + B
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:
ρ v = ψ v p
where the compressibility of the vapor ψ v is a function of the fluid temperature:
ψ v = 1 R T
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]:
ρ l v = ψ l v p + ρ l 0 α l + ( α l ψ l + α v ψ v ψ l v ) p s a t
The non-condensable gases are assumed to be a perfect gas:
ρ n c = ψ n c p
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:
ρ = α n c ψ n c p + α l v ψ l v p + ρ l s a t α l + ( α l ψ l + α v ψ v ψ l v ) p s a t = ( α n c ψ n c + α l v ψ l v ψ ) p + α l v α l ρ l 0
Equation (41) is the continuous barotropic relation for the compressible three-phase mixture and links pressure and density.

5. 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]:
ρ E t + · ( u ρ E ) = · ( u · Σ ) mechanical · q thermal
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].
D ( ρ e ) D t + D ( ρ K ) D t + · ( u p ) = · q
A heat flux q = α eff e is assumed [61], where the effective thermal diffusivity α eff is the sum of laminar and turbulent thermal diffusivities.
D D t ρ ( e + K ) + · ( u p ) = · ( α eff e )
The energy equation written for a single-fluid mixture can be written for temperature as
( ρ T ) t + · ( u ρ T ) + α l v c v l v + α n c c v n c ( ρ K ) t + · ( u ρ K ) + · ( u p ) = α eff 2 T
where α eff is the mixture thermal diffusivity, as defined in Equation (7).

6. 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.

7. 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].
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:
p ( x l x x r ) = p r p l x r x l x + p l x r + p r x l x r x l
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
u t + u u x = 1 ρ l p x
Being that u = d x d t , it follows that
u t + u u x = t d x d t + d x d t u x = 0 = d 2 x d t 2
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:
d 2 x d t 2 = 1 ρ l p x | x l x r = 1 ρ l p r p l L
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:
ε l v ( t ) = m l v ( t ) m l v , 0 ε n c ( t ) = m n c ( t ) m n c , 0 ε t o t ( t ) = m l v ( t ) + m n c ( t ) m n c , 0 + m l v , 0
where m i ( t ) are the different masses evaluated during the simulation through an integration on the volume:
m i ( t ) = 0 t x L x R α i ( t ) ρ i ( t ) d x
and
m i , 0 = m i ( t = 0 )
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 α n c 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 C F L m a x , 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 y-direction, 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 s a t = 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, α n c and density in the domain are shown in Figure 12, in Figure 12b, it is possible to see a plot distribution of volume fraction α n c 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 C F L 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.
  • Evolution of mass for each phase:
    M i = α i ρ i d y
  • Time evolution of the relative mass error, to verify if mass is conserved during phase change:
    E m a r c h i n g = | ( M l ( t + 1 ) + M v ( t + 1 ) ) ( M l ( t ) + M v ( t ) | M l ( t + 1 ) + M v ( t + 1 )
  • Global mass relative error:
    E g l o b a l = | ( M l ( f ) + M v ( f ) ) ( M l ( 0 ) + M v ( 0 ) ) | M l ( 0 ) + M v ( 0 )
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- α n c . 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:
    α i ¯ = j = 1 n c α i j V j V
  • Global conservation of the volume weighted void fraction:
α ¯ = i = 1 3 α i ¯
where in Equations (56) and (57), n c is number of cells, while V is the total volume of the mesh:
V = j = 1 n c V j
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 lcavitatingTank 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)

8. 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 i n = 301 μ m and an outlet D o u t = 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.
A time-averaged velocity profile was available from [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. Time-averaging on the monitored quantities was performed after 55 μs of simulation for five flow-through times.
In the first operating condition simulated ( Δ p = 55 ), 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 as
m ˙ = i n ρ i ( x , t ) | u i ( x , t ) · S f i |
where ρ i and U i are the mixture density and velocity at the i-th cell of the outlet and
S f i = n i · | S f i |
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.

9. 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.   

Author Contributions

These authors equally contributed to this work. Conceptualization, Methodology, Writing—original draft, Editing: D.A.P.S. and F.P.; Software, D.A.P.S.; Resources, Supervision, Writing—review: Federico Piscaglia. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

Authors gratefully acknowledge the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory (Lemont, U.S.) for the computing resources provided.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Vallet, A.; Burluka, A.A.; Borghi, R. Development of a eulerian model for the atomization of a liquid jet. At. Sprays 2001, 11, 619–642. [Google Scholar] [CrossRef]
  2. Mayer, W.; Smith, J. Fundamentals of Supercritical Mixing and Combustion of Cryogenic Propellants. In Liquid Rocket Thrust Chambers; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2012; pp. 339–367. [Google Scholar] [CrossRef]
  3. Smith, J.; Schneider, G.; Suslov, D.; Oschwald, M.; Haidn, O. Steady-state high pressure LOx/H2 rocket engine combustion. Aerosp. Sci. Technol. 2007, 11, 39–47. [Google Scholar] [CrossRef]
  4. Paravan, C.; Galfetti, L.; Bisin, R.; Piscaglia, F. Combustion processes in hybrid rockets. Int. J. Energetic Mater. Chem. Propuls. 2019, 18, 255–286. [Google Scholar] [CrossRef]
  5. Piscaglia, F.; Giussani, F.; Montorfano, A.; Hélie, J.; Aithal, S. A MultiPhase Dynamic-VoF solver to model primary jet atomization and cavitation inside high-pressure fuel injectors in OpenFOAM. Acta Astronaut. 2019, 158, 375–387. [Google Scholar] [CrossRef]
  6. Martinez, J.; Piscaglia, F.; Montorfano, A.; Onorati, A.; Aithal, S.M. Influence of spatial discretization schemes on accuracy of explicit LES: Canonical problems to engine-like geometries. Comput. Fluids 2015, 117, 62–78. [Google Scholar] [CrossRef]
  7. Hélie, J.; Piscaglia, F. On the Influence of Nozzle-generated Coherent Structures on the Primary Atomization of Fuel Jets. In Proceedings of the ICLASS 2021, International Conference on Liquid Atomization and Spray Systems, Edinburgh, Scotland, 29 August–2 September 2021. [Google Scholar]
  8. Piscaglia, F. Developments in Transient Modeling, Moving Mesh, Turbulence and Multiphase Methodologies in OpenFOAM. In Proceedings of the Keynote Lecture at the 4th Annual OpenFOAM User Conference 2016, Cologne, Germany, 11–13 October 2016. [Google Scholar]
  9. Piscaglia, F.; Montorfano, A.; Onorati, A. Development of a non-reflecting boundary condition for multidimensional nonlinear duct acoustic computation. J. Sound Vib. 2013, 332, 922–935. [Google Scholar] [CrossRef]
  10. Cailloux, M.; Helie, J.; Reveillon, J.; Demoulin, F.X. Large Eddy Simulation of a Cavitating Multiphase Flow for Liquid Injection. J. Phys. Conf. Ser. 2015, 656, 012081. [Google Scholar] [CrossRef]
  11. Kolev, N.I. Multiphase Flow Dynamics 3; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  12. Reid, B.; Hargrave, G.; Garner, C.; Wigley, G. An investigation of string cavitation in a true-scale fuel injector flow geometry at high pressure. Phys. Fluids 2010, 22, 031703. [Google Scholar] [CrossRef]
  13. Brennen, C.E. Cavitation and Bubble Dynamics; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]
  14. Schmidt, D.P.; Corradini, M.L. The internal flow of diesel fuel injector nozzles: A review. Int. J. Engine Res. 2001, 2, 1–22. [Google Scholar] [CrossRef]
  15. Ishimoto, J.; Sato, F.; Sato, G. Computational Prediction of the Effect of Microcavitation on an Atomization Mechanism in a Gasoline Injector Nozzle. J. Eng. Gas Turbines Power 2010, 132, 082801. [Google Scholar] [CrossRef]
  16. Ishimoto, J.; Oike, M.; Kamijo, K. Two-Dimensional Numerical Simulation of Boiling Two-Phase Flow of Liquid Nitrogen. Jpn. Soc. Aeronaut. Space Sci. Trans. 2001, 43, 114–121. [Google Scholar] [CrossRef]
  17. Deshpande, M.; Feng, J.; Merkle, C.L. Numerical Modeling of the Thermodynamic Effects of Cavitation. J. Fluids Eng. 1997, 119, 420–427. [Google Scholar] [CrossRef]
  18. Örley, F.; Trummler, T.; Hickel, S.; Mihatsch, M.S.; Schmidt, S.J.; Adams, N.A. Large-eddy simulation of cavitating nozzle flow and primary jet break-up. Phys. Fluids 2015, 27, 086101. [Google Scholar] [CrossRef]
  19. Edelbauer, W. Numerical simulation of cavitating injector flow and liquid spray break-up by combination of Eulerian—Eulerian and Volume-of-Fluid methods. Comput. Fluids 2017, 144, 19–33. [Google Scholar] [CrossRef]
  20. Mithun, M.G.; Koukouvinis, P.; Gavaises, M. Numerical simulation of cavitation and atomization using a fully compressible three-phase model. Phys. Rev. Fluids 2018, 3, 064304. [Google Scholar] [CrossRef]
  21. Piscaglia, F.; Giussani, F.; Hèlie, J.; Lamarque, N.; Aithal, S. Vortex Flow and Cavitation in Liquid Injection: A Comparison between High-Fidelity CFD Simulations and Experimental Visualizations on Transparent Nozzle Replicas. Int. J. Multiph. Flow 2021, 138, 103605. [Google Scholar] [CrossRef]
  22. Giussani, F.; Piscaglia, F.; Saez-Mischlich, G.; Hélie, J. A three-phase VOF solver for the simulation of in-nozzle cavitation effects on liquid atomization. J. Comput. Phys. 2020, 406, 109068. [Google Scholar] [CrossRef]
  23. Ahmed, A.; Duret, B.; Reveillon, J.; Demoulin, F.X. Numerical simulation of cavitation for liquid injection in non-condensable gas. Int. J. Multiph. Flow 2020, 127, 103269. [Google Scholar] [CrossRef]
  24. Trummler, T.; Schmidt, S.J.; Adams, N.A. Investigation of condensation shocks and re-entrant jet dynamics in a cavitating nozzle flow by Large-Eddy Simulation. Int. J. Multiph. Flow 2020, 125, 103215. [Google Scholar] [CrossRef] [Green Version]
  25. Wermelinger, F.; Rasthofer, U.; Hadjidoukas, P.E.; Koumoutsakos, P. Petascale simulations of compressible flows with interfaces. J. Comput. Sci. 2018, 26, 217–225. [Google Scholar] [CrossRef]
  26. Schillaci, E.; Antepara, O.; Balcázar, N.; Rigola Serrano, J.; Assensi, O. A numerical study of liquid atomization regimes by means of conservative level-set simulations. Comput. Fluids 2019, 179, 137–149. [Google Scholar] [CrossRef]
  27. Coutier-Delgosha, O.; Fortes-Patella, R.; Reboud, J.L. Evaluation of the Turbulence Model Influence on the Numerical Simulations of Unsteady Cavitation. J. Fluids Eng. 2003, 125, 38–45. [Google Scholar] [CrossRef]
  28. Huang, B.; Ducoin, A.; Young, Y.L. Physical and numerical investigation of cavitating flows around a pitching hydrofoil. Phys. Fluids 2013, 25, 102109. [Google Scholar] [CrossRef]
  29. Huang, B.; Young, Y.L.; Wang, G.; Shyy, W. Combined Experimental and Computational Investigation of Unsteady Structure of Sheet/Cloud Cavitation. J. Fluids Eng. 2013, 135, 071301. [Google Scholar] [CrossRef]
  30. Reboud, J.L.; Stutz, B.; Coutier, O. Two-phase flow structure of cavitation: Experiments and modelling of unsteady effects. In Proceedings of the 3rd International Symposium on Cavitation, Grenoble, France, 7–10 April 1998. [Google Scholar]
  31. Koukouvinis, P.; Naseri, H.; Gavaises, M. Performance of turbulence and cavitation models in prediction of incipient and developed cavitation. Int. J. Engine Res. 2017, 18, 333–350. [Google Scholar] [CrossRef]
  32. Huang, H.; Sukop, M.; Lu, X.C. Multiphase Lattice Boltzmann Methods: Theory and Application; Wiley: Hoboken, NJ, USA, 2015; pp. 1–373. [Google Scholar] [CrossRef]
  33. Liu, G.; Liu, M. Smoothed Particle Hydrodynamics: A Meshfree Particle Method; World Scientific Publishing Company: Singapore, 2003. [Google Scholar] [CrossRef]
  34. Reynolds, A.J. Thermo-fluid Dynamic Theory of Two-phase Flow. By M. ISHIL. Eyrolles 1975. 248 pp. J. Fluid Mech. 1976, 78, 638–639. [Google Scholar] [CrossRef]
  35. Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y.J. A Front-Tracking Method for the Computations of Multiphase Flow. J. Comput. Phys. 2001, 169, 708–759. [Google Scholar] [CrossRef] [Green Version]
  36. Du, J.; Fix, B.; Glimm, J.; Jia, X.; Li, X.; Li, Y.; Wu, L. A simple package for front tracking. J. Comput. Phys. 2006, 213, 613–628. [Google Scholar] [CrossRef]
  37. Wallis, G.B. Critical two-phase flow. Int. J. Multiph. Flow 1980, 6, 97–112. [Google Scholar] [CrossRef]
  38. Baer, M.; Nunziato, J. A Two-Phase Mixture Theory for the Deflagration-to-Detonation Transition (ddt) in Reactive Granular Materials. Int. J. Multiph. Flow 1986, 12, 861–889. [Google Scholar] [CrossRef]
  39. Saurel, R.; Abgrall, R. A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows. J. Comput. Phys. 1999, 150, 425–467. [Google Scholar] [CrossRef]
  40. Sethian, J.A.; Smereka, P. Level Set Methods for Fluid Interfaces. Annu. Rev. Fluid Mech. 2003, 35, 341–372. [Google Scholar] [CrossRef]
  41. Sussman, M.; Puckett, E.G. A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows. J. Comput. Phys. 2000, 162, 301–337. [Google Scholar] [CrossRef]
  42. Ménard, T.; Tanguy, S.; Berlemont, A. Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet. Int. J. Multiph. Flow 2007, 33, 510–524. [Google Scholar] [CrossRef]
  43. Wang, Z.; Yang, J.; Koo, B.; Stern, F. A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves. Int. J. Multiph. Flow 2009, 35, 227–246. [Google Scholar] [CrossRef]
  44. Hirt, C.W.; Nichols, B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. [Google Scholar] [CrossRef]
  45. Niedźwiedzka, A.; Schnerr, G.H.; Sobieski, W. Review of numerical models of cavitating flows with the use of the homogeneous approach. Arch. Thermodyn. 2016, 37, 71–88. [Google Scholar] [CrossRef]
  46. Shukla, R.K.; Pantano, C.; Freund, J.B. An interface capturing method for the simulation of multi-phase compressible flows. J. Comput. Phys. 2010, 229, 7411–7439. [Google Scholar] [CrossRef]
  47. Miller, G.H.; Puckett, E.G. A High-Order Godunov Method for Multiple Condensed Phases. J. Comput. Phys. 1996, 128, 134–164. [Google Scholar] [CrossRef]
  48. Saurel, R.; Abgrall, R. A Simple Method for Compressible Multifluid Flows. SIAM J. Sci. Comput. 1999, 21, 1115–1145. [Google Scholar] [CrossRef]
  49. Plesset, M.S. Cavitating Flows; Technical Report; California Intstitute of Technology: Pasadena, CA, USA, 1969. [Google Scholar]
  50. Plesset, M.S. The dynamics of cavitation bubbles. J. Appl. Mech. 1949, 16, 277–282. [Google Scholar] [CrossRef]
  51. Plesset, M.S.; Prosperetti, A. Bubble dynamics and cavitation. Annu. Rev. Fluid Mech. 1977, 9, 145–185. [Google Scholar] [CrossRef]
  52. Rayleigh, L. VIII. On the pressure developed in a liquid during the collapse of a spherical cavity. Lond. Edinb. Dublin Philos. Mag. J. Sci. 1917, 34, 94–98. [Google Scholar] [CrossRef]
  53. Schmidt, D.P.; Rutland, C.J.; Corradini, M.L. A fully compressible, two-dimensional model of small, high-speed, cavitating nozzles. At. Sprays 1999, 9, 255–276. [Google Scholar] [CrossRef]
  54. Altimira, M.; Fuchs, L. Numerical investigation of throttle flow under cavitating conditions. Int. J. Multiph. Flow 2015, 75, 124–136. [Google Scholar] [CrossRef]
  55. Gopalan, S.; Katz, J. Flow structure and modeling issues in the closure region of attached cavitation. Phys. Fluids 2000, 12, 895–911. [Google Scholar] [CrossRef]
  56. Adams, N.A.; Schmidt, S.J. Shocks in Cavitating Flows. In Bubble Dynamics and Shock Waves; Delale, C.F., Ed.; Springer: Berlin/Heidelberg, Germany, 2013; pp. 235–256. [Google Scholar] [CrossRef]
  57. Egerer, C.P.; Hickel, S.; Schmidt, S.J.; Adams, N.A. Large-Eddy Simulation of Turbulent Cavitating Flow in a Micro Channel. Phys. Fluids 2014, 26, 085102. [Google Scholar] [CrossRef]
  58. Schnerr, G.H.; Sezal, I.H.; Schmidt, S.J. Numerical investigation of three-dimensional cloud cavitation with special emphasis on collapse induced shock dynamics. Phys. Fluids 2008, 20, 040703. [Google Scholar] [CrossRef]
  59. Rahbarimanesh, S.; Brinkerhoff, J.; Huang, J. Development and Validation of a Homogeneous Flow Model for Simulating Cavitation in Cryogenic Fluids. Appl. Math. Model. 2017, 56, 584–611. [Google Scholar] [CrossRef]
  60. Yu, H.; Goldsworthy, L.; Brandner, P.; Garaniya, V. Development of a compressible multiphase cavitation approach for diesel spray modelling. Appl. Math. Model. 2017, 45, 705–727. [Google Scholar] [CrossRef]
  61. The OpenFOAM Foundation. OpenFOAM User Guide. Available online: https://cfd.direct/openfoam/user-guide (accessed on 24 July 2022).
  62. Issa, R.I. Solution of the implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys. 1986, 62, 40–65. [Google Scholar] [CrossRef]
  63. Ferziger, J.H.; Perić, M. Computational Methods for Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar]
  64. Suhas, V.P. Numerical Heat Transfer and Fluid Flow; Hemisphere Publishing Corporation: London, UK, 1980. [Google Scholar]
  65. de Villiers, E.; Gosman, A.; Weller, H. Large Eddy Simulation of Primary Diesel Spray Atomization. SAE Trans. 2004, 113, 193–206. [Google Scholar] [CrossRef]
  66. Brackbill, J.U.; Kothe, D.B.; Zemach, C. A continuum method for modeling surface tension. J. Comput. Phys. 1992, 100, 335–354. [Google Scholar] [CrossRef]
  67. Popinet, S. An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 2009, 228, 5838–5866. [Google Scholar] [CrossRef]
  68. Gueyffier, D.; Li, J.; Nadim, A.; Scardovelli, R.; Zaleski, S. Volume-of-Fluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional Flows. J. Comput. Phys. 1999, 152, 423–456. [Google Scholar] [CrossRef]
  69. Fedkiw, R.P.; Aslam, T.; Merriman, B.; Osher, S. A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method). J. Comput. Phys. 1999, 152, 457–492. [Google Scholar] [CrossRef]
  70. Weller, H.G. The Development of a New Flame Area Combustion Model Using Conditional Averaging; Thermo-Fluids Section Report TF/9307; Department of Mechanical Engineering, Imperial College of Science Technology and Medicine: London, UK, 1993. [Google Scholar]
  71. Weller, H.G. A new Approach to VOF-Based Interface Capturing Methods for Incompressible and Compressible Flow; Technical Report; TR/HGW/04; OpenCFD Ltd.: Bracknell, UK, 2008. [Google Scholar]
  72. Wardle, K.E.; Weller, H.G. Hybrid Multiphase CFD Solver for Coupled Dispersed/Segregated Flows in Liquid-Liquid Extraction. Int. J. Chem. Eng. 2013, 2013, 128936. [Google Scholar] [CrossRef]
  73. Peng Karrholm, F. Numerical Modelling of Diesel Spray Injection, Turbulence Interaction and Combustion. Ph.D. Thesis, Chalmers University of Technology, Gothenburg, Sweden, 2008. [Google Scholar]
  74. Weller, H.G.; Tabor, G.; Jasak, H.; Fureby, C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 1998, 12, 620–631. [Google Scholar] [CrossRef]
  75. Bell, I.H.; Wronski, J.; Quoilin, S.; Lemort, V. Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp. Ind. Eng. Chem. Res. 2014, 53, 2498–2508. [Google Scholar] [CrossRef]
  76. Nedderman, R.M. One-Dimensional Two-Phase Flow. BY G. B. WALLIS. McGraw Hill, 1969. 408pp. £7. 18s. Cocurrent Gas-Liquid Flow. Edited by E. RHODES AND D. S. SCOTT. Plenum Press, 1969. 698 pp. $27.50. J. Fluid Mech. 1970, 42, 428–430. [Google Scholar] [CrossRef]
  77. Chung, T.J. Computational Fluid Dynamics; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar] [CrossRef]
  78. Li, Y.H. Equation of State of Water and Sea Water. J. Geophys. Res. 1967, 72, 2665–2678. [Google Scholar] [CrossRef]
  79. Jasak, H. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. Ph.D. Thesis, Imperial College London, University of London, London, UK, 1996. [Google Scholar]
  80. Yu, H.; Goldsworthy, L.; Brandner, P.; Li, J.; Garaniya, V. Modelling thermal effects in cavitating high-pressure diesel sprays using an improved compressible multiphase approach. Fuel 2018, 222, 125–145. [Google Scholar] [CrossRef]
  81. Greenshields, C.; Weller, H. Notes on Computational Fluid Dynamics: General Principles; CFD Direct Ltd.: Reading, UK, 2022. [Google Scholar]
  82. Winklhofer, E.; Kull, E.; Kelz, E.; Morozov, A. Comprehensive Hydraulic and Flow Field Documentation in Model Throttle Experiments under Cavitation Conditions. In Proceedings of the ILASS-Europe 2001, 17 International Conference on Liquid Atomization and Spray Systems, Zurich, Switzerland, 2–6 September 2001; Available online: https://www.ilasseurope.org/ICLASS/ilass2001/ILASS2001.pdf (accessed on 24 July 2022).
  83. Duret, B.; Canu, R.; Reveillon, J.; Demoulin, F.X. A pressure based method for vaporizing compressible two-phase flows with interface capturing approach. Int. J. Multiph. Flow 2018, 108, 42–50. [Google Scholar] [CrossRef]
  84. Daru, V.; Le Quéré, P.; Duluc, M.C.; Le Maître, O. A numerical method for the simulation of low Mach number liquid–gas flows. J. Comput. Phys. 2010, 229, 8844–8867. [Google Scholar] [CrossRef]
  85. Kadioglu, S.Y.; Sussman, M.; Osher, S.; Wright, J.P.; Kang, M. A second order primitive preconditioner for solving all speed multi-phase flows. J. Comput. Phys. 2005, 209, 477–503. [Google Scholar] [CrossRef]
  86. Koren, B.; Lewis, M.; van Brummelen, E.; van Leer, B. Riemann-Problem and Level-Set Approaches for Homentropic Two-Fluid Flow Computations. J. Comput. Phys. 2002, 181, 654–674. [Google Scholar] [CrossRef]
  87. Zhao, H.; Quan, S.; Dai, M.; Pomraning, E.; Senecal, P.; Xue, Q.; Battistoni, M.; Som, S. Validation of a Three-Dimensional Internal Nozzle Flow Model Including Automatic Mesh Generation and Cavitation Effects. J. Eng. Gas Turbines Power 2014, 136, 1–11. [Google Scholar] [CrossRef]
  88. Nicoud, F.; Ducros, F. Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor. Flow Turbul. Combust. 1999, 62, 183–200. [Google Scholar] [CrossRef]
Figure 1. Available methods for the solution of two-phase problems. In this work, an algebraic-type VoF is applied to capture the interface.
Figure 1. Available methods for the solution of two-phase problems. In this work, an algebraic-type VoF is applied to capture the interface.
Fluids 07 00289 g001
Figure 2. Solution algorithm employed for a single iteration in a time step in the three-phase single-fluid VOF solver.
Figure 2. Solution algorithm employed for a single iteration in a time step in the three-phase single-fluid VOF solver.
Fluids 07 00289 g002
Figure 3. Oscillating water set-up: geometry and fluid properties.
Figure 3. Oscillating water set-up: geometry and fluid properties.
Fluids 07 00289 g003
Figure 4. Oscillating water set-up, initial conditions: (a) pressure profile; (b) air distribution; (c) fluid properties.
Figure 4. Oscillating water set-up, initial conditions: (a) pressure profile; (b) air distribution; (c) fluid properties.
Fluids 07 00289 g004
Figure 5. Oscillating water pipe. Evolution of void fraction α n c for the first oscillation; time evolution from right to left ((ac) and (df)) They were simulated 5 grids. Legend: () 512 cells; (⋯⋯) 256 cells; (- - -) 128 cells; (Fluids 07 00289 i001) 64 cells; (Fluids 07 00289 i002) 32 cells.
Figure 5. Oscillating water pipe. Evolution of void fraction α n c for the first oscillation; time evolution from right to left ((ac) and (df)) They were simulated 5 grids. Legend: () 512 cells; (⋯⋯) 256 cells; (- - -) 128 cells; (Fluids 07 00289 i001) 64 cells; (Fluids 07 00289 i002) 32 cells.
Fluids 07 00289 g005
Figure 6. Oscillating water pipe: (a) evolution of the relative error ε t o t ; (b) evolution of the velocity magnitude in the core region. Simulations were carried out at CFL max =0.1. Legend: () 512 cells; (⋯⋯) 256 cells; (- - -) 128 cells; (Fluids 07 00289 i001) 64 cells; (Fluids 07 00289 i002) 32 cells.
Figure 6. Oscillating water pipe: (a) evolution of the relative error ε t o t ; (b) evolution of the velocity magnitude in the core region. Simulations were carried out at CFL max =0.1. Legend: () 512 cells; (⋯⋯) 256 cells; (- - -) 128 cells; (Fluids 07 00289 i001) 64 cells; (Fluids 07 00289 i002) 32 cells.
Fluids 07 00289 g006
Figure 7. Oscillating water pipe 512 × 1 cells: (a) evolution of the velocity magnitude; (b) distribution of the time-steps in time. Legend: () approximated solution; (⋯⋯) CFL max = 0.1 ; (- - -) CFL max = 0.01 ; (Fluids 07 00289 i001) CFL max = 0.001 .
Figure 7. Oscillating water pipe 512 × 1 cells: (a) evolution of the velocity magnitude; (b) distribution of the time-steps in time. Legend: () approximated solution; (⋯⋯) CFL max = 0.1 ; (- - -) CFL max = 0.01 ; (Fluids 07 00289 i001) CFL max = 0.001 .
Fluids 07 00289 g007
Figure 8. Oscillating water case: evolution of the left liquid–gas interface for different values of CFL max ; time evolution from right to left ((ac) and (df)). Grid size 512 × 1 cells. Legend: () CFL max = 0.001 ; (⋯⋯) CFL max = 0.01 ; (- - -) CFL max = 0.1 .
Figure 8. Oscillating water case: evolution of the left liquid–gas interface for different values of CFL max ; time evolution from right to left ((ac) and (df)). Grid size 512 × 1 cells. Legend: () CFL max = 0.001 ; (⋯⋯) CFL max = 0.01 ; (- - -) CFL max = 0.1 .
Fluids 07 00289 g008
Figure 9. Oscillating water case. Evolution of the relative error ε t o t for different CFL max . Grid size 512 × 1 cells. Legend: () CFL max = 0.001 ; (⋯⋯) CFL max = 0.01 ; (- - -) CFL max = 0.1 .
Figure 9. Oscillating water case. Evolution of the relative error ε t o t for different CFL max . Grid size 512 × 1 cells. Legend: () CFL max = 0.001 ; (⋯⋯) CFL max = 0.01 ; (- - -) CFL max = 0.1 .
Fluids 07 00289 g009
Figure 10. Test cases for solver verification with phase change, geometry, boundary and initial conditions. Geometry type: (a) open; (b) closed.
Figure 10. Test cases for solver verification with phase change, geometry, boundary and initial conditions. Geometry type: (a) open; (b) closed.
Fluids 07 00289 g010
Figure 11. Evolution of phase masses in time. Legend: () liquidColumn; (- - -) cavitatingTank.
Figure 11. Evolution of phase masses in time. Legend: () liquidColumn; (- - -) cavitatingTank.
Fluids 07 00289 g011
Figure 12. Initial flow conditions; (a) pressure; (b) α n c ; (c) density.
Figure 12. Initial flow conditions; (a) pressure; (b) α n c ; (c) density.
Fluids 07 00289 g012
Figure 13. Temporal evolution of the mass for the i-th phase: (a) Liquid mass, (b) Vapour mass, (c) Air mass. Legend: () grid A (1 × 640 cells); (- - -) grid B (1 × 1280 cells).
Figure 13. Temporal evolution of the mass for the i-th phase: (a) Liquid mass, (b) Vapour mass, (c) Air mass. Legend: () grid A (1 × 640 cells); (- - -) grid B (1 × 1280 cells).
Fluids 07 00289 g013
Figure 14. Temporal evolution of the void fraction profiles along the y-axis for (aj) t = 0.05 s to t = 0.5 s. Grid A: α l , α v , α n c . Grid B: - - - α l , α v , - - - α n c .
Figure 14. Temporal evolution of the void fraction profiles along the y-axis for (aj) t = 0.05 s to t = 0.5 s. Grid A: α l , α v , α n c . Grid B: - - - α l , α v , - - - α n c .
Fluids 07 00289 g014
Figure 15. Evolution in time of (a,d) marching error; (b,e) volume-weighted void fractions; (c,f) sum of the volume-weighted void fractions. Grid A: (ac); grid B (df). Legend: α l (), α v (), α n c ().
Figure 15. Evolution in time of (a,d) marching error; (b,e) volume-weighted void fractions; (c,f) sum of the volume-weighted void fractions. Grid A: (ac); grid B (df). Legend: α l (), α v (), α n c ().
Fluids 07 00289 g015
Figure 16. Temporal evolution of the mass for the i-th phase: (a) Liquid, (b) Vapour, (c) non-condensable gases. Legend: () grid A (1 × 640 cells); (- - -) grid B (1 × 1280 cells).
Figure 16. Temporal evolution of the mass for the i-th phase: (a) Liquid, (b) Vapour, (c) non-condensable gases. Legend: () grid A (1 × 640 cells); (- - -) grid B (1 × 1280 cells).
Fluids 07 00289 g016
Figure 17. Void fraction profiles over the y-axis for (aj) t = 0.05 s to t = 0.5 s. Legend: for the grid A α l (), α v (), α n c (), for the grid B α l (- - -), α v (), α n c (- - -).
Figure 17. Void fraction profiles over the y-axis for (aj) t = 0.05 s to t = 0.5 s. Legend: for the grid A α l (), α v (), α n c (), for the grid B α l (- - -), α v (), α n c (- - -).
Fluids 07 00289 g017
Figure 18. Evolution time of (a,d) marching error; (b,e) volume-weighted void fractions; (c,f) sum of the volume-weighted void fractions. Grid A: (ac), Grid B: (df). Legend: α l (), α v (), α n c ().
Figure 18. Evolution time of (a,d) marching error; (b,e) volume-weighted void fractions; (c,f) sum of the volume-weighted void fractions. Grid A: (ac), Grid B: (df). Legend: α l (), α v (), α n c ().
Fluids 07 00289 g018
Figure 19. Reference geometry for the test case simulated.
Figure 19. Reference geometry for the test case simulated.
Fluids 07 00289 g019
Figure 20. Polyhedral mesh with 5 M cells used for the simulations. Refinement regions (labelled as II and III) are applied in proximity to the nozzle region.
Figure 20. Polyhedral mesh with 5 M cells used for the simulations. Refinement regions (labelled as II and III) are applied in proximity to the nozzle region.
Fluids 07 00289 g020
Figure 21. Velocity profiles at position 1 of Figure 19 for: (a) Δp = 55 bar; (b) Δp = 60 bar. (c) Mass flow rate (MFR) at different operating conditions. Legend: () simulation; (■) experiments.
Figure 21. Velocity profiles at position 1 of Figure 19 for: (a) Δp = 55 bar; (b) Δp = 60 bar. (c) Mass flow rate (MFR) at different operating conditions. Legend: () simulation; (■) experiments.
Fluids 07 00289 g021
Figure 22. Code validation against experimental results: (a) numerical α v prediction; (b) vapor measurement.
Figure 22. Code validation against experimental results: (a) numerical α v prediction; (b) vapor measurement.
Fluids 07 00289 g022
Figure 23. Temporal evolution ( t ( 30 , 90 ] μs) of the vapor fraction at the early instants of the simulation (incipient cavitation, Δ p = 60 bar). Iso-contour lines are plotted for α v = 0.5 .
Figure 23. Temporal evolution ( t ( 30 , 90 ] μs) of the vapor fraction at the early instants of the simulation (incipient cavitation, Δ p = 60 bar). Iso-contour lines are plotted for α v = 0.5 .
Fluids 07 00289 g023
Table 1. Global error.
Table 1. Global error.
No. Cells6401280
E g l o b a l % 0.0440.011
Table 2. Relative error of mass in the cavitatingTank test case for the two grids tested.
Table 2. Relative error of mass in the cavitatingTank test case for the two grids tested.
No. Cells6401280
E g l o b a l % 0.0450.040
Table 3. Fluid properties.
Table 3. Fluid properties.
Fluid Properties
ParameterValue
ρ 998.16 kg
p 2000 Pa
B3.07 × 109 Pa
n T 1.75
μ l 20.5 × 10−3 kg/ms
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Palomino Solis, D.A.; Piscaglia, F. Toward the Simulation of Flashing Cryogenic Liquids by a Fully Compressible Volume of Fluid Solver. Fluids 2022, 7, 289. https://doi.org/10.3390/fluids7090289

AMA Style

Palomino Solis DA, Piscaglia F. Toward the Simulation of Flashing Cryogenic Liquids by a Fully Compressible Volume of Fluid Solver. Fluids. 2022; 7(9):289. https://doi.org/10.3390/fluids7090289

Chicago/Turabian Style

Palomino Solis, Daniel Angel, and Federico Piscaglia. 2022. "Toward the Simulation of Flashing Cryogenic Liquids by a Fully Compressible Volume of Fluid Solver" Fluids 7, no. 9: 289. https://doi.org/10.3390/fluids7090289

APA Style

Palomino Solis, D. A., & Piscaglia, F. (2022). Toward the Simulation of Flashing Cryogenic Liquids by a Fully Compressible Volume of Fluid Solver. Fluids, 7(9), 289. https://doi.org/10.3390/fluids7090289

Article Metrics

Back to TopTop