Pierre Cordesse
email: pierre.cordesse@polytechnique.edu

Angelo Murrone

Thibaut Ménard

Marc Massot

Validation strategy of reduced-order two-fluid flow models based on a hierarchy of direct numerical simulations

Keywords: two-phase flow, reduced-order model, diffuse interface model, direct numerical simulation

come     Validation strategy of reduced-order two-fluid flow models based on a hierarchy of direct numerical simulations

Introduction

Rocket engines safety is one of the main priorities given to space industry. In a cryogenic combustion chamber, the multi-scale and multi-physical phenomena are very complex and their interaction still animate a research domain. In particular, the primary atomization plays a crucial part in the way the engine works, thus must be thoroughly studied to understand its impact on high frequencies instabilities. The latter have been encountered in the past and can lead to critical damages of the rocket. Even though experimentations must be conducted to enable simulation validation and to understand the observed physical phenomena, predictive numerical simulations are mandatory, at least as a complementary tool to understand the physic and even more to conceive new combustion chambers and predict instabilities they may generate in a given configuration.

In sub-critical condition, downside a air-assisted coaxial injector engenders three two-phase flow topologies: at the injector exit, the two phases are separated by a smooth interface. Downstream a polydisperse spray of droplets is carried by the gaseous phase. In between, shear stress caused by strong velocity gradients tears the liquid core apart and ligaments are formed. This process is called primary atomization. The ligaments get thinner and thinner until they break into droplets during the secondary atomization process. In this mixed region, the subscale physics and the topology of the flow are very complex. The typical range of flow numbers in sub-critical cryogenic cylindrical jet in real configuration are: for the liquid Reynolds number, Re L = 1e 5 , for the hy- drodynamic Weber number, W e H = 1e 5 . On the Mascotte test bench at ONERA, the flows numbers are both an order of magnitude less. As for the ambient pressure, p, it reaches p = 25 bar.

Direct numerical simulations (DNS) in real configuration of such engines are still out of reach, CPU needs being too high due to the high Reynolds and Weber numbers, predictive numerical tools using reduced-order models must be developed. However great care must be taken on the choices of these models in order to both have solid mathematics properties and lead to predictive simulations after a validation process.

The strategy retained to perform a numerical simulation of the primary atomization from the injection to the combustion of the spray is to couple a Kinetic-Based-Moment Methods (KBMM) to describe the disperse flow as in [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods[END_REF] with a diffuse interface model for the separated phases and the mixed zone. Among the hierarchy of diffuse interface model, we focus on the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-to-Detonation Transition (DDT) in reactive granular materials[END_REF] accounting for full disequilibrium of the phases. In the present work, we propose a validation of these reduced-order model relying on DNS on an hierarchy of specific test cases. In [START_REF] Cordesse | Comparative study of jet atomization simulations: direct numerical simulations and diffuse interface models coupled with kinetic-based moment methods[END_REF], we have started with an air-assisted water atomization using a coaxial injector, which also provides experimental results from the LEGI test bench. The comparison has shown good agreements in terms of liquid core length and important CPU gains between the seven equation model implemented in the CEDRE code [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF] and the DNS results from the ARCHER code. It has also shown the limits of diffuse interface models to capture complex liquid structures such as ligaments, rings or deformed droplets and encourages to add a sub-scale description of the interface dynamics through geometric variables such as the interfacial area density, the mean and Gaussian curvatures as proposed in [START_REF] Cordesse | Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus[END_REF].

We propose here a second test case, an air-assisted water atomization using a planar injector, which reproduces in terms of Weber and Reynolds number the liquid sheet flowing out a swirling atomizer used in agricultural applications [START_REF] Belhadef | Pressure-swirl atomization: Modeling and experimental approaches[END_REF], Re L = 1.5e 4 and W e R = 4.0e 2 and p = 1 bar. The DNS results obtained thanks to the code ARCHER [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows[END_REF] denotes that this test case offers an atomization regime, which makes it complementary with the first test case in order to eventually validate our reduced-order models on a cryogenic coaxial injection.

In the first part, we will recall the mathematical properties of the models solved. Then, we explain the numerical methods implemented in the two code, CEDRE and ARCHER. Finally, we present the results obtained by both simulations and compare them.

Mathematical modelling

At the top of the hierarchy of diffuse interface models [START_REF] Drui | A hierarchy of simple hyperbolic two-fluid models for bubbly flows[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF] stands the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-to-Detonation Transition (DDT) in reactive granular materials[END_REF], also called the seven equation model, accounting for full disequilibrium of the phases. This first-order non-linear non-conservative system of partial differential equations is composed of a mass, momentum and energy equation for each phase and a seventh equation on the volume fraction α to reconstruct the interface. The extended form proposed in [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] introduces an interfacial pressure, p I , and an interfacial velocity, v I , that need to be closed. From this seven-equation model, the instantaneous relaxation of the pressures and the velocities leads to the five-equation model [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF] and relaxing also instantaneously the temperatures, one obtains the compressible multi-species Navier-Stokes equations referred as the fourequation model. These three models define a hierarchy of diffuse interface model and each of them are hyperbolic and appropriate for reactive two-phase flow or interface problems.

In cryogenic applications, the pressure of the phases may be considered to relax instantaneously, but due to the strong velocity and temperature gradients at the interface, it is unrealistic to assume hydrodynamic and thermal instantaneous relaxations. Whereas single velocity models are sufficient in the separated zone with an adapted resolution mesh, it is not the case in the mixed region and thus fail at predicting the atomization process. Therefore the Baer-Nunziato model appears as the best two-velocity model candidate for the present study.

Nonetheless, stemming from rational thermodynamics, the macroscopic set of equations can not be derived from physics at small scale of interface dynamics and thus require closure of interfacial pressure and velocity as well as to postulate the thermodynamics. The theory on the existence of a supplementary conservative equation to first-order order system of partial differential equation including non-conservative terms [START_REF] Cordesse | Mathematical modelling of multi-phase flow using entropy symmetrization[END_REF][START_REF] Cordesse | Supplementary conservative law for non-linear systems of PDEs with nonconservative terms: application to the modelling and analysis of complex fluid flows using computer algebra[END_REF] has brought about supplementary conservative equations together with constraints on the interfacial quantities and the definition of the thermodynamics of the mixture. In this work, we use a non-miscible fluid thermodynamic, where the mixture entropy is defined by the mass averaged sum of the phasic entropies with no mixing effect, together with an closure proposed in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF] obtained through a discrete element method approach [START_REF] Saurel | A multiphase model with internal degrees of freedom : application to shock-bubble interaction[END_REF]. To the convection part described in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF], we let the two phases relax towards an equilibrium state thanks to sources terms. Only mechanical, and hydrodynamic relaxations are accounted for in the present study. They are modelled as in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF] through a pressure difference term and a drag term. Each term is driven by a characteristic time, which can be either finite to account for disequilibrium, or in the limit of a zero characteristic time, to account for an instantaneous relaxation and kill all the disequilibrium.

The direct numerical simulation solves the incompressible Navier-Stokes equations and takes into account surface tension. Then, the interface is implicitly derived by the zero of a level set of a scalar function and its motion is captured by the transport of the level set function at the hydrodynamic speed, combined with a projection technique [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompress-10[END_REF]. To guarantee the conservation of mass, we also transport the volume of fluid function (VOF), defined as the liquid volume fraction expressed in terms of the level set function. The benefits of a VOF formulation coupled with a Level Set function are to conserve mass and to have access to geometrical properties of the interface.

Numerical methods

The numerical methods employed to solve the Baer-Nunziato model are implemented in the multiphysics computational fluid dynamics software CEDRE [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF] working on general unstructured meshes and organized as a set of solver [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF]. The solver SEQUOIA is in charge of the diffuse interface model whereas the solver SPIREE takes care of the droplets once atomized.

A Strang splitting technique is applied on a multi-slope HLLC with hybrid limiter solver [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF][START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] to achieve a time-space second-order accuracy on the discretized equations. The issue encountered when discretizing the nonconservative terms is tackled in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF] by assuming (1) the interfacial quantities p I and u I to be local constants in the Riemann problem, (2) the volume fraction to vary only across the interfacial contact discontinuity u I . As a result, the non conservative terms vanish, u I and p I are determined locally by Discrete Equation Method (DEM) [START_REF] Saurel | A multiphase model with internal degrees of freedom : application to shock-bubble interaction[END_REF] at each time step and stay constant during the update. Thus, phases are decoupled the system splits into two conservative sub-systems to which we apply the multi-slope HLLC with hybrid limiter solver.

Depending on the application, the relaxations are assumed either instantaneous or finite in time. In the present test case, it is reasonable to assume a instantaneous pressure relaxation but one need to consider a finite velocity relaxation since the interface dynamic is mainly driven by the shear stress induced by a high velocity difference between the phases at the injection.

To obtain the relaxed pressure, one needs to solve an simple ODE. Since the characteristic time is taken to be infinitely small, the problem reduces to apply an iterative procedure as a Newton method to solve a second order equation and obtained a single equilibrium pressure. Detailed of the equation can be found in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF]. As for the velocities, since we want to account for finite relaxation time, the associated ODE takes the form:

∂ t v d - A o v v d = 0, with A o = α o 1 ρ o 1 + α o 2 ρ o 2 α o 1 ρ o 1 α o 2 ρ o 2 , (1)
where v d is the slip velocity, v is the characteristic relaxation time and superscript o denotes the state before relaxation and α k , ρ k denote respectively the volume fraction and the partial density of phase k. A first numerical approach is to fix a remaining slip velocity ratio target at each computational time step ∆t. It defines the characteristic relaxing time as

v A o = ln v d (∆t) v o d ∆t. (2)
An instantaneous velocity relaxation is in practice also possible and manipulating the ODE leads to a unique relaxed velocity, which is the mass weighted average of the two velocities before relaxing.

As for the DNS, as mentioned by Rudman [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF], it is recommended to solve Navier-Stokes equations in their conservative form to ensure consistance between mass, VOF, and momentum fluxes. Because a staggered grid is adopted, VOF and velocity have not the same control volume making it difficult to stay consistant. To avoid the finest grid for the VOF transport proposed in [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF] to preserve consistency, the method developed in [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows[END_REF] is used here and permit to reduce time computation.

The Navier-Stokes equations are implemented in the code ARCHER [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows[END_REF]. A projection method is employed and physical properties (viscosity and density) are expressed in term of both VOF and Level Set. The temporal integration is performed through a second-order Runge-Kutta scheme. The discretization of convective term is achieved with WENO 5 scheme [START_REF] Jiang | Efficient Implementation of Weighted {ENO} Schemes[END_REF]. For the viscosity term, we retain the method presented by Sussman [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. The Ghost-Fluid [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF] method is employed to take into account surface tension, treated as a pressure jump.

To compute correctly mass fluxes at the interface (consequently momentum fluxes) and to preserve consistency between both flow and interface solvers, VOF and Level set transport are performed with a CLSVOF algorithm [START_REF] Sussman | A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows[END_REF] at each step of Runge Kutta scheme.

Description of the configuration

To attest the reduced-order model introduced in Section 1 with direct numerical simulation, we propose an air-assisted flat liquid sheet at high Reynolds and Weber number, inspired from a swirling injector used for agriculture [START_REF] Belhadef | Pressure-swirl atomization: Modeling and experimental approaches[END_REF]. This test case contributes to the validation strategy based on an hierarchy of direct numerical simulation test cases that eventually wants to reproduce an air-assisted coaxial atomization of liquid in cryogenic rocket engines and started in [START_REF] Cordesse | Comparative study of jet atomization simulations: direct numerical simulations and diffuse interface models coupled with kinetic-based moment methods[END_REF].

Air-assisted flat liquid sheets have been widely experimentally investigated such as in [START_REF] Dumouchel | On the experimental investigation on primary atomization of liquid streams[END_REF]. The behaviour of such flow is mainly driven by four parameters. The ratio V G /V L and the momentum flux ratio M determine the breakup regime. Then the width of liquid injection, D L and the relative gaseous Weber number, W e R , influence the breakup length, the breakup frequency and the liquid sheet vibrating frequency [START_REF] Dumouchel | On the experimental investigation on primary atomization of liquid streams[END_REF]. M and W e R are defined in Equation ( 3).

M = ρ G V 2 G ρ L V 2 L , W e R = ρ G (V G -V L ) 2 D L σ . (3)
Other classical flow parameters exist such as the liquid Reynolds number, Re L = ρ L V L D L /µ L , the gas Reynolds number,

Re G = ρ G V G (D G -D L )/µ G , the liquid Weber Number, W e L = ρ L V 2 L D L /σ , the gas Weber Number, W e G = ρ G V 2 G (D G -D L )/σ or the aerodynamic Weber number, W e H = ρ G V 2
G D L /σ, but they are preferred in cylindrical jet applications.

The present configuration has been obtained by adjusting the flow parameters in order to optimize the flapping amplitude and to guarantee a wide liquid core. The momentum flux ratio is fixed to M = 1.88 and the relative gaseous Weber number to W e R = 403. No Break-up regimes in the parameter spaces Re L -W e exists for planar jet, as what exists for coaxial jet [START_REF] Lasheras | Liquid Jet Instability and Atomization in a Coaxial Gas Stream[END_REF]. However, high relative velocity, as it is the case here, usually leads to strong sinusoidal oscillations of the liquid flow inducing an high spray angle and short sheet breakup length, called "flag-effect" [START_REF] Dumouchel | On the experimental investigation on primary atomization of liquid streams[END_REF]. The DNS results attests these predictions and clearly shows an atomized spray.

The simulated domain described in Figure 2 shows the boundary conditions of the simulation. 1 is the liquid injection plan, 2 is the gas injection plan, 3 and 3 are periodic plans, 4 , 5 and 6 are outflow plans. The liquid injection height, D L , the gas injection height, D G , the x-y-z length of the box, L x , L y , L z , are given in Table 1. Figure 3 defines the velocity profile of the gas and the liquid, which are symmetric with respect to the x-axis and y-invariant. The gas velocity profile V G given in Equation ( 4) is typical for turbulent pipe flow [START_REF] Schlichting | Boundary -Layer Theory[END_REF].

V G = V avg G 7 6 2|z| -D L D G -D L 1 6 + V L . (4)
The average gas velocity V avg G is 65 m/s. An offset equal to V L = 1.5 m/s ensures the continuity of the velocity profile at the injection plan. The domain is initially filled with a liquid sheet in the x -y plan, as thick as the liquid slit as shown on Figure 2. Table 2 states the fluid properties in terms of density ρ, surface tension coefficient, σ, and viscosity, µ. The ARCHER simulations are performed on a Cartesian mesh 512 × 128 × 512 with a cell size equal to ∆x = 3.125 10 -5 m, so a total of 101 millions of faces whereas CEDRE mesh is composed of tetrahedral cells with a total of 148k faces.

In order to compare the results of the DNS to the results obtained with CEDRE, one must consider the fact that the DNS solver is incompressible, thus there is no acoustic impacting the liquid sheet and its density is constant. To restrain the acoustic role in the CEDRE compressible solver, we have enlarged the computational domain by a factor 5 in the x and z direction and meshed it with a very coarse mesh to avoid any reflecting waves as shown in Figure 4a. The minimum cell size is ∆x = 1.0 10 -4 m and can be found along the liquid slit as one notices on Figure 4b.

Furthermore the use of a compressible thermodynamics in CEDRE through the Stiffened-Gas equation of state makes 2 and in practice, the liquid density almost stays constant as shown on Figure 5.

V L V G D L 2 D G 2 z x
The convection time t conv of the system, defined as

t conv = L x V L , (5)
equals t conv = 10.7 ms. The CEDRE simulation has run t sim = 27 ms. which corresponds to approximately three convective times. The simulation information are summarized in Table [START_REF] Cordesse | Comparative study of jet atomization simulations: direct numerical simulations and diffuse interface models coupled with kinetic-based moment methods[END_REF]. The total CPU cost is defined as the prod- uct of the number of processors, N proc , times the computational time for a given simulation time, t sim . The CPU difference is notable, ARCHER is 33 more costly than CEDRE.

Results and discussion

Before comparing the results from CEDRE to the DNS, we propose to analyze qualitatively and quantitatively the results obtained with CEDRE.

On Figure 6, the instantaneous liquid sheet axial length is plotted with respect to the simulation time t sim . The liquid sheet length is defined with respect to a liquid volume fraction isovalue. Left blue y-axis encompasses the whole range of the liquid sheet at the isovalue, α l = 0.95. Initially at a length L lc = 35 mm, the liquid sheet rapidly gets smaller and breaks at t sim = 2.5 ms to reach an average value of L lc = 0.8 mm emphasized on the right y-axis where two isovalues, α l = 0.91 and α l = 0.99 are used to capture the liquid sheet length. The transient regime lasts around t = 3 ms during which the initial liquid sheet breaks, the downside part being carried away by the fast flowing air. Then the planar jet starts flapping steadily.

A sensor placed near the injection slit, at (x, y, z) = (2, 0, 0) mm, has recorded the liquid volume fraction over

Many experimental correlations exist for this longitudinal instability, and as observed for example in [START_REF] Lozano | Longitudinal instabilities in an air-blasted liquid sheet[END_REF], the frequency is highly depending on velocity profile at the interface at the edge of the slit. It means that this frequency depends mainly on the mesh resolution and it is very unlikely to find the same frequency for different meshes, and thus to compare the results of CEDRE simulation with the DNS.

The two-phase flow model offers a non-instantaneous velocity relaxation of the phases. This disequilibrium is relevant in strong velocity gradient regions such as close to the liquid sheet injection plan. Physically, the liquid is expected to be accelerated by the gas, but not instantaneously. Therefore it is physically wrong to assume an instantaneous relaxation time for the velocities of the phase. The maximum slip velocity norm over time is plotted in Figure 9a with the averaged slip velocity norm in Figure 9b, together with the averaged liquid volume fraction isolines. Interestingly, the slip velocity is not negligible in the flapping zone and is in average 0.3 m/s. Numerically, after the velocity relaxing procedure, the remaining slip velocity is only ten percent of the slip velocity before relaxing procedure. It means that in this region, the slip velocity before relaxation is about 3 m/s. These results attest the robustness of the solver able to handle velocity slip but also the interest of a two velocity model.

We now propose to compare the liquid sheet obtained with ARCHER and CEDRE at given simulation time. The Level Set function of the DNS permits an exact reconstruction of the interface whereas for the diffuse interface model, the interface lays in the region where the volume fraction varies from α ≈ 0 to α ≈ 1. Consequently, in Figure 10 we have superimposed the solved interface of the DNS to a volume rendering of the liquid volume fraction and a single liquid volume fraction isosurface. α H2O In Figure 10, we distinguish two regions. Close to the injector slit, the diffuse interface model matches the DNS results. The CEDRE mesh used in this region prevents the interface from diffusing too much such that the interface captured by the DNS lays in the volume rendering of the liquid volume fraction of CEDRE. The breakup length of the liquid sheet matches qualitatively. Downstream the liquid core of the planar jet in the mixed region, the DNS shows complex liquid structures that rapidly become a dense spray. The diffuse interface model is ineffective in capturing these effects and we only see a red cloud corresponding to a low liquid volume fraction isovolume. The reasons of such a discrepancy between CEDRE and the DNS is twofold. First the CEDRE mesh is too coarse compared to the DNS mesh. Second, the volume fraction, alone, is not sufficient to describe such complex interface dynamics. We would need to add a subscale modelling of the interface dynamics as proposed in [START_REF] Cordesse | Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus[END_REF] where a supplementary geometric variable, the interfacial density area, is solved. Finally, this comparison gives an interesting interpretation of the diffuse interface models. The volume fraction is not enough to reconstruct the whole dynamic of the interface but attests to the presence or the absence of liquid.

Conclusions

In the present work, we proposed to pursue the evaluation of reduced-order models to perform predictive simulations of the primary atomization in cryogenic rocket engines. As a baseline comparator, we relied on the DNS results of an hierarchy of specific test cases. We originally started with an air-assisted water atomization using a coaxial injector, in [START_REF] Cordesse | Comparative study of jet atomization simulations: direct numerical simulations and diffuse interface models coupled with kinetic-based moment methods[END_REF], which also provides experimental results from the LEGI test bench. Here we proposed an air-assisted water atomization using a planar injector representative of swirling atomizers found in agricultural applications [START_REF] Belhadef | Pressure-swirl atomization: Modeling and experimental approaches[END_REF].

The reduced order model was a diffuse interface model offering an non-instantaneous velocity relaxation and was implemented using a Strang splitting technique applied on a multislope HLLC solver [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium twophase flows[END_REF][START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] in the CEDRE code of ONERA [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF]. The DNS results were obtained using the incompressible Navier-Stokes equations together with a CLSVOF [START_REF] Sussman | A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows[END_REF] thanks to the code ARCHER [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows[END_REF]. This test case offers an atomization regime with a dense spray, which makes it complementary with the first test case in order to eventually validate our reduced-order models on a cryogenic coaxial injection.

The comparison has shown good agreements in terms of liquid sheet length near the slit, together with an important CPU However, it has shown the limits of diffuse interface models to capture the complex liquid structures from the liquid sheet in the mixed region due to the coarse CEDRE mesh, but also due to the limitation of the volume fraction to describe the dynamic of the interface. It suggests to add a sub-scale description of the interface dynamics through geometric variables such as the interfacial area density, the mean and Gaussian curvatures as proposed in [START_REF] Cordesse | Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus[END_REF].

In future works, we will first aim at defining a physical dynamic velocity relaxation time. Then, we will add the subscale description of the interface dynamics to obtain a better solution in the mixed region and compare it to the DNS results.

Figure 1 :

 1 Figure 1: Jet atomization flow topologies

Figure 3 :

 3 Figure 3: Injector schematic and velocity profiles

10 th

 10 International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019 (a) CEDRE mesh (b) Refined DNS box

Figure 4 :

 4 Figure 4: CEDRE mesh of the configuration

Figure 5 :

 5 Figure 5: Instantaneous liquid density ρ L over time, maximum value, minimum value.

Figure 6 :

 6 Figure 6: Liquid sheet length L lc over time at isovalue: α 2 = 0.95, α 2 ∈ [0.91, 0.99].

10 th

 10 International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019

Figure 7 :

 7 Figure 7: Liquid volume fraction at sensor (x, y, z) = (2, 0, 0) mm over time.

(a) 23

 23 .70 ms (b) 23.75 ms (c) 23.80 ms (d) 23.85 ms (e) 23.90 ms (f) 23.95 ms (g) 24.00 ms (h) 24.05 ms (i) 24.10 ms (j) 24.15 ms (k) 24.20 ms (l) 24.25 ms

Figure 8 :

 8 Figure 8: Instantaneous liquid volume fraction near the injector, α l = 0.1 α l = 1.0, probe at (2, 0, 0) mm

10 th

 10 International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019 (a) Maximum value over time (b) Time average

Figure 9 :

 9 Figure 9: Contour of slip velocity norm over time, v d = 0 m/s v d = 0.65 m/s, and averaged liquid volume fraction isolines, range [0 : 1] by 0.1 increment. Box size is L x /4 × L z /4.

  10 th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019

Table 1 :

 1 10 th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019 Test case configuration -dimensions units D L D G L x L y L z

		z		5
		2
				3	3	6
	1	y 0				x
		2
				4
	Figure 2: Test case configuration and initial boundary con-
	ditions
	(mm)	1	16	16	4	16

Table 2 :

 2 Physical properties of liquid and gasPhase ρ (kg/m 3 ) p (M P a) σ (N/m) µ (P a.s)

	Liquid	1000	0.1	0.01	0.0001
	Gas	1	0.1	0.01	0.0001

Table 3 :

 3 Simulation costs comparison for 9 ms t sim [ms] N proc CPU cost [h]

	CEDRE	9.0	700	8.54 10 3
	ARCHER	9.0	552	2.83 10 5

th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil,[START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompress-10[END_REF]

Acknowledgments

The support of the French Government Space Agency (CNES) and the French Aerospace Lab (ONERA) through a PhD grant for P.Cordesse, the help of L.H. Dorey (ON-ERA) and M. Théron (CNES) are gratefully acknowledged. We also thank A. Remigi, CNRS CORIA UMR 6614. Simulations have been successfully conducted using the CEDRE computational fluid dynamic software on the ONERA cluster.


domain: spi phys

text:
Angelo Murrone
email: angelo.murrone@onera.fr

Aymeric Boucher
email: aymeric.boucher@onera.fr

Pierre Cordesse
email: pierre.cordesse@onera.fr

A FIVE EQUATION MODEL FOR THE SIMULATION OF THE TWO-PHASE FLOW IN CRYOGENIC COAXIAL INJECTOR

Keywords: liquid rocket engines, numerical simulation, diffuse interface modelling, coupling, primary atomization

In this paper, we propose models and methods for the simulation of two-phase flows in Liquid Rocket Engines (LRE) under subcritical conditions. The numerical strategy consists into coupling models dedicated to different topologies. Actually, we propose a five equation diffuse interface model for the treatment of the dense "separated two-phase flow" near the injector and an Eulerian kinetic based model for the "dispersed two-phase flow" in the chamber. We derive a novel formulation of the 5 equation system to build a robust HLLC type scheme. Then we use a fully Eulerian coupling strategy to take into account for primary atomization. We first run classical test cases in order to validate the numerical methods. Then a simulation on a test case representative to one coaxial injector is performed under subcritical conditions.

INTRODUCTION

To help the development and enhancement of launcher propulsion systems and to ensure their reliability, great care must be taken to study all the physical processes involved in the combustion chamber. Concerning liquid propulsion, this is a major issue which concerns a large variety of multiscale phenomena such as evaporation, combustion, turbulence, primary atomization and break-up droplets. Moreover high frequency (HF) instabilities resulting from the complex interaction between propellants injection, flame dynamics and acoustic modes have been encountered in the recent past and have to be studied in a more comprehensive way. In this study, we focus on injection and primary atomization which play a key role to control the combustion under transient operating conditions and have also an influence on HF instabilities. As a result, one of the long term objectives for the launcher industry is to better understand under which conditions the atomization process tends to generate or amplify these HF instabilities, and under which conditions it has the opposite effect. We expect that CFD tools combined with experiments will eventually help to achieve this goal.

Figure 1 represents the different phenomena occurring at the exit of one coaxial cryogenic injector within a liquid rocket engine (LRE) operating in the subcritical regime. Under such conditions, the oxidizer (oxygen) exits the injector in liquid form whereas the fuel (hydrogen) is in gaseous form, which consequently entails a twophase flow. Initial breakup of the bulk liquid, which is called primary atomization, occurs near the injector exit due to the strong difference of velocity between the two phases. The subsequent disintegration into child droplets due to secondary fragmentation occurs downstream. Following the classification of Ishii [START_REF] Ishii | Thermo-fluid dynamics of twophase flow. 2 nd Edition[END_REF] this two-phase flow covers the entire range of volume fraction. In the dense region of the liquid jet, the atomization results from interaction between liquid oxygen (LOx) and gaseous hydrogen (GH 2 ) phases. Ligaments start growing from the liquid core because of Kelvin-Helmholtz and Rayleigh-Taylor interfacial instabilities [START_REF] Marmottant | On spray formation[END_REF]. The ligaments thus formed are unstable because disruptive forces exceed the liquid surface tension and viscous forces. This results in a spray of small LOx droplets with final stable sizes, mainly spherical, which are dispersed by the turbulent gas flow and finally vaporized to feed the combustion with GH 2 . Thus the gas phase is a mixture of multiple chemical species, namely hydrogen, vaporized oxygen and combustion products. Eventually, the resulting high-enthalpy combustion products exhaust through a nozzle at supersonic speed, thereby providing the required thrust.

Figure 1: schematic representation of a coaxial cryogenic injector

A lot of experimental and theoretical studies were conducted on the atomization process but the details of the breakup cascade are still not fully understood. There are no sufficient results to provide parameters like the expansion angle, the penetration depth or the droplet size distribution. Regarding experimentation, the main reason lies in the droplets cloud surrounding the liquid core region which block the access to optical rays. Moreover, it is difficult to view phenomena because the length and time characteristics are very small. Nevertheless, with recent improvements in optical ray and x-ray imaging techniques, a detailed analysis of gas liquid interface has been provided in [START_REF] Marmottant | On spray formation[END_REF]. For the MASCOTTE cryogenic injector [START_REF] Vingert | MASCOTTE: a research facility for high pressure combustion of cryogenic propellants[END_REF], the spray was recently investigated with a high speed camera in a backlighting optical configuration [START_REF] Murrone | Atomization of cryogenic rocket engines coaxial injectors[END_REF]. On the other hand, experiments concerning diesel injection are gaining importance [START_REF] Wang | Ultrafast x-ray study of dense liquid jet flow dynamics using structure tracking velocimetry[END_REF], [START_REF] Crua | High-speed microscopic imaging of the initial stage of diesel spray formation and primary breakup[END_REF], [START_REF] Osta | Study of the effects of the injector length/diameter ratio on the surface properties of turbulent liquid jets in still air using x-ray imaging[END_REF]. We expect that in a near future, spatial resolution at the sub-micrometer level may provide databases for cryogenic rocket engines but also for Direct Diesel Injection.

Concerning numerical simulation, computational resources give us the possibility to view the atomization process more accurately than before and DNS are gaining importance [START_REF] Fuster | Simulation of primary atomization with an octree adaptive mesh refinement and VOF method[END_REF], [START_REF] Desjardins | Detailed numerical investigation of turbulent atomization of liquid jets[END_REF], [START_REF] Herrmann | A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure[END_REF], [START_REF] Herrmann | On simulating primary atomization using the refines level set grid method[END_REF]. A lot of works are thus based on the ARCHER software [START_REF] Menard | Coupling level set/vof ghost fluid methods : validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF] which combines both Level Set and VOF methods. In [START_REF] Duret | Improving primary atomization modeling through DNS of twophase flows[END_REF], the DNS of a two-phase flow provides very promising results. Nevertheless, the finest scales of two-phase flows are not known a priori, contrary to the Kolmogorov scale, and moreover cannot be resolved even if the mesh size is below the micrometer scale. Other recent works of DNS can also be mentioned [START_REF] Bouali | DNS analysis of small scale turbulence scalar interaction in evaporating two phase flows[END_REF] but it seems that DNS is not mature enough for industrial configurations.

Leaving DNS aside, we can find in the literature various numerical works based on the RANS or LES formalisms where the primary atomization is not taken into account. In this kind of approach, the parameters giving the initial properties of the spray such as droplets size and velocity distribution are given as presumed input data in the numerical simulation. Consequently we cannot rely on these models to simulate complex unsteady phenomena in which primary atomization plays a key role (as it is for instance highly suspected with HF instabilities in LRE). We need to derive numerical models, within the RANS or LES frameworks, valid from the primary atomization of the liquid jet, to the computation of the spray without having to impose its properties as input parameters.

This quest has been initiated in the seminal work of [START_REF] Vallet | Contribution à la Modélisation de l'Atomisation d'un Jet Liquide Haute Pression[END_REF], [START_REF] Vallet | Modélisation Eulerienne de l'atomisation d'un jet liquide[END_REF], [START_REF] Vallet | Development of an Eulerian model for the atomization of a liquid jet[END_REF] and was pursued by many subsequent works, most of them applied to Direct Diesel Injection in internal combustion engines but also sometimes to cryogenic rocket engines [START_REF] Pourouchottamane | Modélisation des brouillards denses pour la combustion cryotechnique[END_REF], [START_REF] Jay | Modélisation de la combustion diphasique au moyen de bilans d'aire interfaciale et de surface de flamme[END_REF], [START_REF] Jay | Combined surface density concepts for dense spray combustion[END_REF], [START_REF] Meyers | Modélisation de la combustion cryotechnique avec prise en compte de l'atomisation primaire du jet d'oxygène liquide[END_REF], [START_REF] Demoulin | A new model for turbulent flows with large density fluctuations: application to liquid atomization[END_REF]. All these works involve a sophisticated surface density equation (see [START_REF] Morel | On the surface equations in two phase flows and reacting single phase flows[END_REF] for theoretical analysis of such an equation) that contains source terms for creation and destruction of interface area. These source terms are closed within the RANS framework of turbulence, thanks to sub-grid models derived from theoretical, experimental or DNS results [START_REF] Duret | Improving primary atomization modeling through DNS of twophase flows[END_REF]. A similar surface density equation is used in an LES context in [START_REF] Chesnel | Simulation aux Grandes Échelles de l'Atomisation[END_REF], [START_REF] Chesnel | Subgrid analysis of liquid jet atomization[END_REF], [START_REF] Chesnel | Large eddy simulation of liquid jet atomization[END_REF] together with accurate tracking interface methods. Another interesting work in this context is the recent one presented in [START_REF] Devassy | Atomization modelling of liquid jets using a two-surface-density approach[END_REF] which deals with the liquid jet atomization under Direct Diesel engine conditions, but without considering combustion. In the latter, the authors retain a 7 equation diffuse interface model [START_REF] Baer | A twophase mixture theory for the deflagration-todetonation transition (DDT) in reactive granular materials[END_REF] together with relaxation procedures for pressure, velocity [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] and temperature [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF], which they use for both the liquid core and the spray droplets. A new atomization model is also formulated using two surface density equations and the closure terms are based on RANS simulations, the description of the turbulence being necessary for both liquid and gas phases. Using two surface density equations is a novelty since all previous works necessitated a unique surface density equation. The two surface density equations are transported with the same interface velocity; one for the dense "separated" phases and another one for the "dispersed" phase.

The objective is to improve both the primary atomization and the secondary breakup modelling. Besides, in [START_REF] Andreini | Development of a turbulent liquid flux model for Eulerian-Eulerian multiphase flow simulations[END_REF], a quasi-multiphase Eulerian (QME) solver is implemented in OpenFOAM and applied to a jet in crossflow. Comparisons with experiments and DNS of [START_REF] Herrmann | A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure[END_REF], [START_REF] Herrmann | On simulating primary atomization using the refines level set grid method[END_REF] are presented. The innovative second order closure for the slip velocity is based on transport equations for momentum, volume fraction, surface density and liquid flux. It is therefore an intermediate approach between fully multiphase and mixture approaches since it holds for a large range of liquid volume fraction going from dense "separated" to dilute "dispersed" flow.

Based on the above-mentioned works, the idea of coupling the "RANS equations plus surface density equation" strategy with a Lagrangian description for the spray was introduced by work presented in [START_REF] Beau | Modélisation de l'atomisation d'un jet liquide[END_REF], which has led to the development of the dedicated ELSA (Eulerian Lagrangian Spray Atomization) code. The idea is to reach a better overall accuracy by allowing the spray to be described with its own variables (velocity, size, temperature), rather than describing it only through the surface density equation and the mean velocity and temperature of the whole subgrid liquid-gas mixture. Following this strategy, the subsequent works of [START_REF] Lebas | Modélisation Eulérienne de l'Atomisation Haute Pression -Influences sur la Vaporisation et la Combustion Induite[END_REF], [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF] made the source terms depend on whether the topology of the twophase flow is "dispersed" or "separated", and they also improved the definition of the equilibrium Weber number. However it seems that the strong coupling between Eulerian and Lagrangian methods induces some difficulties [START_REF] Demoulin | Toward using direct numerical simulation to improve primary break-up modeling[END_REF]. In [START_REF] Herrmann | A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure[END_REF], an efficient parallel multi-scale coupling procedure between an Eulerian level set method tracking interface and a Lagrangian description of small scales has been applied successfully to a turbulent liquid jet under Diesel engine conditions, but again with some remaining numerical difficulties. Indeed, several difficulties may arise when dealing with Eulerian -Lagrangian methods. The first one is due to the lack of robustness of the Lagrangian methods in the case of a strong two-way coupling. The second one is due to the statistical interpretation of numerous numerical particles produced in the atomization area. It seems also that taking into account for the compressibility of the fluid with the tracking interface method is not an easy task. In the Eulerian-Eulerian methods, the major drawbacks are that numerics are possibly tricky when solving the spray equations (discontinuities, weak hyperbolicity...), and also that the description of the phenomena at the droplet scale may be more complex than in the Lagrangian formalism. But one major advantage against interface tracking and Lagrangian methods is that models are general, fully compressible and adapted to the two-way coupling simulations. Moreover, Eulerian-Eulerian methods are also well adapted to parallel and time implicit computations. This is why we follow in this paper a fully Eulerian coupling strategy between "separated" and "dispersed" two-phase flow models for the simulation of primary atomization, introduced in [START_REF] Murrone | Atomization of cryogenic rocket engines coaxial injectors[END_REF], [START_REF] Murrone | Eulerian coupling of two-phase flow models for the Large Eddy Simulation of the atomization in cryogenic combustion chamber[END_REF], [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF]. It combines a kinetic-based model for the description of the spray with a diffuse interface model to describe both the "separated" and "mixed" two-phase flows. This work has been conducted in the framework of the multi-physics CEDRE platform [START_REF] Refloch | CEDRE Software[END_REF], [START_REF] Murrone | Numerical Modeling of Dispersed Two-Phase Flows[END_REF] developed at ONERA and applied to the cryogenic MASCOTTE facility [START_REF] Vingert | MASCOTTE: a research facility for high pressure combustion of cryogenic propellants[END_REF] representative of one coaxial injector of a liquid rocket engine.

In a previous works [START_REF] Murrone | Eulerian coupling of two-phase flow models for the Large Eddy Simulation of the atomization in cryogenic combustion chamber[END_REF], we have proposed a coupling strategy between a 4 equation diffuse interface model and an Eulerian kinetic model for the spray. The fully Eulerian coupling between "separated" and "dispersed" two-phase flow solvers has been applied to the simulation of the MASCOTTE [START_REF] Gicquel | Caractérisation expérimentale en combustion d'un brouillard LOX-GH2 issu d'un injecteur coaxial (campagne d'essais 1997 à 1 MPa)[END_REF], [START_REF] Vingert | MASCOTTE: a research facility for high pressure combustion of cryogenic propellants[END_REF], [START_REF] Vingert | Sélection de cas tests Mascotte : -cas subcritique LOX/GH2 10 bar.cas supercritique LOX/GH2 60 bar[END_REF], [START_REF] Murrone | Atomization of cryogenic rocket engines coaxial injectors[END_REF] test facility on the 10-bar operating point corresponding to cryogenic rocket engines under transient operating conditions as depicted in Figure 2.

In this paper, we focus on the improvement of the diffuse interface model for the simulation of separated two-phase flow in coaxial cryogenic injector in order to realize the same coupling strategy. Instead of using the 4 equation model, we have selected the 5 (two temperatures) equation model which can be obtained by velocity and pressure relaxation of the 7 equation model [START_REF] Baer | A twophase mixture theory for the deflagration-todetonation transition (DDT) in reactive granular materials[END_REF]. The two temperatures of the 5 equation model [START_REF] Murrone | A five equation reduced model for compressible two-phase flow problems[END_REF] provide an accurate description of the liquid temperature and is able to avoid spurious pressure oscillations of the 4 equation model due to mixing of hot gas with cold liquid in the diffuse interface.

The derivation of the 5 equation model is presented in the next section 2. Moreover, we give some details on the HLLC schemes [START_REF] Boucher | Développement d'un modèle pour la simulation de l'atomisation dans les chambres de combustion des moteurs fusées à ergols liquides -Rapport d'activités[END_REF] we have implemented in the code. We also present briefly the coupling strategy together and the atomization model we have used. Finally in section 3, we present the numerical results obtained on validation test cases and on a coaxial injector configuration in subcritical conditions. A brief conclusion is addressed at the end of the paper. Finally, we draw the principal perspectives.

MODELS AND METHODS

For the simulation of liquid gas interface, we point out that diffuse interface modelling are gaining importance [START_REF] Saurel | Simple and efficient relaxation method for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF], [START_REF] Lemartelot | Liquid and liquid-gas flows at all speeds[END_REF]. The coupling is necessary because we want to deal with transition from dense to dilute concentration of the liquid phase and interface and kinetic models present fundamental differences. More specifically, in the dilute limit, diffuse interface models involve acoustic waves propagation which are not correct. Moreover, in kinetic models, the volume occupied by liquid is neglected and particles are considered as incompressible. In the framework of diffuse interface modelling, different levels of description are available. Modelling of two-phase flows is typically based on averaging procedures [START_REF] Ishii | Thermo-fluid dynamics of twophase flow. 2 nd Edition[END_REF]. In their most general form, these averaging techniques produce models characterized by two different velocities and pressures for each phase supplemented by one or several topological equations. This is namely the 7 equation model [START_REF] Baer | A twophase mixture theory for the deflagration-todetonation transition (DDT) in reactive granular materials[END_REF], [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF]. As usual, a delicate balance between the complexity of the model and its performance has to be found. In the dense region of the two-phase flow, we propose to use a 5 equation model providing two different temperatures, one for the liquid and another one for the gas.

The 7 equation model

If we extend this 7 equation model to the case of multi-species fluid, mass conservation equation may be readily replaced by mass fractions equations and the convective part of the (5+n 1 +n 2 ) equations system with relaxation pressures, velocities and temperatures source terms can be written under the form below. Moreover, the two terms = -, = -in the system stand respectively for the heat and mass transfer.

+ = + ⨂ + ∇ = ∇ + - + ! " + " + = ! • ∇ + ! • - + $ - + + %" + ! 2 ' + = + ⨂ + ∇ = -∇ - - -! " + " + = -! • ∇ -! • - -$ - --%" + ! 2 ' + ! • ∇ = $ - + ) + * Eq. 1
The notations are classical. Firstly, / are the volume fractions of each phase ( + = 1), / the phase densities, / the vector velocities, / the pressures and " / = 0 / + / /2 the specific total energies, with 0 / the specific internal energies. On the other hand, ! , stand for the interfacial velocity and pressure for which we have chosen the following expressions:

! = ∑ / / / /3 ∑ / / /3 , = 4 / / /3
Eq. 2

The 5 equation model

The general 7 equation model contains relaxation parameters and $ which determine the rates at which velocities and pressures of the two phases reach equilibrium. Here we are interested in situations where relaxation times are small compared with the other physical characteristic times. From an asymptotic analysis, one can derive a five equation reduced model with a single velocity and pressure. The original formulation includes an equation for the volume fraction and one for the total mixture energy. This (3+n 1 +n 2 ) equation system [START_REF] Murrone | A five equation reduced model for compressible two-phase flow problems[END_REF] extended here to the multispecies case will be referred to as the ''(e,α) formulation''. The system writes:

+ = 0 + = 0 + ⨂ + ∇ = 0 " + " + = 0 + ⋅ ∇ = 7 -7 7 + 7 Eq. 3
Now, we establish a formulation of this reduced model using two equations of transport for the internal energies instead of an equation for the volume fraction and another for the total mixture energy. If the different formulations are equivalent for smooth solutions, we expect the second one to be more adapted for non-conservative term. This second formulation of the 5 equation system will be referred to as the ''(ε 1 ,ε 2 ) formulation''. Thus we set = 9 /0 and $ = $ 9 /0 where ɛ tends to zero and we look for limit equations of the 7 equation model. In this way, we use an asymptotic expansion in terms of ɛ and we try to establish the governing equations when 0 → 0. This analysis can be performed directly on the conservative form of the system. However it is more convenient to work with the set of variables / / , / , 0 / , / , / and to use the quasi-linear form of the equations as in the following. Using the notation ; / /; = ∂ = + / • ∇ for the material derivative, from the momentum and the mass conservation equations, it is easily seen that the velocities / obey: / / ; / / ; + ∇ / / = ∇ / + /> -/ Eq. 4

From these last equations, we deduce the ones for the kinetic energy / /2 of each phase:

/ / ; / / /2 ; + / • ∇ / / = / • ∇ / + / • / @ -/ Eq. 5

Then from the total energy equations, we get the equations for the specific internal energies:

/ / ; / 0 / ; + / / / = ! -/ • ∇ / + $ / @ -/ + ! -/ • / @ -/ Eq. 6

Now, we suppose a binary law state 0 / = 0 / / , / and introduce the coefficients C / , ) / for the partial derivatives:

C / = 0 / / D E F , ) / = 0 / / D G F Eq. 7
After some calculations, the transport equations for the pressures write: / ; / / ; + / I / / = I / ! -/ • ∇ / + $I / / @ -/ + ) / / ! -/ • / @ -/ Eq. 8

HLLC schemes for the 5 equation model

In this section, we propose to resolve the model in a Finite Volume framework. A lot of works have been recently dedicated to liquid-gas interface problems. These types of two-phase flows exhibit strong gradients of variables and large variation of sound speed near the interface between the phases. In this context, it is mandatory to derive accurate but also very robust schemes. The class of HLLC schemes based on the approximated resolution of the non-linear Riemann problem seems to be a good candidate. For the 5 equation model, we have presented two formulations of the system. The first one is the ''(e,α) formulation'' and the one second is the ''(ε 1 ,ε 2 ) formulation''. The last one is very interesting and gives us the possibility to examine the most adapted closure of non-conservative terms satisfying total energy conservation. These shock computational difficulties due to non-conservative character of the model has to be examined carefully. This could help us to maintain positivity for crucial variables such as density, pressure or volume fraction. Then an extension of the multi-slope technique for this two-phase flow model and their specific variables has been implemented. When we are interested in two-phase flows, we have to deal with high density ratios, strong gradients and also discontinuous solutions. Nowadays, the MUSCL technique remains a good compromise between accuracy and robustness. For this reason, we have chosen the multi-slope approach for general unstructured meshes to cope with strong discontinuities of the flow.

N = 4 -O ⋅ P QRSTUK 3 -V W N Eq. 19
The conservative part of the flux read respectively for the ''(e,α)" and the ''(ε 1 ,ε 2 ) formulations":

O ⋅ P = X Y Y Y Z [ Q [ Q [ Q + \ " + [ Q [ Q ] _ , O ⋅ P = X Y Y Y Z [ Q [ Q [ Q + \ 0 [ Q 0 [ Q ] _ Eq. 20
On the other hand, the non-conservative terms are given by the following expressions:

W = X Y Y Y Y Z 0 0 0 0 - I I + I ] _ , W = X Y Y Y Y Y Z 0 0 0 I I + I I I + I ] _ Eq. 21
If we consider the vector K constant in the cell, we can proceed to integration;

N = 4 -O + W[ Q, QRSTUK 3 ⋅ P Eq. 22
Then the conservative part of the flux at the interface is evaluated following the HLLC scheme described below. First, we consider the ''(e,α) formulation'' of the model and the only one non conservative equation for the volume fraction writes:

N = 4 -g [ Q t + W[ u Q, hP QRSTUK 3
Eq. 25

Then, according to the sign of the waves P ` , P a , P b , we are able to compute:

[ Q t = ,`[Q,`, , * P b , ,a * P b , ,a [ a , [ u Q, = [ Q,`, P b , P b , [ a Eq. 26
Then combining, the two components of the mass fractions (Eq. 23) in the star region, the volume fractions read:

, * = v 1 + v `, v `= ` ` , * , * , ,a * = v a 1 + v a , v `= a a ,a * ,a * Eq. 27
Then, let us suppose isentropic rarefaction across the waves P ` , P a and we can deduce the real phasic density: /, * -/,`= * -7 u /,` , /,a * -/,a = a *a 7 u /,a Eq. 28

Now, we consider the second ''(ε 1 ,ε 2 ) formulation'' of the 5 equation model. We have to deal with the non-conservative equations for the internal energy of each phase:

N / / 0 / = 4 g-/ / 0 / [ Q w + W[ u Q, hP Eq. 29 QRSTUK 3
Then, according to the sign of the waves P ` , P a , P b , we are able to compute:

/ / 0 / [ Q w = / / `0/,`[`,
/ / * 0 /, * P b , / / a * 0 /,a * P b , / / a 0 /,a [ a , [ u Q, = [ Q,`, P b , P b , [ a Eq. 30

Now we have to determine the values of internal energies in the star region 0 /, * , 0 /,a * which are the only unknowns. We have implemented and tested different solutions presented below.

Coupling strategy and atomization model

The classical source terms between the diffuse interface model and the kinetic model are the drag force (velocity relaxation), the heat transfer (temperature relaxation) and the mass transfer. n • , ' ƒ , " …SE Eq. 32

The correlation of Schiller-Naumann is used for the drag force and the Abramzon-Sirignano model is used for the heat and mass transfer. Now we describe the source terms s(u) for the spray solver and then focus on the primary atomization. The first component traduces the growth of mass due to primary atomization and its decrease by vaporization. If we now look at the source term S(U) affected to the carrier phase. The first ng components of gas species concern combustion and evaporation phenomena. The first term includes the evaporation of liquid oxygen droplets in the "dispersed" phase which creates gas oxygen in the carrier phase. The component number ng +1 standing for the transport equation of mass liquid fraction concerns primary atomization. It transforms the dense LOx into "dispersed" LOx.

r [ = † r G ‡ r G ‡ … ‡ r G ‡ ˆ ‡ r ‰ ‡ Š = X Y Y Y Z " S}{ -‹ ~ " …SE r G ‡ ~+ ‹ ~n• r G ‡ ℎ ~+ ‹ ~ n • • ~+ ' ƒ ‹ " S}{ ] _ P • = X Y Y Y Y Y Y r GŽ •• r GŽ ' ⋮ r GŽ "" r GŽ • P … P U ] ^= X Y Y Y Y Y Y Y " -+ ‹ ~ " …SE " ⋮ " Q- -" S}{ -r G ‡ … ‡ -r G ‡ ˆ ‡ ] ^
Eq. 33

The model used to describe the mass transfer between solvers accounting for primary atomization reads in:

" S}{ = ˜ S}{ S}{ ˜ Eq. 34
where ˜ is the liquid mass in a given volume of control, S}{ is the characteristic frequency of the primary atomization process, and S}{ ˜ is an efficiency function. We assume the atomization frequency to be directly connected to the strength of the velocity gradient, which is the only information locally available in the 5 equation model. This could be estimated using several approaches, amongst which the Q criterion, the vorticity or the resolved strain tensor, all being based on the velocity gradient. In this study we have chosen to use the latter approach: }™|š = ›2; ; = 4 ; oe ; oe oe

; oe = 1 2 % [ • ž oe + [ • oe ž ' Eq. 35
The efficiency function reads in [START_REF] Murrone | Atomization of cryogenic rocket engines coaxial injectors[END_REF]. It is designed to make sure that when some LOx mass is transferred from the fluid towards the spray in a given volume of control the corresponding vanishing volume in the fluid is actually negligible.

Otherwise the gas would experience some unphysical expansion in the volume control, which obviously has to be avoided, and the dispersed hypothesis made for the spray would not be respected. In other words, we use the numerical diffusion which spreads the interface over several mesh elements in order to carry out the mass transfer in a smooth way. As this point with this 5 equation model, except for the temperature, the properties of the created droplets resulting from the primary atomization have to be presumed. Actually, properties are estimated based on the instability analysis from Villermaux [START_REF] Marmottant | On spray formation[END_REF]. In the latter work, the drop size and velocity distributions of the spray are estimated as a function of the injected propellant properties (density ratio, inlet velocities, vorticity thickness...). Consequently, the knowledge of the steady operating conditions of the MASCOTTE configuration under study permits to derive an overall mean droplet diameter subsequent to the primary atomization process and a corresponding mean droplet velocity. The direction given to the droplet velocity in each mesh cell has been set to that of the fluid, which may be actually a rough approximation. Finally, the temperature of the created droplet is the LOx dense temperature of the interface diffuse model. This is the first advantage of the 5 equation model.

The second one is the correct description of the thermodynamic mixture between liquid and gas.

Therefore the spurious pressure and temperature oscillations encountered in the 4 equation model at the interfaces are now prevented.

RESULTS

Water-air shock tube

The length of the domain is 1 m and initially the interface is located at x = 0.7 m. The tube is initially filled with high pressure liquid water on the left side and with air on the right side. This test problem consists of a classical shock tube with two fluids and admits an exact solution. The initial condition consists in a pressure discontinuity between p = 10 9 Pa in the liquid side and p = 10 5 Pa in the gas side. The right and left chambers contain nearly pure fluids: the volume fraction of the gas in the water chamber is 10 -8 and inversely the water volume fraction is 10 -8 in the gas chamber. Figure 3 displays for the different numerical methods the volume fraction, the mixture density, the pressure and the velocity. The exact solutions are represented on these curves. The results are shown at time 229µs and seem to be of comparable accuracy with respect to the exact solution. Figure 3 plots the results at the order 1 in space. For the monophasic variables such as the temperatures, we emphasize that different schemes could give very different results and don't have the same robustness. Nevertheless, the formulation with two internal energies seems very promising compared to the one with the volume fraction and total energy transport equations. For the simulations, a Stiffened Gas equation of state [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow model[END_REF] has been implemented for each phase. It takes into account for attractive and repulsive molecular effects.

/ = Ÿ / -1 / 0 / -/ -Ÿ / ¡ / / = Ÿ / -1 z …/ / -¡ / 0 / = z …/ / + ¡ / / Eq. 36
The Stiffened Gas law is used for both gas and liquid. Parameters are summarized below: In order to reduce the mesh size and solution time, the simulation is performed with a 60° sector geometry. The geometry and the mesh are presented in figure 4. As seen in the Figure 4, the mesh is of course refined near the injection where the atomization takes place. Far from validation with experiments results, we want to show the enhancement of the coupling strategy with the 5 equation interface diffuse model. The law state for the LOx taking into account for liquid compressibility which has been used in the coaxial cryogenic injector test case reads:

, = 1 = 1 9 1 + £ 9 ¤-¤ 9 1 + 9 -9 ¥ . 38
The different parameters for the oxygen are summarized below: The capacity of the coupling strategy is illustrated in the following Figures. In Figure 5, we present an instantaneous visualization of the volume fraction of the liquid jet (iso-surface) as well as the particles generated by atomization. The physical time simulated corresponds to 74ms.

In Figure 6, we present two instantaneous isovalues field of the volume fraction at times 65ms and 74ms. This time interval corresponds to twice the convective characteristic convection time of the liquid core. We compare results obtained with atomization (upper side of the axis) and without atomization. (Color map legend: min = dark, max = white, with cutoff below = 0:01).

Moreover, quantitative result such as length of the liquid core is plotted. In the Figure 7, we present the liquid core representation thanks to iso-values of liquid volume fraction between [0.95±0.04]. The length of the liquid core renormalized by the injector diameter is plotted over time for two times the characteristic convection time (4,5ms) of the liquid core. We can observe instabilities along the axial direction. Comparison is made for the simulations with and without atomization. The length of the liquid core is also very similar when taking into account or not for atomization. This confirms that the atomization procedure has been carefully implemented. The results seem very promising.

CONCLUSION

We have successfully implemented a 5 equation interface diffuse model for the simulation of the dense "separated two-phase flow" within a coaxial LRE injector in sub-critical conditions. This new level of modelling allows us to perform accurate and robust simulations with a correct description of thermal transfer in the interface between liquid and gas. Moreover, the 5 equation model has been coupled with a spray kinetic solver following the fully Eulerian strategy proposed in our previous works. So, the atomization process has been taken into account and we have performed a difficult simulation showing very promising results.

Concerning on going and future works, we intend to develop this Eulerian coupling strategy with a hierarchy of interface diffuse models and kineticbased moment method (KBMM).

Figure 2 :

 2 Figure 2: LES of a subcritical combustion LOx/GH2 (10-bar operating MASCOTTE 10 facility) (Liquid oxygen in blue, hydrogen in red and flame in yellow)

  0 /, * = 0 /,` x xs 7 0 /, * = 0 /, * , * x xs y /, * = /,`, 0 /, * = 0 /, * , * x xs z 0

  4.4, z …/ = 4000, ¡ / = 6. 10 ¢ Eq. 37

Figure 3 :

 3 Figure 3: results of the 5 equation model with different HLLC schemes on the water air shock tube

  1.71 10 §ª z … = 951«/© /• Eq. 39On the other hand a classical "perfect" Stiffened gas law state is used for the hydrogen GH 2 :Ÿ = 1.4 ¡ = 0 z … = 10112 «/© /• Eq.[START_REF] Vallet | Contribution à la Modélisation de l'Atomisation d'un Jet Liquide Haute Pression[END_REF]

Figure 4 :

 4 Figure 4: geometry and mesh of the coaxial injector

Figure 5 :Figure 7 :

 57 Figure 5: instantaneous visualization of the jet with droplets produced by atomization

In these last equations, we have introduced the phasic sound speeds 7 / and the acoustic impedances I / . Then, the notations 7 / , I / stand for the same quantities evaluated at the interfaces:

I / = / 7 / , I / = / 7 / Eq. 9

Now, we perform the asymptotic analysis introducing the following expansion in term of 0 for velocities and pressures but also for the other variables:

/ = 9 + 0 / , / = 9 + 0 / Eq. 10

The mass conservation equations at order 0 read:

The equations for velocities at order 0 read:

= 9 ∇ / 9 + 9 /> -/ Eq. 12

Then we can deduce at the order 0 the equations for internal energies and pressures: Then, combining the above equations, we can get the pressures fluctuations at the order 1:

$ 9 -= 9 9 I 9 -I 9 9 I 9 + 9 I 9 9 Eq. 15

Now, if we use the equations at order 0 for internal energies, and we introduce the expression of the pressure fluctuations, we get: Eq. 16

Now using the mass conservation, internal energy equations read: Eq. 17

Finally, the 5 equation model referred to as the ''(ε 1 ,ε 2 ) formulation'' can be written under the form below. This new formulation will be used for the design of numerical fluxes.


Dr, Senior scientist, Fabrice Poirion

Quentin Mercier

A multiobjective optimization based approach for RBDO

In this paper we present a novel algorithm in order to solve multiobjective design optimization problems of a sandwich plate when the objective functions are not smooth and when uncertainty is introduced into the material properties. The algorithm is based on the existence of a common descent vector for each sample of the random objective functions and on an extension of the stochastic gradient algorithm. It will be shown that a chance constraint optimization problem such as a RBDO problem can be written as a multiobjective optimization problem. Chance constraint optimization problems yields optimal designs for a fixed given level of probability for the constraint. However in real life problem it is not realistic to introduce a given probability because it is not known. It is more efficient to solve the problem for a whole range of probability in order to obtain an overview of the probability level appearing in the constraint effect on the solution. We show in this paper how to transform a chance constraint optimization problem into a multiobjective optimization problem and we give an illustration on simple examples.

Manufacturers are ever looking for designing products with better performance, higher reliability at lower cost and risk. One way to address these antagonistic objectives is to use multiobjective optimization approaches. But real world problems are rarely described through a collection of fixed parameters and uncertainty has to be taken into account, may it appear in the system description itself or in the environment and operational conditions. Indeed the system behavior can be very sensitive to modifications in some parameters [START_REF] Papadimitriou | Effects of structural uncertainties on tmd design : A reliability-based approach[END_REF]; [START_REF] Matthies | Uncertainties in probabilistic numerical analysis of structures and solids-stochastic finite elements[END_REF]; Arnaud and Poirion (2014). This is why uncertainty has to be introduced in the design process from the start. Optimization under uncertainty has known important advances since the second-half of the 20th century [START_REF] Dantzig | Linear programming under uncertainty[END_REF]; [START_REF] Bellman | Decision-making in a fuzzy environment[END_REF] and various approaches have been proposed including robust optimization, where only the bounds of the uncertain parameters are used, and stochastic optimization where uncertain parameters are modeled through random variables with a given distribution and where the probabilistic information is directly introduced in the numerical approaches. In that context the uncertain multi objective problems is written in terms of the expectation of each objective. In our paper we shall focus on this last interpretation of the optimization problem.Considering single objective stochastic optimization problems, a large variety of numerical approaches [START_REF] Sahinidis | Optimization under uncertainty: State-of-the-art and opportunities[END_REF]; [START_REF] Roy | Recent advances in engineering design optimisation : challenges and future trends[END_REF] can be found in the literature. Two main distinct approaches exist, one based on stochastic approximations such as the Robbins Monro algorithm and the various stochastic gradient approaches [START_REF] Robbins | A stochastic approximation method[END_REF]; [START_REF] Ermoliev | Stochastic quasigradient methods and their application to systems optimization[END_REF]; [START_REF] Ermoliev | Numerical Techniques for Stochastic Optimization[END_REF], the second one based on scenario approaches [START_REF] Shapiro | Handbooks in OR & MS[END_REF]; [START_REF] Nemirovski | Scenario approximations of chance constraints[END_REF], the latter being more frequently applied for chance con- strained problems. Again two directions can be found, a robust approach and a scenario based approach used to calculate an estimate of the mean objective function [START_REF] Fliege | Stochastic multiobjective optimization: Sample average approximation and applications[END_REF]; [START_REF] Bonnel | Stochastic optimization over a Pareto set associated with a stochastic multi-objective optimization problem[END_REF]; [START_REF] Mattson | Pareto frontier based concept selection under uncertainty, with visualization[END_REF]. RBDO and more generally chance constraint problems are numerically difficult to solve and moreover their solution is obtained for a single given of chance level. An interesting situation would be to construct the solutions for a whole range of probability levels. Writing the probabilistic constraint as the expectation of particular random function, we show that if the original chance constraint problem is replaced by a stochastic multiobjective optimization problems, the Pareto solution set of the new problem contains the solutions for all levels of probability.

1. OPTIMIZATION OF UNCERTAIN OB-JECTIVES Let (Ω, A , P) be an abstract probabilistic space, and W : Ω → R d a random vector. We denote µ the distribution of the random variable W and W its image space W (Ω). Let W 1 , ...,W k , ... independent copies of the random variable W which will be used to generate independent random samples with distribution µ. Consider m convex functions f i : R n × W → R, i = 1, ...m depending on uncertain parameters modeled trough random vector W (ω). In this paper we shall consider the following optimization problem :

min x∈R n {E[ f 1 (x,W (ω))], ..., E[ f m (x,W (ω))]} . (1)
More precisely we want to construct the associated Pareto set: multiobjective optimization is based on the notion of Pareto optimal and weak Pareto optimal solutions. Consider m convex functions f i : R n → R, i = 1, ...m and the unconstrained optimization problem

min x∈R n { f 1 (x), ..., f m (x)} .
(2)

A solution x * of problem (2) is Pareto optimal if no point x such that f i (x) ≤ f i (x * ) ∀i = 1, ..., m and f j (x) < f j (x * ) for an index j ∈ {1, ..., m} exists. It is weakly Pareto optimal if no point x such that f i (x) < f i (x * ) ∀i = 1, ..., m exists. A complete review on multiobjective optimization can be found in [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. Before going on with the algorithm description that will be used to solve the previous problem we shall recall definitions of some notions appearing in the context of nonsmooth analysis and multiobjective optimization. Throughout the paper the standard inner product on R n will be used and denoted •, • , the norm being denoted • .

COMMON DESCENT DIRECTION

The algorithm presented in the next section is based on the existence and construction of a descent direction. We first recall its definition.

Definition 1 A vector d is called a descent direc- tion if ∃t 0 > 0 such that f (x + td) < f (x) for all t ∈ [0,t 0 ].
For smooth functions it is well known that the opposite direction of the gradient is a descent vector.

In the nonsmooth convex or nonconvex context not all elements of the subdifferential is a descent vector. There exist several techniques to construct such a descent vector: proximal bundle methods [START_REF] Kiwiel | Methods of Descent for Nondifferentiable Optimization[END_REF]; [START_REF] Wilppu | New multiple subgradient descent bundle method for nonsmooth multiobjective optimization[END_REF]; [START_REF] Mäkelä | Mathematical Modeling and Optimization of Complex Structures[END_REF], quasisecant methods Bagirov et al. ( 2013), or gradient sampling methods [START_REF] Burke | Approximating subdifferentials by random sampling of gradients[END_REF][START_REF] Burke | A robust gradient sampling algorithm for nonsmooth, nonconvex optimization[END_REF].

Considering now m functions f 1 , ..., f m we show that there exists a vector d which is a descent direction for each function. Its construction is based on properties of the following convex set C :

Lemma 1 Let C be the convex hull of either 1. the gradients ∇ f i (x) of the objective functions when they are differentiable, 2. or the union of the subdifferentials The existence of the common direction d and its construction is given by the next theorem:

∂ f i (x), i = 1, ...,
Theorem 1 Let C be the convex set defined in Lemma 1 and p * its minimum norm element. Then either we have 1. p * = 0 and the point x is Pareto stationary or 2. p * = 0 and the vector -p * is a common descent direction for every objective function.

We have now the sufficient materials to present the SMSGDA (Stochastic Multi Descent Algorithm) algorithm.

THE SMGDA ALGORITHM

As written problem (1) is a deterministic problem but the objective function expectations are seldom known. A classical approach, the sample average approximation (SAA) method, is to replace each expectancy by an estimator built using independent samples w k of the random variable W , Bonnel and Collonge (2014); [START_REF] Fliege | Stochastic multiobjective optimization: Sample average approximation and applications[END_REF]. The algorithm we propose does not need to calculate the objective function expectancy and is based only on the construction of a common descent vector. Let ω be given in Ω and consider the deterministic multiobjective optimization problem:

min x∈R n { f 1 (x,W (ω)), f 2 (x,W (ω)), ..., f m (x,W (ω))]}
(3) Following theorem 1 there exists a descent vector common to each objective function f k (x,W (ω)), k = 1, ..., m at point x.

The common descent vector depends on x and ω and therefore will be considered as a random vector denoted d(ω) defined on the probability space (Ω, A , P).

The algorithm

We give now the successive steps of the algorithm that we propose.

1. Choose an initial point x 0 in the design space, a number N of iterations and a σ -sequence t k :

∑t k = ∞ ; ∑t 2 k < ∞,
2. at each step k draw a sample w k of the random variable W k (ω), 3. construct the common descent vector d(w k ) using theorem 1 and the gradient sampling approximation method, 4. update the current point :

x k = x k-1 +t k d(w k ).
The last step of the algorithm defines a sequence of random variables on the probability space (Ω, A , P) through the relation

X k (ω) = X k-1 (ω) -t k d(X k-1 (ω),W k (ω)) (4)
Theorem 2 [START_REF] Mercier | A stochastic multi gradient descent algorithm[END_REF] ) Under a set of assumptions, 1. the sequence of random variables X k (ω) defined by relation (4) converges in mean square towards a point X * of the Pareto set:

lim k→+∞ E[ X k (ω) -X 2 ] = 0.
2. The sequence converges almost surely towards X * .

P ω ∈ Ω, lim k→∞ X k (ω) = X = 1.

SOLVING RELIABILITY PROBLEM USING A STOCHASTIC MULTIOBJECTIVE FORMULA-

TION

Introducing probabilistic constraints is a rather natural way to take into account the notion of risk in an optimization process. Let us consider the following chance constraint problem:

argmin x∈X ad {E[ f (x, ξ (ω))] | P[g(x, ξ (ω)) ≥ 0] ≥ p 0 }.
(5) Here, X ad is a feasible closed convex set of the set of control variables X, g : X × R → R represents a physical or structural quantity. In this formulation failure occurs when g(x, ξ (ω)) becomes positive, and p 0 denotes the level of risk one is ready to accept. Such a problem is rather difficult to solve. The reason is twofold: first it is very difficult to check whether a given chance constraint is satisfied at a given point x or not. Typically Monte-Carlo simulation is the only way to estimate the probability of violating the constraint, but becomes too costly when p 0 approaches unity. The second reason comes from the fact that the feasible set of problem ( 5) can be nonconvex even if the set X ad is convex as well as function g. Several developments can be found in the literature in order to overcome those difficulties: transforming the problem into a combinatorial problem by discretizing the probability distribution [START_REF] Dentcheva | Concavity and efficient points of discrete distributions in probabilistic programming[END_REF], using convex approximation [START_REF] Nemirovski | Scenario approximations of chance constraints[END_REF] or sample average approximations [START_REF] Luedtke | A sample approximation approach for optimization with probabilistic constraints[END_REF]; [START_REF] Pagnoncelli | Sample average approximation method for chance constrained programming: theory and applications[END_REF].

The general formulation of a reliability based design optimization (RBDO) problem is the following

min x∈R n { E[ f (x,W (ω))]} s.t. P{g(x,W (ω)) ≥ 0} ≤ α. (6)
P{g(x,W (ω)) ≥ 0} represents the probability of failure of the design x, W (ω) is a random variable modeling the uncertainties and α represents the threshold of failure probability authorized. In most cases, α is a very small positive number. Let us remark that P{g(x,W (ω)) ≥ 0} =

E[I R + (g(x,W (ω)))] def = E[G(x,W (ω))]
, where I R + denotes the indicator function of R + : I R + (x) = 1 when x is positive and is equal to 0 otherwise. We replace the RBDO problem (6) by the following stochastic multiobjective optimization problem

min x∈R n { E[ f (x,W (ω))], E[G(x,W (ω))]} . (7)
It is clear from the definition of the Pareto front that for a given value of α, the point (x α , α), where x α is a solution of the RBDO problem for the given threshold, lies on the Pareto front of the problem (7). Therefore solving problem (7), one obtains directly the set of solutions of the RBDO problem ( 6) for all values of probability of failure α ∈ [0, 1]. However the SMGDA algorithm cannot be used directly to solve this last problem since the second objective function is not locally Lipschitz. We use the mollifier introduce in (Andrieu et al., 2011) in order to render this objective smooth. More precisely we introduce a smooth non-negative even function h r (x) = 1 r h(x/r) such that it reaches its maximum for x = 0 and 1 r

+∞ -∞ h x r dx = 1
We have then:

I r R + (z) ≈ 1 r +∞ 0 h( z -y r )dy,
which yields the following expression :

I r R + (g(x,W (ω))) = 1 r +∞ 0 h g(x,W (ω)) -y r dy.
This approximation can be differentiate with respect to x and it can be checked that :

∇I r R + (g(x,W (ω))) = - 1 r h g(x,W (ω)) r ∇g(x,W (ω)).

ILLUSTRATION

The above development is applied to a simple example used as a reference test cas in several publications [START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive Kriging surrogate models[END_REF]. The problem is to minimize the cross-sectional area b × h of a rectangular column submitted to a compression load F while avoiding buckling, which occurs when the compression load is higher than the critical Euler force. Failure will occur when:

F - π 2 Ebh 3 12L 2 ≥ 0, ( 8
)
where E is the Young's modulus of the column material and L its length. In their paper the authors considered the parameters E and L uncertain and modeled as lognormally distributed random variables. The critical Euler force is also considered uncertain through the appearance of a multiplicative lognormally distributed random variable k in its expression. The RBDO problem is then written

Argmin h,b b × h (9)
under the constraints 2016) that an analytic solution of the RBDO problem ( 9) and ( 10) exists and is given by:

h -b ≤ 0 and P F - k(ω)π 2 E(ω)bh 3 12L(ω) 2 ≥ 0 ≤ α (
b * = h * = 12F π 2 exp(λ k +λ E -2λ L +Φ -1 (α) √ ζ 2 k +ζ 2 E +4ζ 2 L 1 4 , (11)
where ζ . = ln(1 + δ 2 . ) and λ . = ln(µ . )-1/2ζ 2 . , µ and δ being the mean and covariance of the lognormal random variables, and Φ the cumulative distribution function of the standard normal distribution.

We consider now the stochastic multiobjective optimization problem min

b,h {b × h, E[G(b, h,W (ω))]} ; h -b ≤ 0, (12) where W (ω) = (k(ω), E(ω), L(ω)) and G(b, h,W (ω)) = ĨR + (F -k(ω)π 2 E(ω)bh 3 12L(ω) 2
). The SMGDA algorithm is used in order to construct the Pareto front. At the same time the analytic solutions (h * , b * ) of the RBDO problem (9) are constructed for a set of values of α from which the analytic Pareto front (b * × h * , α) is obtained. In their paper [START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive Kriging surrogate models[END_REF] used an adaptative Kriging surrogate model in order to solve the RBDO problem for the specific values α = .05. Figure 1 shows the comparison of the two Pareto fronts, the solution found by [START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive Kriging surrogate models[END_REF] being represented by a green triangle. Almost every solution proposed by SMGDA stick to the analytic Pareto front, and thus can be consider as good. During the optimization process SMGDA does not require the costly calculation of the failure probability and is still able to converge, even for very small values. In this way we obtain the solutions of the original RBDO problem for all probability (or failure) levels, allowing to evaluate the impact of this level on the design parameter solutions. The SMGDA algorithm used is completely parallelizable and the numerical problem of evaluating the probability constraint is decoupled from the optimization procedure. It becomes a postprocessing procedure where any adequate method can be used. research, 212, 345-351. Arnaud, R. and Poirion, F. (2014). "Optimization of an uncertain aeroelastic system using stochastic gradient approaches." Journal of Aircraft, 51(3), 1061-1069.

Bagirov, A., Jin, L., Karmitsa, N., Nuimat, A. A., and Sultanova, N. (2013). "Subgradient method for nonconvex nonsmooth optimization." Journal of Optimization Theory and Applications, 157, 416-435.

Figure 1 :

 1 Figure 1: Pareto front of the reliability problem

  .,Cohen, G., and Vàzquez-Abad, F. (2011). "Gradient-based simulation optimization under probability constraints." European Journal of Operational

  13th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13 Seoul, South Korea, May 26-30, 2019

  International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13 Seoul, South Korea, May 26-30, 2019

	m when they are nondifferentiable but
	convex or
	3. the union of the Clarke's subdifferentials
	∂ f

i (x), i = 1, ..., m if they are nonconvex. Then there exists a unique vector p * = Argmin p∈C ||p|| such that ∀p ∈ C : p T p * ≥ p * T p * = ||p * || 2 .

13th

Table 1 :

 1 10)The chosen parameters for the distributions are recalled in table 1. The value of the compression force is chosen as F = 1.4622 × 10 6 N. With 13th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13 Seoul, South Korea, May 26-30, 2019 Probabilistic model for the column problem

	parameter distribution mean :µ covariance : δ
	k	lognormal	0.6	10
	E (MPA)	lognormal	10 4	5
	L (mm)	lognormal	3	1
	these specific probability distributions it is shown in
	Moustapha et al. (

13th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13 Seoul, South Korea, May 26-30, 2019