Benoît Péden
email: peden@cerfacs.fr

Julien Carmona

Pierre Boivin

Schmitt Thomas

Bénédicte Cuenot

Nicolas Odier

Thomas Schmitt

Numerical assessment of Diffuse-Interface method for air-assisted liquid sheet simulation

Keywords: Two-phase flow, airblast, diffuse interface, hyperbolic systems

  L'archive ouverte pluridisciplinaire

Introduction

Liquid injection and atomization are complex processes which are at the heart of many industrial systems. In the context of combustion, these phenomena play a key role as they drive the fuel vapour distribution in a combustion chamber [START_REF] Réveillon | Analysis of weakly turbulent dilute-spray flames and spray combustion regimes[END_REF][START_REF] Luo | Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor[END_REF], and therefore have a major impact on the flame regime and burner efficiency, as well as pollutant emissions [START_REF] Beck | Identification of droplet burning modes in lean, partially prevaporized swirl-stabilized spray flames[END_REF]. In addition, fuels are blends of numerous molecules, their exact composition being usually unknown. For example Shastry et al. [START_REF] Shastry | Numerical study of multicomponent spray flame propagation[END_REF] highlighted the preferential evaporation effect on the structure and propagation speed of spray flames. For these complicated problems numerical simulation represents a convincing approach to supplement experiments. However, atomization is numerically challenging due to the large range of physical scales involved in this multi-physics problem, from the surrounding gas flow to the molecular scale at the interface between liquid and gas phases, considered as a discontinuity for viscosity and density, and consequently requiring specific numerical treatment.

For this reason, a strategy often used in complex industrial systems is to directly inject a train with prescribed droplet velocities and diameter size distribution [START_REF] Senoner | Eulerian and Lagrangian Large-Eddy Simulations of an evaporating two-phase flow[END_REF][START_REF] Sanjosé | Fuel injection model for Euler-Euler and Euler-Lagrange Large-Eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF]. Such spray injection models have allowed successful spray flame simulations in both academic and realistic combustion chambers [START_REF] Shum-Kivan | Experimental and numerical analysis of a turbulent spray flame structure[END_REF][START_REF] Tillou | Evaluation and Validation of Two-Phase Flow Numerical Simulations Applied to an Aeronautical Injector Using a Lagrangian Approach[END_REF][START_REF] Schiavo | On the impact of fuel injection angle in Euler-Lagrange large eddy simulations of swirling spray flames exhibiting thermoacoustic instabilities[END_REF][START_REF] Shastry | Large Eddy Simulations of complex multicomponent swirling spray flames in a realistic gas turbine combustor[END_REF].

However, they require input parameters which are not always known and they suffer from a lack of physics of atomization, which prevent the spray formation to adapt to the flow.

Several approaches have been proposed in the literature to numerically predict the primary atomization of liquid jets or sheets by modelling the interface between liquid and gaseous phases, which can be split into two main categories: sharp and diffuse interface methods. Sharp interface methods consider the interface between the phases as a discontinuity in a way that flow variables experience a jump in their values. The material interface is tracked by the Front-Tracking methods [START_REF] Glimm | Front tracking for hyperbolic systems[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Tryggvason | A Front-Tracking Method for the Computations of Multiphase Flow[END_REF] or captured using methods such as the Volume-Of-Fluid (VOF) method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF], or the Level-Set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. These approaches are able to finely predict the evolution of the interface topology and its primary breakup [START_REF] Pairetti | On shear layer atomization within closed channels: Numerical simulations of a cough-replicating experiment[END_REF][START_REF] Carmona | A comparative study of direct numerical simulation and experimental results on a prefilming airblast atomization configuration[END_REF][18][START_REF] Janodet | A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids[END_REF]. However they require significant mesh resolution in the interface vicinity [START_REF] Pairetti | Mesh resolution effects on primary atomization simulations[END_REF] which makes them unsuitable for prediction in complex industrial systems. In addition they are usually restricted to incompressible flows, i.e., not directly applicable to combustion. Diffuse interface methods (DIM) do alleviate these issues, with considering the interface as a transition region between the phases instead of a discontinuity, thus relaxing the mesh resolution constraints. To do so, they rely on unified equations of state describing both phases. Two sub-groups of DIM arise from the literature: the Phase-Field (PF) and the Multi-Fluid (MF) methods.

PF methods are based on the early works of Van der Waals [START_REF] Van Der Waals | Over de Continuiteit van den Gas-en Vloeistoftoestand[END_REF], Korteweg [START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais continues et sur la théorie de la capillarité dans l'hypothèse d'une variation continue de la densité[END_REF] and Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] to describe the Gibbs potential in the diffuse interface region. The concept is to track the phases with a dedicated field function and to solve each phase independently. Numerous Lattice-Boltzmann numerical formulations rely on this approach, allowing interesting applications for atomization in a incompressible context [START_REF] He | A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability[END_REF][START_REF] Zheng | Lattice Boltzmann interface capturing method for incompressible flows[END_REF][START_REF] Wang | A mass-conserved diffuse interface method and its application for incompressible multiphase flows with large density ratio[END_REF][START_REF] Liang | Axisymmetric lattice Boltzmann model for multiphase flows with large density ratio[END_REF]. For compressible finitevolume or finite-element approaches, the second gradient theory [START_REF] Jamet | The Second Gradient Method for the Direct Numerical Simulation of Liquid-Vapor Flows with Phase Change[END_REF] recently allowed to handle high density ratio in academic configurations [START_REF] Nayigizente | Development of an interface thickening method for the direct numerical simulation of compressible liquid-vapor flows in the framework of the second gradient theory[END_REF]. MF methods rather consider the interface region as a mixture of phases, obeying a single set of equations. Several equilibrium hypotheses have been proposed to reduce the number of equations, from the general 7-equation formulation of Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-todetonation transition (ddt) in reactive granular materials[END_REF], to a reduced 5-equation formulation assuming pressure equilibrium [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Allaire | A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF][START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF], and down to a so-called homogeneous 4-equation formulation, assuming the mechanical and thermal equilibrium between phases [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. Because they involve a only set of conservative equations to describe both phases, without tracking the interface, DIM relax the resolution constraint of sharp methods.

Convex thermodynamic models [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] are used to close the equation set, such as the Stiffened-Gas [START_REF] Métayer | Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF] or cubic equations of state, as recently proposed in several works [START_REF] Pelletier | A multifluid Taylor-Galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states[END_REF]. Multi-Fluid methods allow the simulation of compressible, high density ratio industrial configurations at a reasonable computational cost.

They were used in liquid rocket engine applications [START_REF] Le Touze | A compressible twophase flow framework for Large Eddy Simulations of liquid-propellant rocket engines[END_REF] using cubic [START_REF] Pelletier | A multifluid Taylor-Galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states[END_REF][START_REF] Gaillard | A diffuse interface Lox/hydrogen transcritical flame model[END_REF][START_REF] Gioud | Injection modelling in LOX/GCH4 rocket engines with a diffuse interface method[END_REF][START_REF] Tudisco | Numerical investigations of thermodynamic states formed near a reacting, supercritical CH4/O2 flame[END_REF] or Stiffened-Gas equations of state [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF][START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF][START_REF] Chiapolino | Extended Noble-Abel Stiffened-Gas Equation of State for Sub-and-Supercritical liquid-gas systems far from the critical point[END_REF][START_REF] Deng | Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet[END_REF], or in Diesel injection setups with interesting results [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of high-pressure fuel injection and application to ECN Spray A[END_REF][START_REF] Ma | Large-eddy simulations of transcritical injection and auto-ignition using diffuse-interface method and finite-rate chemistry[END_REF][START_REF] Yang | Real-fluid injection modeling and LES simulation of the ECN Spray A injector using a fully compressible two-phase flow approach[END_REF]. Contrary to PF approach, existing MF works have mainly neglected surface tension forces, which play a key role in atomization [START_REF] Lefebvre | Atomization and Sprays[END_REF]. Primary and secondary atomization and directly the droplets distribution in the combustion chamber are a consequence of the surface tension forces.

The present work intends to assess the accuracy and robustness of the homogeneous model (four-equation model [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]) to predict liquid-gas flow dynamics. A Noble-Abel Stiffened-Gas (NASG) equation of state is used to close the system and the Continuum Surface Force formalism [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] extended to a compressible formulation [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF] is considered to account for surface tension forces. As recommended in [START_REF] Saurel | Diffuse-Interface Capturing Methods for Compressible Two-Phase Flows[END_REF], the vertex-centered HLLC-MUSCL Riemann solver [START_REF] Cheng | AVIP: a low temperature plasma code[END_REF] is used for a robust and accurate prediction of the high density ratio faced in the targeted applications. Finally, characteristic boundary conditions which are required for compressible flows, are adapted to NASG thermodynamics. The model is first validated on a series of academic test cases, and finally applied to liquid film dynamics in an airblast atomizer.

The paper is structured as follows. The transport equations and closure models are presented in Sec. 2, then the numerical scheme is described in Sec. 3. Validations for the characteristic boundary conditions, and the surface tension implementation are proposed in Sec. 4. Finally, the results for the target application are application are presented in Sec. 5. Average and instantaneous results are compared to existing correlations, and the effect of pressure injection is investigated.

Model description

The 4-equation model

Initially proposed to simulate evaporating, cavitating and boiling flows [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF][START_REF] Lund | A Hierarchy of Relaxation Models for Two-Phase Flow[END_REF][START_REF] Lemartelot | Steady one-dimensional nozzle flow solutions of liquid-gas mixtures[END_REF][START_REF] Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF], the homogeneous 4-equation model assumes a mechanical and thermal equilibrium between the phases. Mechanical and thermal equilibriums are always true in pure phase media. In a multi-phase mixture, they mean that no slip velocity exists between phases, both phases have the same temperature. This implies that the region where both phases coexist is small. Species composing the mixture occupy their own volume and share the same pressure whether they are in liquid or gas phase which is similar to Dalton's law [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF]. With this hypotheses, the 4-equation model considers a mono-component liquid (k = 1) and a multi-component gas mixture (k = [2, . . . , N ]), distinguishing the liquid species and its vapour as two different species. The equilibrium leads to the following variable definitions

               T = T k , ∀k, P = P k , ∀k, v = 1/ρ = N k=1 Y k v k , e = N k=1 Y k e k . (1)
The resulting 4 equations write

∂ρ ∂t + ∇•(ρu) = 0, (2a)
∂ρu ∂t + ∇•(ρu ⊗ u + P I) = 0, (2b)
∂ρE ∂t + ∇•((ρE + P )u) = 0, (2c)
∂ρY k ∂t + ∇•(ρY k u) = 0, (2d)
where Y k represents the mass fraction of species k, ρ, u, P and E respectively denote the mixture density, the velocity, the pressure and the total mixture energy (E = e + u 2 /2). The internal energy writes e. As the 4-equation model governed by Sys. (2) does not contain any non-conservative term, it is fully conservative.

Thermodynamic closure

A consistent Equation of State (EoS) must be used to correctly relate density, temperature and pressure of liquid and gas phases. This EoS should be able to reproduce the thermodynamic behaviour of species on a large range of temperatures and pressures.

As suggested by Le Métayer et al. [START_REF] Métayer | Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF], in the context of the 4-equation model using a convex EoS for each phase allows to preserve the hyperbolicity of the Sys. (2) which simplifies its resolution. The liquid and its vapour are thus transported through two distinct equations [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] and mass-transfer between phases occurs assuming an instantaneous local thermodynamic equilibrium. When the conditions are far from the saturation, no phase change occurs and the local thermodynamic equilibrium can be neglected. This is the case in this work, i.e., no evaporation/condensation are considered.

The Noble-Abel Stiffened-Gas equation of state

Cubic equations of state [START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF][START_REF] Peng | A New Two-Constant Equation of State[END_REF] are good candidates to describe both phases but they loose convexity in the spinodal region of the Clapeyron diagram. The Noble-Abel Stiffened-Gas (NASG) equation of state [START_REF] Métayer | The Noble-Abel Stiffened-Gas equation of state[END_REF] extends the popular Stiffened-Gas formulation which considers covolumes and expresses as for a given constituent k = [1, . . . , N ] as

P k = ρ k (γ k -1)C v,k T k 1 -ρ k b k -P ∞,k , (3)
where γ k , C v,k , b k and P ∞,k are respectively the Laplace coefficient, the isovolume specific heat capacity, the covolume and the constant representing attractive effects in liquid phase: the gaseous phase relaxes to the ideal gas EoS with P ∞,k = 0 and b k = 0.

Internal energy, enthalpy, entropy and the speed of sound read

e k = P k + γ k P ∞,k P k + P ∞,k C v,k T k + e 0 k , (4)
h k = γ k C v,k T k + b k P k + q k , (5)
s k = C v,k ln T γ k k (P k + P ∞,k ) γ k -1 + s 0 k , (6)
c 2 k = γ k (P k + P ∞,k ) ρ k (1 -ρ k b k ) , (7)
where the two added constants e 0 k and s 0 k respectively represent the reference energy and entropy of species k. The NASG EoS parameters are determined according the methodology proposed in Boivin et al. [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF]. In order to avoid small time-steps induced by the CFL constraint in the liquid phase, these parameters are derived assuming a modified sound speed in the liquid phase, equal to the one of the gaseous phase [START_REF] Carmona | Modélisation des phénomènes diphasiques dans des injecteurs aéronautiques de type Airblast[END_REF].

Extension of NASG equation of state to multi-component gas mixtures

Extension of Sys. (1) to multi components gas has been proposed in Chiapolino et al. [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], and is reminded here for clarity. Writing Eq. ( 3) for each species with ρ k = 1/v k , the mixture temperature is easily retrieved

T (v, P, Y k ) = v - N k=1 Y k b k N k=1 Y k (γ k -1)C v,k P +P ∞,k . (8)
The mixture pressure is then obtained by rewriting Eq. ( 8), yielding a quadratic form

P = b + √ b 2 + 4ac 2a , (9)
with

         a = v -Y 1 b 1 , b = N k=1 Y k (γ k -1)C v,k T -P ∞,1 (v -Y 1 b 1 ), c = N k>1 Y k (γ k -1)C v,k P ∞,1 T. (10)
Note that a proof of the pressure positivity was presented by Chiapolino et al. [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], and that the formula is only valid if a single liquid component is present (here, k = 1).

Finally density and internal energy are computed through Eqs. ( 3) and ( 4) according to Sys. ( 1)

1 ρ = N k=1 Y k (γ k -1)C v,k T P + P ∞,k + b k , (11) e
= N k=1 Y k P + γ k P ∞,k P + γ k P ∞,k C v,k T + e 0 k . (12)

Jacobian flux matrix & boundary conditions

Characteristic boundary conditions [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] are required in a compressible context to be consistent with the system of waves travelling inside the computational domain. The Locally One-Dimensional Inviscid flow (LODI) hypothesis introduced by Thompson [START_REF] Thompson | Time dependent boundary conditions for hyperbolic systems[END_REF] is considered. Following Okong'o and Bellan [START_REF] Okong'o | Consistent Boundary Conditions for Multicomponent Real Gas Mixtures Based on Characteristic Waves[END_REF], characteristic formulations must be adapted when using non-perfect equation of state. Pelletier et al. [START_REF] Pelletier | A multifluid Taylor-Galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states[END_REF] have derived the characteristic formulation for a 4-equation model, closed with a cubic equation of state. A similar exercise is conducted here, for the NASG equation of state. The computation of coefficients Λ and Γ involved in the Jacobian fluxes matrix as well as the speed of sound are detailed in Appendix A,

Λ = ∂P ∂ρe ρY k = 1 ρ Cpβ α -αT , (13)
Γ = ∂ρe ∂ρY k P,Y i̸ =k = h k - ρ Cp αρ k , (14)
with α and β respectively the isobaric dilatation and isotherm compressibility coefficients. The mean isobaric heat capacity is defined as

Cp = k Y k C p k .

Additional terms: surface tension & viscous term 145

Surface tension is accounted for following the methodology proposed by Brackbill et al. [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], taking the liquid volume fraction α l as the continuous color function. The resulting force F σ yields

F σ = σκ∇α l , (15)
with σ the surface tension coefficient and κ the local curvature defined as

κ = -∇• ∇α l |∇α l | . ( 16
)
The viscous stress tensor and the energy flux read

τ = - 2 3 ∇•(u)I + 2S, (17)
q = -λ∇T, (18)
with S = 1 2 (∇u + ∇u t ). The mixture dynamic viscosity µ and the thermal diffusion coefficient λ are defined as in Deng and Boivin [START_REF] Deng | Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet[END_REF]. Using this formulation, coefficients continuity is conserved across the interface,

µ = α l µ l + α v µ v , (19)
λ = α l λ l + α v λ v , (20)
where µ l and λ l are the dynamic viscosity coefficient and the thermal diffusion coefficient of the liquid phase and µ g and λ g are the ones of vapour phase. The volume fractions are defined as α k = ρY k /ρ k . In this study since only one liquid species and a multicomponent gas mixture are considered, α l = α 1 and 19) and ( 20), µ l and λ l are assumed constant and independent on temperature. The gaseous viscosity µ g depends on temperature and its behaviour is modelled with Sutherland's law [START_REF] Sutherland | The viscosity of gases and molecular force[END_REF]. The thermal diffusion coefficient of gaseous phase is computed through the Prandtl number P r, µ g and the mean gaseous mass specific calorific capacity Cp,g ,

α g = N k=2 α k . In Eqs. (
λ g = µ g Cp,g P r , (21)
with

Cp,g = N k=2 Y k C p,k N k=2 Y k . (22)
The complete system of equations with non-Euler terms, taking the surface tension force following [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF], reads

∂ρ ∂t + ∇•(ρu) = 0, (23a)
∂ρu ∂t + ∇•(ρu ⊗ u + P I) = F σ + ∇•(µτ ), (23b)
∂ρE ∂t + ∇•((ρE + P )u) = u • F σ + ∇•(µτ • u) -∇• q, ( 23c
)
∂ρY k ∂t + ∇•(ρY k u) = 0. (23d)

Numerical scheme

A vertex-centered Riemann solver scheme has been developed at CERFACS [START_REF] Cheng | AVIP: a low temperature plasma code[END_REF], and added in the compressible, massively parallel and unstructured Large-Eddy-Simulation (LES) solver AVBP [START_REF] Schönfeld | Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP[END_REF][START_REF] Gicquel | High performance parallel computing of flows in complex geometries[END_REF], to numerically handle with robustness high density gradients. Initially intended for plasma applications, it is adapted here to solve the Sys. (2), in conjunction with the NASG equation of state. Appendix B presents the numerical solver HLLC [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] along with the second order MUSCL reconstruction. The Sweeby [START_REF] Sweby | High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws[END_REF] slope limiter is used for the MUSCL procedure. The temporal integration is performed with a three-step Runge-Kutta (RK) algorithm, as in [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF]. The non-Euler terms from Sys. [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] are computed with specific methodologies. In the capillary term the divergence is computed with a 4∆ methodology presented in [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF][START_REF] Segui-Troth | Multi-physics coupled simulations of gas turbines[END_REF]. The viscous terms are computed with a Finite Element 2∆ operator [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF][START_REF] Segui-Troth | Multi-physics coupled simulations of gas turbines[END_REF] to dissipate small wavelength diagonal modes that could appear. The validation of the Eulerian part of Sys. ( 2) is presented in appendix Appendix B.3.

Canonical validations

In the following, we successively validate:

• the surface tension implementation,

• boundary conditions. The theoretical pressure drop across the interface predicted by the Laplace law is

P l -P v = ∆P = σ r = 14.6 Pa. (24)
Azimuthally averaged pressure drop measured in our simulation is shown in Fig. 1, and leads to a 0.7 % relative error to the theoretical value.

Oscillating droplet test case

The oscillation of an elliptic viscous droplet is then considered, as shown in Fig. The initial species profiles are given by

Y O 2 (x) = Y O 2 ,min + (Y O 2 ,max -Y O 2 ,min ) Φ(x), (25a)
Y N 2 (x) = Y N 2 ,min + (Y N 2 ,max -Y N 2 ,min ) Φ(x), (25b)
Y H 2 O (x) = 1 -Y N 2 (x) -Y O 2 , (25c)
with

Φ(x, y) = 1 2 1 + tanh ϕ(x, y) 2ϵ , (25d)
ϕ(x, y) = 1 - x 2 a 2 - y 2 b 2 , ( 25e
)
where 2a and 2b are respectively the ellipse width and height, and ϵ is the interface width. The evolution of kinetic energy in the domain is shown in Fig. 3, evidencing its periodic evolution, damped by viscous forces. The normalized time t ⋆ is defined as

Coefficients H 2 O (liquid) O 2 N 2 C p,k [J.kg -1 .K -1 ] 4185 919 1045 C v,k [J.kg -1 .K -1 ] 4180 650 742 P ∞,k [Pa] 1.026×10 7 0 0 q k [J.kg -1 ] -7.8998×10 5 -2.104×10 3 -2.355×10 3 b k [m 3 .kg -1 ] 9.2×10 -4 0 0 2a 2b L = 10 mm
t ⋆ = t ρ l r 3 eq σ . ( 26
)
Theoretical studies have been proposed to predict the expected frequency of 195 a 2D oscillating droplet. The work by Rayleigh [START_REF] Rayleigh | On the Capillary Phenomena of Jets[END_REF], extended by Fyfe et al. [START_REF] Fyfe | Surface tension and viscosity with lagrangian hydrodynamics on a triangular mesh[END_REF] to account for surface tension predicts a pulsation

ω 2 = (n 3 -n) σ (ρ l + ρ g )r 3 eq , ( 27
)
where n = 2 is the most unstable mode. The frequency reads

f = ω 2π , (28)
where r eq is the equilibrium droplet radius. The numerical frequency is 73.03 Hz while the theoretical one computed with Eq. ( 28) is 71.61 Hz, resulting in 2 % 200 relative error and falls in the range of previously reported simulations [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF][START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF].

Boundary conditions validation

The Navier-Stokes Characteristic Boundary Condition (NSCBC) formalism for this two-phase framework is first assessed on the evacuation of an acoustic wave, superimposed to an

H 2 O-O 2 mixture at rest, with Y H 2 O l = 1-Y O 2 = 0.98,
which results in a liquid volume fraction α l = 7.55 %. Thermodynamic conditions are P 0 = 1 atm and T 0 = 300 K, and no artificial stabilization is applied.

Gaussian forward and backward acoustic waves are superimposed. The left and right boundary conditions are outlet non-reflecting characteristic conditions [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF].

Both HLLC-MUSCL [START_REF] Cheng | AVIP: a low temperature plasma code[END_REF] and a Two-step Taylor Galerkin (TTGC) scheme [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] are applied with CFL number 0.7. The use of the TTGC scheme is another validation as its formulation includes the Jacobian matrix and thus the thermodynamic specific terms, see Appendix A. Results are shown in Fig. 4, and evidence the fact that characteristic boundary conditions are able to properly handle acoustics without any reflection in this two-phase context.

The evacuation of a purely gaseous bubble from a pure liquid domain is also considered. This corresponds to the exit of an entropic wave. The test case consists in a pure N 2 -O 2 bubble, initialized in a pure H 2 O l domain of 1 m long, discretized by 100 mesh nodes, with an initial velocity u 0 = 10 m.s -1 . The thermodynamic conditions are P 0 = 1 atm, T 0 = 300 K. The initial profile of species across the interface follows Eq. ( 25) but the Φ function is modified and reads

Φ(x) = 1 2 1 + tanh √ 6(x + 0.3) 4∆x - 1 2 1 + tanh √ 6(x + 0.1) 4∆x , ( 29
)
where ∆x is the cell size.

As shown in Fig. 5, the considered characteristic boundary condition allows to correctly evacuate the gaseous bubble from the liquid domain without reflecting waves.

Water droplet advection to non-reflecting boundary condition

This section focuses on the convection of a two-dimensional droplet and aims to investigate the interaction between a curved interface and the outlet The mesh consists in a 1000 × 1000 grid cells.

The results presented in Fig. 6 demonstrate that the droplet is evacuated without producing any reflecting waves.

Application: 2D Airblast atomizer

The intent of this section is to assess the implemented methodology on a practical configuration of interest for aircraft engines. Airblast atomization consists in a low-velocity liquid jet sheared by a faster gaseous co-flow. Such 240 a configuration leads to the development of liquid ligaments which break into droplets [START_REF] Marmottant | Fragmentation of stretched liquid ligaments[END_REF]. It is of particular interest since liquid atomization is a key parameter for combustion quality in aeronautical or automotive applications. For this study, a configuration experimentally studied at the LEGI lab [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF][START_REF] Ben Rayana | Assisted atomization of a liquid layer: investigation of the parameters affecting the mean drop size prediction[END_REF][START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF][START_REF] Matas | Experimental and analytical study of the shear instability of a gas-liquid mixing layer[END_REF][START_REF] Jerome | Vortices catapult droplets in atomization[END_REF][START_REF] Matas | Influence of Gas Turbulence on the Instability of an Air-Water Mixing Layer[END_REF] is considered, where a liquid film injected at low velocity along a horizontal wall is sheared by 245 a high-speed air stream.

Experimental and theoretical works [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF][START_REF] Ben Rayana | Assisted atomization of a liquid layer: investigation of the parameters affecting the mean drop size prediction[END_REF][START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF][START_REF] Matas | Experimental and analytical study of the shear instability of a gas-liquid mixing layer[END_REF] have shown that a primary instability develops due to a Kelvin-Helmholtz instability between the liquid and the gas, while a transverse Rayleigh-Taylor instability leads to the development of liquid ligaments. A linear stability analysis allows to predict the dominant axial wavelength λ axi for the primary instability [START_REF] Marmottant | Fragmentation of stretched liquid ligaments[END_REF][START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF][START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF][START_REF] Matas | Experimental and analytical study of the shear instability of a gas-liquid mixing layer[END_REF] as

λ axi = C axi ρ l ρ g δ g , (30)
where C axi is a constant, ρ l and ρ g are respectively liquid and gas densities.

The gas vorticity thickness δ g is computed using the experimental correlation established in [START_REF] Matas | Experimental and analytical study of the shear instability of a gas-liquid mixing layer[END_REF]

δ g = 6H g Re g , ( 31
)
where Re g is the gaseous flow Reynolds number defined as

Re g = ρ g U g H g µ g . ( 32
)
The velocity of the axial instability can be computed thanks to the Dimotakis velocity U c [START_REF] Dimotakis | Two-dimensional shear-layer entrainment[END_REF] which has shown to be in very good agreement with measured wave velocities [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF][START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF]]

U c = √ ρ g U g + √ ρ l U l √ ρ g + √ ρ l , (33)
and finally the dominant axial frequency is

f axi = U c λ axi . ( 34
)
Several numerical studies have investigated this experimental configuration, relying on incompressible VOF methods in 2D [START_REF] Fuster | Simulation of primary atomization with an octree adaptive mesh refinement and VOF method[END_REF][START_REF] Fuster | Instability regimes in the primary breakup region of planar coflowing sheets[END_REF][START_REF] Agbaglah | Numerical simulation of the initial destabilization of an air-blasted liquid layer[END_REF][START_REF] Desjardins | Direct Numerical and Large-Eddy Simulations of primary atomization in complex geometries[END_REF][START_REF] Bozonnet | Stability of an air-water mixing layer: focus on the confinement effect[END_REF] and in 3D [START_REF] Chiodi | A Numerical Parametric Study on the Air-Blast Atomization of a Planar Liquid Layer[END_REF]. Transversal instabilities are not taken into account in 2D simulations, but results showed a good quantitative agreement with experiments.

The computational domain considered in the present study consists in a 20 H l × 10 H l rectangular box, sketched in Fig. 7a. The numerical domain is discretized with a cell size equal to δ g /n (with n = 2, 4 and 8) near the injection and shear zone (13 H l × 6 H l ), and the mesh is stretched in both vertical and axial directions outside of these regions. The geometric parameters considered in our numerical simulation are the same as in the numerical work of Fuster et al.

[86], i.e., liquid and gas enter in the domain from left with heights respectively equal to H l = 6.4 mm and H g = 10 mm, as shown in Fig. 7a. The considered species are liquid water H 2 O l for liquid flow and respectively 22 % and 78 % in mass of oxygen O v 2 and nitrogen N v 2 for the gaseous flow, and the surface tension coefficient is σ = 69 × 10 -3 N.m -1 . Following Bozonnet et al. [START_REF] Bozonnet | Stability of an air-water mixing layer: focus on the confinement effect[END_REF], no splitter separating plate is considered and continuous profiles are imposed at the inlet, as illustrated in Fig. 7b. The inlet velocity profile is given by Eqs. [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] to [START_REF] Métayer | Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF]:

u(y) = U l erf H l -y δ l + U i 1 + erf y -H l δ l × erf y δ l , (35)
for y ≤ H l .

u(y) = U g erf y -H l δ g + U i 1 -erf y -H l δ l × erf H g + H l -y δ g + U sm 1 + erf y -(H g + H l ) δ g , (36)
for

H l < y ≤ H g + H l . u(y) = U sm 1 -erf y -(H g + H l ) δ g + U cf erf y -(H g + H l ) δ g , (37)
for y > H g + H l , where δ g is the gas vorticity thickness, δ l the liquid vorticity thickness, U i = δ l (U g µ g /δ g + U l µ l /δ l )/(µ l + µ g ) is the interfacial velocity defined in [START_REF] Otto | Viscous instability of a sheared liquid-gas interface: Dependence on fluid properties and basic velocity profile[END_REF], U cf = 0.2 U g is the co-flow velocity and U sm = (U cf +U g )/2 is the gaseous phase mixing layer velocity. Liquid and gaseous dynamic viscosities are respectively denoted µ l and µ g .

At the initial state the domain is at rest with temperature T = 300 K. The air flow is fixed at U g = 20 m.s -1 and liquid velocity is computed to reach a dynamic pressure ratio M = (ρ g U 2 g )/(ρ l U 2 l ) = 16 resulting in U l = 0.1735 m.s -1 at P = 1 bar. The simulation parameters are given in Tab. 2, and species properties are the ones considered in the validation test-cases Tab. 1, for which the liquid water species has a modified sound speed, resulting in a slightly modified density (ρ l = 939 kg.m -3 at 1 bar). Imposed inlet velocity values are kept constant in the study. An air co-flow is added on top of the high-speed air Phase

U [m.s -1 ] µ [Pa.s] σ [N.m -1 ] Liquid 0.1735 1×10 -3 69×10 -3
Gas 20 1.8×10 -5 stream similarly as [START_REF] Chiodi | A Numerical Parametric Study on the Air-Blast Atomization of a Planar Liquid Layer[END_REF], to limit gaseous recirculations in this region. Inlet and outlet boundary conditions are characteristic non-reflecting boundary conditions, validated above in Sec. 4.2 for the current framework.

An instantaneous density field is shown in Fig. 8a, evidencing the development of liquid ligaments. In order to quantitatively assess the results, the liquid core length predicted by the correlation of Raynal [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF] is reported with a thick orange line. This correlation provides the intact liquid core length L thanks to Eq. ( 38),

L = 12 H l √ M , (38)
and is retrieved from previous numerical works relying on sharp interface methods [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF][START_REF] Agbaglah | Numerical simulation of the initial destabilization of an air-blasted liquid layer[END_REF][START_REF] Desjardins | Direct Numerical and Large-Eddy Simulations of primary atomization in complex geometries[END_REF]. Mean profiles of liquid volume fraction α l = 0.1, 0.5 and 0.9 are given in Fig. 8b for the three mesh resolutions (n = 2, 4 and 8). The intact liquid core more likely corresponds to the iso-contour α l = 0.9, and Fig. 8b shows that increasing mesh resolution allows to obtain a fair agreement with the experimental correlation of Raynal Eq. [START_REF] Pelletier | A multifluid Taylor-Galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states[END_REF]. Figure 8c shows the α l = 0.5 isocontour evolution for the resolution n = 8 during 400 ms, corresponding to the passage of approximately 10 waves. The liquid sheet oscillates above the experimental correlation of Raynal [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF] and the film height decreases following the slope of Eq. [START_REF] Pelletier | A multifluid Taylor-Galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states[END_REF] given by the thick black line. The accurate prediction of this liquid core length shows the ability of the diffuse interface method to predict the liquid entrainment within the two-phase mixing layer.

The normalized amplitude spectrograms (A * ) of the interface height for P = 1 bar with different mesh resolutions are given in Fig. 9a (∆x = δ g /n with n given on Fig. 9a), at position x = 2 H l . These spectra have been obtained

with the Welch method [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF] to extract the most amplified frequency. The signal

] (a) 0 1 2 3 0 1 2 x/H l [-] y/H l [-] n = 2 0 1 2 3 x/H l [-] n = 4 0 1 2 3 x/H l [-] n = 8 Raynal α l = 0.1 α l = 0.5 α l = 0.9 (b)
x/H l [-] is divided into ten segments with a 60 % overlap as in Bozonnet et al. [START_REF] Bozonnet | Stability of an air-water mixing layer: focus on the confinement effect[END_REF], and the liquid volume fraction is integrated across vertical lines to determine the liquid height. This process is repeated for multiple axial positions ranging from 0 to 13 times the liquid injection height H l . Afterwards, the signals undergo

y/H l [-] 0 0.5 1 1.5 2
Fast Fourier Transform (FFT), and the results are averaged. Despite an increase of numerical noise the main frequency remains the same for the three meshes, meaning that independence from mesh resolution is reached. The signal is processed for a time t so that the number of generated waves is larger than 30 (τ = f ∆t). Bozonnet et al. [START_REF] Bozonnet | Stability of an air-water mixing layer: focus on the confinement effect[END_REF] showed that even τ ≈ 10 is enough to obtain a good estimate of the frequency. The value found is f = 24 Hz which is equivalent to a Strouhal number St = f δ g /U g = 6.5 × 10 -4 . The value is compared in Tab. 3 with previous studies at the same operating point and falls in the range of reported Strouhal numbers meaning that the essential instability mechanisms are captured by the method.

The resolution n = 2 is adopted for the rest of the study to keep reasonable computational costs. Further simulations are finally conducted to investigate new operating ranges, hardly reachable with both experiments and incompressible solvers. Simulations are performed with various density ratios by adapting the pressure imposed at the outlet. The densities for each pressure are given in Tab. 4. 30) and [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

The gas vorticity thickness δ g is adapted to the operating points in Eqs. [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] and [START_REF] Métayer | Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF] according to Eq. [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF]. Simulation parameters in Tab. 2 are kept constant. Results for the main frequency at various pressures are given in Fig. 9b and are compared to the one obtained from correlation Eq. ( 34) (black solid line)

where a prefactor C axi = 3.2 is considered. The same methodology as in Fig. 9a 340 has been used, however for P = 7 and 8 bar raw FFT are considered because of a gas recirculation at the end of the liquid cone length. This phenomenon has also been experimentally observed by Raynal [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF], Ben Rayana [START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF] for large momentum ratios, typically M > 100. Note that several values for C axi can be found in literature. Raynal [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF] proposed C axi = 4, Ben Rayana [START_REF] Rayana | Contribution à l'étude des instabilités interfaciales liquidegaz en atomisation assistée et tailles de gouttes[END_REF] proposed 345 C axi = 1.5 to 1.8. It has been shown that the different methodologies considered to experimentally measure the frequency can be a reason for these different values of C axi [START_REF] Matas | Experimental and analytical study of the shear instability of a gas-liquid mixing layer[END_REF]. Despite this uncertainty, linear stability analysis [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF] predicts a linear relation between f and U c ρ g /ρ l /δ g , that is perfectly found in our simulations as pressure increases (Fig. 9b). It is to the authors' knowledge a rare confirmation of this theoretical prediction on a practical case.

Normalized amplitude spectrograms over a distance x/H l = 8 are given in Fig. 10. For each operating pressure point the main axial frequency is clearly visible near the injection. The frequency decreases downstream of a length approximately equal to the intact liquid cone length given by Eq. ( 38)

highlighted with a solid black line in Fig. 10.

This finally shows that the sheared liquid-gas interface dynamics is well reproduced by the present method. The intact liquid core length is correctly retrieved, and the frequency evolution with pressure theoretically predicted by linear stability analysis is confirmed by our numerical approach. This methodology can therefore be considered with confidence for pressurized airblast configurations in future studies.

Conclusions

This paper presents first a validation of implemented methods for handling boundary conditions in conjunction with the NASG EoS. The non-reflectivity is assessed for both 1D and 2D two-phase test-cases. Capillary terms associated with the liquid-gas interface have been validated on simple configurations, with satisfactory results.

The methodology is then assessed on an air-assisted atomization configuration, for which experimental, theoretical, and numerical results can be found in the literature. Our results show the ability of the diffuse interface method to predict the intact liquid core length, as well as the interface dynamics. Furthermore, the methodology allows to investigate new operating ranges, and the effect of pressure increase on the primary interfacial wave frequency is studied. Our results show a very good agreement with the theoretical prediction of Raynal [START_REF]Instabilité et entraînement à l'interface d'une couche de mélange liquide-gaz[END_REF]. This allows to bridge the gap between usually incompressible methods considered for primary atomization, and compressible methods considered for prediction of combustion. Future work will be devoted to the development of a phase-change methodology, along with temperature-variable heat capacities, to properly handle reactive conditions. In addition, a coupling with a disperse spray Lagrangian model will be investigated, to allow a proper spray prediction in a realistic pressurized and reactive environment.

Jacobian matrix reads

A x (U) = ∂F(U) ∂U =               0 1 0 0 • • • 0 -u 2 u(2 -Λ) Λ -Λ(Γ 1 -e c ) • • • -Λ(Γ N -e c )
-uH H -Λu where the L 2 -error is defined as

795 L 2 (ε ρ ) = i V i ( ρi -ρ i ) 2 i V i , (

4. 1 .

 1 Surface tension validation 4.1.1. Laplace law Surface tension, taken into account thanks to Eq. (15), is validated in this subsection thanks to several test cases. The Laplace law is first assessed, considering a H 2 O liquid droplet of radius 5 mm, in an atmosphere of a gaseous O 2 -N 2 mixture. The domain is a square of 50 mm side. The mesh consists in a 500 × 500 grid cells, and the surface tension coefficient is σ = 73 × 10 -3 N.m -1 .

2 .Figure 1 :

 21 Figure 1: Pressure difference profiles for simulation (red dashed line) and Laplace (black solid line).

Figure 2 :Figure 3 :

 23 Figure 2: Oscillating droplet test case initialization. Dimensions 2a and 2b are respectively equal to 1.7 mm and 1 mm.

Figure 4 :

 4 Figure 4: Superimposed backward and forward acoustic waves with fully non-reflecting boundary conditions. Red solid line: HLLC-MUSCL RK3 scheme [53], blue dashed line: TTGC scheme [76].

Figure 5 :

 5 Figure 5: Air bubble in a liquid H 2 O environment transport to the outlet boundary condition.

Figure 6 :

 6 Figure 6: 2D liquid H 2 O droplet advection to the outlet boundary condition with. P = P -P 0 .The red circle is α l = 0.5 isocontour. Note that the colorbar is rescaled for the last instant.

  3 x δ g /n ∆x = 1.3 y δ g /n (a) Sketch and mesh discretization of the configuration.

  Inlet velocity profile. The blue dashed curve represents the gaseous flow while the solid red curve represents the liquid flow.

Figure 7 :

 7 Figure 7: 2D Airblast configuration.

Figure 8 :

 8 Figure 8: (a) Instantaneous density field. The thick orange line represents the intact cone length of Raynal [78]. (b) Liquid volume fraction mean profiles (α l = 0.1, 0.5 and 0.9) for the three resolutions used (n = 2, 4 and 8) compared to the correlation of Raynal Eq. (38) (black solid line). (c) Comparison of the liquid cone length at P = 1 bar. The thick black line represents the the correlation length of Eq. (38). The blue area represents 40 filled isocontour of α l = 0.5 over 400 ms with n = 8 resolution.

U

  Normalized amplitudes A * for P = 1 bar. c ρ g /ρ l /δ g [Hz] Main frequency as a function of the pressure and the quantity Uc ρg/ρ l /δg.

Figure 9 :

 9 Figure 9: Results of the spectral analysis. (a) Normalized amplitude A * of the Fast Fourier Transform (FFT) at P = 1 bar for different mesh resolutions at position x = 2 H l and τ ≥ 30. (b) Main frequency evolution with pressure. The solid line represents the frequencies obtained with correlations Eqs. (30) and[START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

Figure 10 :

 10 Figure 10: Amplitude normalized spectrograms (A * ) of interface height for P = 3 to 8 bar over a distance x/H l = 8. The solid black line represents the intact cone length from Eq. (38).

Figure B. 11 :

 11 Figure B.11: Discretization at all the cell nodes. The grey area represents the nodal volume V i . The hatched area is the intersection of a cell τ containing node i with the nodal volume V i .The edges of V i (f ) are defined from the center of cell τ and the middle of its edges with node i (j and h). The virtual nodes k and l are used for the MUSCL procedure.

Figure B. 14 :

 14 Figure B.14: Comparison of the schemes after one convective time with 256 × 256 resolution.

B. 10 )

 10 with ρ the exact solution and V i the node volumes, is shown in Fig. B.15 at CFL = 0.7 for HLLC-MUSCL. A second-order accuracy is recovered, with the slope limiter parameter β having an influence only for coarse meshes.

Table 1 :

 1 NASG thermodynamic coefficients for H 2 O (liquid), O 2 , N 2 species. Determined with NIST reference[START_REF] Linstrom | The NIST Chemistry WebBook: A Chemical Data Resource on the Internet[END_REF] with Boivin et al.[START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF] methodology.

Table 2 :

 2 Simulation parameters: U is the velocity, µ is the dynamic viscosity and σ is the surface tension coefficient.

Table 3 :

 3 Comparison of previous experimental and numerical values of Strouhal numbers with the present study for the case M = 16 and P = 1 bar.

	Strouhal (×10 -4 )

Table 4 :

 4 Liquid and vapour densities for each pressure and the associated momentum ratio M .

	Pressure [bar]	3	4	5	6	7	8
	Liquid density ρ l [kg.m -3 ]	941.6 942.8 943.9 945.1 946.2 947.3
	Vapour density ρ g [kg.m -3 ]	3.4	4.5	5.7	6.8	7.9	9.1
	Momentum ratio M [-]	48	64	80	95	111	127

Acknowledgments

BP is respectively grateful to Jérôme Dabas and Thibault Gioud for fruitful discussions about signal processing and Multi-Fluid methods. This work was supported by the French National Research Agency through the project MADIRAN (grant ANR-21-CE50-0043) and benefited from GENCI-CINES HPC resources (grant 2022 -A0112B10157).

Appendix A. Characteristic boundary conditions

Well applied boundary conditions are of major interest in compressible flows in order to avoid for example pressure drift and acoustic waves in the computational domain. A proper derivation of the coefficients from the Jacobian matrix is essential as is the speed of sound.

One can recall the one-dimensional 4-equation model So Eq. (A.5) becomes

Knowing that d(ρe) = ρde + edρ,

Finally one can write

. Isobaric dilatation and isotherm compressibility coefficients

The Isobaric dilatation α and isothermal compressibility β coefficients re-

In a multicomponent NASG framework the mass-specific volume reads

so that the mixture mass-specific volume derivative at constant composition reads

The mass volume is related to density through So the density derivative reads

Finally one can write

The hyperbolic sound speed reads

The differential of the mixture mass-specific volume at constant composition can be written as so that

The mixture mass specific entropy reads

where Ch. 16).

Appendix B.2. MUSCL reconstruction procedure

A HLLC Riemann solver [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] is considered in this work with a MUSCL reconstruction of variables to compute fluxes at edges. This procedure uses piecewise linear functions to ensure second order accuracy. Cell boundary values are computed as

respectively on the left and on the right sides of node i along a given direction. For stability

i is computed with the limited approach of Sweby [START_REF] Sweby | High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws[END_REF]:

The ∆U n i-1/2 and ∆U n i+1/2 values correspond respectively to ∆ ki and ∆ ij for the node i and ∆ ij and ∆ jl for the node j as shown in

Coefficients

H 2 O (liquid) Air

Results are shown in Fig. B.12, where plots of pressure, density, velocity and liquid volume fraction are depicted at t = 220 µs, and compared to the exact solution. Agreement is seen to be excellent. Consistently with the literature [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF], no overshoot is obtained, and Hugoniot relationships are satisfied.

Appendix B.3.2. Two dimensional Riemann test case

Two dimensional Riemann ideal gas test-cases are then considered. The geometry is a [-0.5, -0.5] × [0.5, 0.5] meter square with a 600 × 600 regular grid.

Specific calorific capacities ratio is constant and equal to γ = 1.4. The initial condition consists in four quadrants with different initial conditions [START_REF] Lax | Solution of Two-Dimensional Riemann Problems of Gas Dynamics by Positive Schemes[END_REF]:

Order of accuracy

The order of accuracy of the scheme is next assessed thanks to the convection of an isentropic Gaussian vortex, defined with the following functions

while the temperature writes

where r = x 2 + y 2 . The pressure finally results from