François Balavoine
email: francois.balavoine@univ-artois.fr

Raphael Romary

Cassoret Bertrand

Cristian Demian

Balavoine François

Romary Raphael

Demian Cristian

Evolution impact of electrical insulation systems in high power electrical machines

I. INTRODUCTION

It is in the 18th century, that the discovery of electrical insulation is assigned by the historians to Stephen Gray. It was with the discovery of electrical conduction and the first experiments that insulating materials are studied. Gray's experiments allow separation of conductive and insulating materials, for example, glass and silk as insulating materials and brass and iron as conducting wire [START_REF] Chipman | An unpublished letter of stephen gray on electrical experiments, 1707-1708[END_REF]. It is from 1884 that research will progress, with the creation of the AIEE (American Institute of Electrical Engineers) for the main aim of advancing "the theory and the practice of electrical engineering and of the allied arts and sciences and the maintenance of a high professional standing among its members." [START_REF] Hickernell | The institute and electrical insulation[END_REF]. There is therefore the creation of research teams to study electrical insulation systems [START_REF] Berberich | Role of aiee in electrical insulation[END_REF]. Nowadays, in electrical machines and typically, the electrical insulation systems are also used to prevent short circuits between conductors and the magnetic circuit. Another purpose is to ensure transmission of the heat released by the conductors to the parts of the magnetic circuit, also the impregnation makes it possible to give a rigidity to the winding, to avoid vibrations despite the magnetic forces that it can be subjected to [START_REF] Stone | Electrical insulation for rotating machines: design, evaluation, aging, testing, and repair[END_REF]. In electrical machines, two types of winding can be identified. The first one, for low power machines, less than 100kW, supplied below 600V is named random winding, mainly made of round wire. The second one, for high-power machines and supplied more than 600V, usually high-voltage, consists of a flat conductor with a generally flat copper plate, is called form-wound coils. The machine windings generally consist of three types of insulation:

• Turn-to-turn insulation, nearest to the conductor, to avoid short circuits between the wires.

• Groundwall insulation, between the conductors and grounded stator or rotor core for separates this, and avoid short-circuits between the winding and the magnetic circuit; • Winding impregnation, used to fill the remaining air spaces between the coils, to solidify it and to give rigidity to the coil heads for mechanical strength [4] [5]. The improvement of insulation is characterised by several things: the thickness, the voltage resistance and the temperature performance, the last point is presented in this article for the turn-to-turn insulation for random winding.

II. TEMPERATURE PERFORMANCES OF THE TURN-TO-TURN INSULATION

A. Old insulation

In the early 1900s, natural fibre materials are mainly used for winding insulation. The most frequently used method of insulation is cotton wrapping. However, other natural fibre materials are used such as linen, silk, wool and later asbestos. This cotton wrapping can be single, but it is usually lined or tripled to prevent some faults during the application, and the important potential difference between the conductors [START_REF] Gray | Electrical machine design[END_REF]. At the beginning of the 1900s, it was observed that the main short-circuits between conductors were caused by the mechanical strength of the insulation (vibrations for example), weakened by the rise in temperature inside the machine. In 1913, the lifespan of insulation was not defined. Steinmetz conducts research to define the temperature performances of some insulating materials [START_REF] Steinmetz | Temperature and electrical insulation[END_REF].

The figure 1 shows the first insulation classes divided into three groups:

• A : Natural fibre materials such as cotton, paper and natural resins. The materials groups A and B are those mainly used for winding, and C for the construction of heaters and rheostats.

For cotton insulation only (column A), the maximum heating temperature can reach 90°C before it deteriorates. The addition of mica or asbestos increases the limit of use to 110°C.

B. Recent insulation

From the 1950s onwards, a considerable amount of research began to replace natural materials by plastic films and nonwoven polymer fibres.

The discoveries made on insulation, made it possible to reduce their thickness, to increase their dielectric qualities, and their temperature resistance.

The new insulation systems require a different technique of deposition on the conductor, enamelling is used. The enamelling allows a homogeneous and calibrates removal of the insulation without default.

The figure 2 shows the most frequently used recent insulation classes for windings, with class H up to 180°C. The differences between the classes are explained by the chemical composition, for example, class F (150°C) is for low-density polyethylene insulation and class H (180°C) for high-density polyethylene [START_REF] Pinel | Fonction isolation dans les matériels électriques[END_REF].

For specific applications, there are materials with a very high temperature resistance, such as polyamide-imide (PAI) with a maximal temperature of 220°C, or polyimide at 240°C. To withstand temperatures of 450°C-500°C, it is possible to use another technical system, with a thin layer coating of ceramic insulation for the conductor and a high temperature cement to consolidate the brittle ceramic insulation layer [9] [10]. But this insulation is fragile and therefore difficult to wind without damaging it.

C. Impact of electrical insulation

The electrical insulation evolution allows the increase of the voltage supply, the dielectric strength and the acceptable temperature. This concerns the high-voltage insulation and the fill factor improvement of copper in the slots [START_REF] Stone | Electrical insulation for rotating machines: design, evaluation, aging, testing, and repair[END_REF] [6] [START_REF] Liwschitz | Electric Machinery, volume II of ACMachines[END_REF]. Tables I andII give an overview of the improvement of electrical insulation.

It can be seen in table I that before 1960s the groundwall insulation was essentially made of micanite tube which is an assembly of thin flakes of mica stuck together with a flexible varnish, the tube thickness depends on voltage supply. The conductor insulation was composed of various papers and cotton. A few turns around the conductors were necessary.

The new electrical insulation presented in the table II, is composed of several types of thin materials, each serves to protect the IM (turn-to-turn insulation, corona effect, mechanical protection). To the nearest of the conductor, the turn-to-turn protection is a base of mica paper with epoxy resin impregnation.

III. INFLUENCE ON MACHINE DESIGN

A. General evolution

The design of electrical machines depends on several factors. The internal dimensions of diameter and length are guided by the speed, namely the number of poles and the power. The dimensions of the stator section are determined by two main factors: firstly, the value of the induction in the air gap, which depends on the supply voltage. The induction in the various parts of the magnetic circuit, together with the quality of the magnetic steels, define the minimum lengths of the teeth and yoke; secondly, the permissible current density as a function of cooling, which determines the surfaces of the slots. So, these two main components of the machine, together with improvements in materials, have made it possible to reduce the size of machines for equivalent power [START_REF] Balavoine | Influence of modern magnetic and insulation materials on dimensions and losses of large induction machines[END_REF].

B. Impact on slot dimensions

Improving electrical insulation systems therefore makes it possible to improve current density because, as can be seen by comparing the figures 1 and 2, there is a possible gain of 90°C if we compare Steinmetz's 110°C class and the current H class 180°C. Therefore, more current can flow through conductors of equivalent size. Then, from the point of view of thicknesses, the tables I and II present two examples of insulation used for high-power asynchronous machines, first in the 1950s for the table I and nowadays for the table II. The gains in thickness are considerable: the thickness of the turn-to-turn insulation is reduced by 50% and for the insulation of the coils combined with that of the bottom of the slot, the thickness is 5 to 6 times less.

C. Refurbishment perspectives

The various gains described in the previous section do not only concern the sizing of new machines. The renovation of electrical machines can also benefit from this development. A large number of old-generation electric motors from the 1950s are still in operation with their original windings. The prospect of renovating the winding with a new winding and new insulation offers several advantages. Firstly, by reducing the thickness of the insulation, the size of the conductors can be increased, and therefore the filling coefficient. With more copper, operating at the same load point reduces Joule losses and heating. Secondly, it is also possible to accept more heat thanks to a better insulation class. This makes it possible either to increase the current to get a little more power, or to improve the longevity of the machine with a stronger winding. This aspect of renovation also has the advantage of having a reduced environmental impact, thanks to reduced Joule losses and higher efficiency at the same operating point.

IV. CONCLUSION

The evolution of electrical insulation systems has been allowed to increase their performances, that is the thickness, the voltage resistance and the temperature performances. This paper presents the evolution of the thermal performance of insulation. The first insulation often based on natural fibres, did not exceed 90°C alone. Nowadays, thanks to the technique of enamelling and research on polymers, it is possible to reach 180°C or even 240°C for more specific insulation. For the design of electrical machines, this improvement is beneficial, it allows the increase of the current density of the winding and therefore reduce its size. There is also a better mechanical strength of the insulation due to the application technical. The better performances of insulation allow them to increase the size of the conductors and consequently to reduce their electrical resistance and improve the efficiency. Another possibility with better insulation is to reduce the size of the machines, which decrease the environmental impacts of the manufacture.

Fig. 1

 1 Fig. 1 Insulation class defined by Steinmetz in 1913

  • B : Asbestos and mica-based materials. • C : Fireproof or heat-resistant mica based materials used in the construction of heaters and rheostats.

Fig. 2

 2 Fig. 2 Recent insulation class for maximum winding temperature

TABLE II

 II

	Recent Electrical Insulation


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

Julien Carmona

Pierre Boivin

Thomas Schmitt

Bénédicte Cuenot

Nicolas Odier

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

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

This study presents the implementation and validation of a second-order accurate solver for the 4-equation multi-fluid method in a cell-vertex context, to handle aeronautical air-assisted liquid sheet configurations. Validations include Laplace tests, droplets oscillations, and a two-dimensional configuration reminiscent of an aeronautical airblast injector. Promising results are obtained in the last case, especially when the pressure is increased, as relevant for aeronautical applications. Being fully compatible of the reactive flow formulation of the cell-vertex solver AVBP, this study paves the way to future monolithic simulations of airblast injectors for aeronautical combustion chambers, including both the multi-phase and the reactive regions.

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 5 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][START_REF] Mukundan | Detailed numerical simulations of primary atomization of airblasted liquid sheet[END_REF][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. A series of pressurized tests are conducted on the last application. It is to the authors knowledge the first time that such airblast applications are conducted with a diffuse interface method, and that pressure effects on the liquid-gas interface dynamics can be investigated under a compressible, eventually reactive formalism.

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] Le Martelot | Steady one-dimensional nozzle flow solutions of liquid-gas mixtures[END_REF][START_REF] Le 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. This model is classified under the category of Multi-Fluid methods, which distinguishes from Phase-Field methods by its approach to interface thickness control. Unlike Phase-Field methods, the simulation does not actively manipulate the interface thickness. This deliberate choice is motivated by the specific applications targeted by the model, where the flow dynamics are sufficiently fast to prevent excessive diffusion of the interface, and by the fact that Sys. (2) corresponds to the usual set of equations solved in combustion solvers [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF].

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 , γ k , C v,k , b k and P ∞,k are respectively the density, 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 + q 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 + q ′ k , (6)
c 2 k = γ k (P k + P ∞,k ) ρ k (1 -ρ k b k ) , (7)
where the two added constants q k and q ′ 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 + P ∞,k C v,k T + q 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

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
)
This formulation of the surface tension force enables a fast and computationally efficient integration of capillary effects into the system of equations. In addition to extending the method to compressible cases, Perigaud and Saurel [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF] formulated the equation system in conservative form, offering improved accuracy in capturing capillary phenomena. However, the implementation of this formulation poses certain challenges, as outlined below:

(i) Modification of the energy equation. A surface tension potential energy E σ is added to the conservative energy variable E. This potential energy term is defined as

ρE σ = σ∥∇α l ∥. (17)
The energy equation then reads,

∂(E + σ∥∇α l ∥) ∂t + ∇• (E + σ∥∇α l ∥ + P ) u - σ ∥∇α l ∥ I - ∇α l ⊗ ∇α l ∥∇α l ∥ . (18)
(ii) Impact on flux formulation and solver: the conservative momentum and energy formulations lead to modifications in the flux formulation, requiring adjustments in the Riemann solver and the Jacobian flux matrix.

Considering that the AVBP solver [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] is primarily used for single-phase gas applications, where the focus is mainly on simplicity and flexibility of implementation, the non-conservative form of the equations is preferred. This choice enables easier integration of the code with the desired level of simplicity and flexibility, aligning with the solver's primary application requirements. The surface tension effects are incorporated into the numerical simulations using a non-conservative formulation. This approach, despite its departure from a fully conservative framework, yields satisfactory results in validation test cases. The implementation of the non-conservative formulation within the AVBP solver has proven to be straightforward and efficient. By adopting this approach, the AVBP solver, which primarily focuses on single-phase gaseous applications, benefits from the simplicity and flexibility offered by the non-conservative treatment of surface tension effects. It allows to successfully incorporate surface tension modeling, enabling the accurate representation of interfacial phenomena and enhancing the reliability of simulations in complex fluid dynamics scenarios. Although a fully conservative formulation is expected to give a more precise prediction of the capillary effects, the results presented later show that this simpler formulation yields reasonable results due to the diffuse nature of the interface.

The viscous stress tensor and the energy flux read

τ = - 2 3 ∇• (u)I + 2S, (19)
q = -λ∇T, (20)
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 , (21)
λ = α l λ l + α v λ v , (22)
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 21) and ( 22), µ 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 , (23)
with

Cp,g = N k=2 Y k C p,k N k=2 Y k . (24)
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, (25a)
∂ρu ∂t + ∇•(ρu ⊗ u + P I) = F σ + ∇•(µτ ), (25b)
∂ρE ∂t + ∇•((ρE + P )u) = u • F σ + ∇•(µτ • u) -∇•q, ( 25c
)
∂ρY k ∂t + ∇•(ρY k u) = 0. (25d)

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] Zheng | Lattice Boltzmann interface capturing method for incompressible flows[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.

Surface tension validation

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 initial species profiles are set with a hyperbolic tangent function, Φ, defined as,

Φ(ϕ(x, y)) = 1 2 1 + tanh 2 ϕ(x, y) ε , (26) ρ
(x, y) = Φ(x, y) ρ l + (1 -Φ(x, y)) ρ g , (27)
Y H 2 O l (x, y) = Φ(x, y) ρ l ρ(x, y) , (28)
Y N 2 (x, y) = 1 -Y H 2 O l (x, y) F + 1 , (29)
Y O 2 (x, y) = F Y N 2 (x, y), (30) (31)
where ϕ(x, y) reads

ϕ(x, y) = r -(x -x c ) 2 + (y -y c ) 2 , (32)
with r the droplet radius, ε the interface width, (x c , y c ) the droplet center coordinates, and F the ratio between O 2 and N 2 in pure air. The densities ρ l and ρ g of pure liquid and pure gas, respectively, are determined with the NASG equation of state. 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 . The theoretical pressure drop across the interface predicted by the Laplace law is

P l -P v = ∆P = σ r = 14.6 Pa. ( 33
)
The results of the simulation is plotted against the theoretical profile obtained with Eq. ( 33) in Fig. 1a, leading to a relative error of 0.5 %. The method employed to compute the error is based on the one adopted by [START_REF] Mirjalili | Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[END_REF], where it is computed in the all domain.

E P = i P (i) -Φ r -(x i -x c ) 2 + (y i -y c ) 2 σ r V i i Φ r -(x i -x c ) 2 + (y i -y c ) 2 σ r V i , (34)
where P (i) is the pressure at node i, Φ(x, y) is the initial hyperbolic tangent, (x i , y i ) are the coordinates of node i and V i is the nodal volume. A mesh convergence study is conducted in Fig. 1b. The droplet radius r is kept at 5 mm and the initial interface width is kept constant at ϵ = 0.4 mm, representing 4 cells in the interface width for the case depicted by Fig. 1a (i.e. the cell size is ∆x = 0.1 mm). The numerical grid has been coarsen and refined to target 1 to 8 cells in the interface width ϵ, by keeping it constant.

The performed simulations have demonstrated a clear relationship between mesh refinement and error reduction. As shown in Fig. 1b, the error level decreases linearly with finer mesh resolutions, indicating a first order convergent behavior.

Oscillating droplet test case

The oscillation of an elliptic viscous droplet is then considered, as shown in Fig. The initial species profiles are the same as in Sec. 4.1.1 but the ϕ function is here replaced by

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
ϕ(x, y) = 1 - x 2 a 2 - y 2 b 2 (35)
where 2a and 2b are respectively the ellipse width and height, and ϵ is the 255 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

t ⋆ = t ρ l r 3 eq σ . (36)
Theoretical studies have been proposed to predict the expected frequency of 260 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 where n = 2 is the most unstable mode. The frequency reads

ω 2 = (n 3 -n) σ (ρ l + ρ g )r 3 eq , (37)
f = ω 2π , (38)
where r eq is the equilibrium droplet radius. The numerical frequency is 73.03 Hz while the theoretical one computed with Eq. ( 38) is 71.61 Hz, resulting in 2 % 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 Eqs. ( 27) to [START_REF] Baer | A two-phase mixture theory for the deflagration-todetonation transition (ddt) in reactive granular materials[END_REF] 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 , ( 39
)
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 boundary condition. A pure liquid H 2 O droplet, initially traveling at a velocity of u 0 = 10 m.s -1 , is placed in a pure N 2 -O 2 domain, with thermodynamic conditions set to P 0 = 1 bar and T 0 = 293 K. The surface tension coefficient for this system is σ = 73 × 10 -3 N.m -1 . The domain is a square of 50 mm side. 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 liquid sheet

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 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 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 , (40)
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 , ( 41
)
where Re g is the gaseous flow Reynolds number defined as

Re g = ρ g U g H g µ g . ( 42
)
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 , (43)
and finally the dominant axial frequency is

f axi = U c λ axi . ( 44
)
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. [START_REF] Fuster | Instability regimes in the primary breakup region of planar coflowing sheets[END_REF], 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. ( 45) to [START_REF] Ma | Large-eddy simulations of transcritical injection and auto-ignition using diffuse-interface method and finite-rate chemistry[END_REF]:

u(y) = U l erf H l -y δ l + U i 1 + erf y -H l δ l × erf y δ l , (45)
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 , (46)
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 , (47)
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 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. 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 -
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. ( 48),

L = 12 H l √ M , (48)
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. ( 48). 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] Yang | Real-fluid injection modeling and LES simulation of the ECN Spray A injector using a fully compressible two-phase flow approach[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 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 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 The gas vorticity thickness δ g is adapted to the operating points in Eqs. ( 46) and ( 47) according to Eq. ( 41). 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. ( 44) (black solid line) where a prefactor C axi = 3.2 is considered. The same methodology as in Fig. 9a 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 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. [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] 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 two-phase flow configurations in a fully compressible vertex-centered framework, utilizing a diffuse-interface Multi-Fluid method and a 4-equation model closed with a Noble-Abel Stiffened-Gas equation of state. The vertex-centered HLLC-MUSCL numerical scheme has been implemented and validated on a number of academic test cases, demonstrating its ability to properly handle contact lines, shocks, and rarefaction waves. Additionally, the method effectively handles acoustics in the context of two-phase flow. Jacobian fluxes have been properly derived to allow the use of non-reflecting characteristic 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 ) -

Figure 1 :

 1 Figure 1: (a) Pressure difference profiles for simulation (red dashed line) and Laplace (black solid line). (b) Error on the Laplace pressure jump for different mesh resolutions. The liquidgas interface width and the droplet radius are kept constant at ε = 0.4 mm and r = 5 mm respectively.

2 .

 2 An elliptic shape droplet of radius a = 1.7 mm and b = 1 mm is initialized in a 10 × 10 mm domain of air (O 2 -N 2 ) at rest. Simulations have been performed on a 750 × 750 cells grid. The thermodynamic coefficients are given in Tab. 1. Thermodynamic features are P = 1 bar, T = 293 K and the surface tension coefficient is σ = 73 × 10 -3 N.m -1 . The liquid and gaseous dynamic viscosities are respectively µ l = 0.854 × 10 -3 Pa.s and µ g = 1.72 × 10 -5 Pa.s.

Figure 2 :

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

Figure 3 :

 3 Figure 3: Time evolution of the kinetic energy in the domain and the corresponding droplet shapes from the simulation (density field).

Figure 4 :

 4 Figure4: Superimposed backward and forward acoustic waves with fully non-reflecting boundary conditions. Red solid line: HLLC-MUSCL RK3 scheme[START_REF] Cheng | AVIP: a low temperature plasma code[END_REF], blue dashed line: TTGC scheme[START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF].

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 [80]. (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. (48) (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. (48). The blue area represents 40 isocontours of α l = 0.5 over 400 ms with n = 8 resolution representing the area swept by the liquid sheet interface.

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. (40) and (44).

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. (48).

  uH H -Λu 2 u(1 + Λ) -uΛ(Γ 1e c ) • • • -uΛ(Γ Ne c ) with H = E + P/ρ, Λ = ∂P ∂ρe ρY k and Γ = ∂ρe ∂ρY k P,ρY i̸ =k . Appendix A.1. Jacobian terms This section proposes the derivation of the Γ and Λ coefficient involved in 770 the jacobian matrix (Eq. (A.2)).From the mixture laws one can write

Figure B. 12 :

 12 Figure B.12: Liquid-gas shock tube results at t = 220 µs at CFL = 0.5.

  (a) Configuration of Eq. (B.4). (b) Configuration of Eq. (B.5). (c) Configuration of Eq. (B.6).

Figure B. 13 :

 13 Figure B.13: Density isocontours (60 lines) on a 600 × 600 grid at (a) t = 0.3 s, (b) t = 0.25 s, (c) t = 0.2 s. The numerical scheme is HLLC-M RK3 with CFL = 0.7. The wave configurations are the following: (a) 4 contact discontinuities, (b) 2 contact discontinuities, 2 rarefaction waves, (c) 1 shock wave, 1 rarefaction wave, 2 contact discontinuities.

Figure B. 14 :

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

  β = 1.0 HLLC-M β = 1.3 HLLC-M β = 1.7

Figure B. 15 :

 15 Figure B.15: Order of the schemes on the vortex test case after one convective time.

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.

	Raynal [80] (exp.) Marmottant and Villermaux [79] (exp.) Ben Rayana [82] (exp.) Matas et al. [83] (exp.) Fuster et al. [88] (VOF simulation -exp.) Agbaglah et al. [89] (VOF simulation) Present study	Strouhal (×10 -4 ) 6 7.5 7.6 5 8 -5 5.3 6.5

Table 4 :

 4 Liquid and vapour densities for each pressure and the associated momentum ratio M . 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.

	Pressure [bar] Liquid density ρ l [kg.m -3 ] Vapour density ρ g [kg.m -3 ] Momentum ratio M [-]	3 941.6 942.8 943.9 945.1 946.2 947.3 4 5 6 7 8 3.4 4.5 5.7 6.8 7.9 9.1 48 64 80 95 111 127
	operating ranges, hardly reachable with both experiments and incompressible solvers.

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

Finally one can write

. Isobaric dilatation and isotherm compressibility coefficients

The Isobaric dilatation α and isothermal compressibility β coefficients respectively read

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 mixture mass specific entropy reads

where C v,k , s 0 k , P ∞,k are specific constants of the NASG thermodynamics framework. The mixture mass specific entropy differential at constant composition reads

(A.28)

Casting Eq. (A.26) into Eq. (A.28) yields

Finally, the hyperbolic 4-equation system sound speed reads

Godunov numerical schemes [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] are usually implemented using a cell-centered framework. Fluxes are estimated on cell faces and conservative variables stored at the cell centre. The AVBP solver features a cell-vertex architecture, which offers an efficient way of handling numerical schemes like Lax-Wendroff or Taylor-Galerkin schemes [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF][START_REF] Lax | Systems of conservation laws[END_REF]. The Godunov scheme implemented in this work relies therefore on a node-centered formalism.

It considers the integral form of the transport equations, written in the compact form as

The integration is made on the nodal volume V i associated to node i, illustrated by the grey area in Fig. B.11. Fluxes are computed at all nodal volume faces, determined from the intersection of V i surface with all cells τ containing the node i, and noted E(i)

In Eq. (B.2), the first sum corresponds to a loop over all cells containing a given node, while the second sum corresponds to a loop over the faces of the nodal volume V i .

Surfaces integrals are computed by projecting the flux into the local basis of the edge as it is less computationally expensive than splitting the flux in 3 components in a fixed basis dimension and solving 3 Riemann problems ( [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF], Ch. 16). 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.

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 U L i = U n i -∆U n i /2 and U R i = U n i + ∆U n i /2 respectively on the left and on the right sides of node i along a given direction. For stability ∆U n i is computed with the limited approach of Sweby [START_REF] Sweby | High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws[END_REF]:

. One dimensional shock-tube

The solver is first validated on a two-phase shock tube configuration without surface tension (Sys. ( 2)). The present study considers a shock tube that is 1 m long and comprised of two chambers, which are separated by an interface , 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]:

(3.0, 0.75, 0.5, 1.0) x ≤ 0, y ≥ 0, (1.0,-0.75, 0.5, 1.0) x ≤ 0, y ≤ 0, (1.0, 0.75,-0.5, 1.0) x ≥ 0, y ≥ 0, (2.0,-0.75,-0.5, 1.0) x ≥ 0, y ≤ 0. 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 A mesh convergence study is conducted to assess the spatial order accuracy of the scheme from 32 × 32 to 2048 × 2048 cells grids. L 2 -error for the density, where the L 2 -error is defined as