Eulerian coupling of two-phase flow models for the large eddy simulation of the atomization in cryogenic combustion chamber

A Large Eddy Simulation (LES) of the primary atomization in cryogenic combustion chamber is proposed. To reach this crucial challenge, different two-phase flow models are necessary because the simulation includes multiscale phenomena. A fully Eulerien coupling strategy has been implemented in a multiphysics platform and a primary atomization model has been suggested. This strategy combines a diffuse interface model for the atomization of the liquid jet and a kinetic based model for the combustion of the spray. Concerning numerical methods, a time splitting and a new specific Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) technique have been used for robustness and precision of the computation. The configuration of the simulation is a test cryogenic facility corresponding to one coaxial injector.


INTRODUCTION
To ensure the reliability of launcher propulsion systems, phenomena occurring in rocket engines have to be studied in a more comprehensive way. Concerning liquid propulsion, high-frequency instabilities have been encountered in the recent past for cryogenic rocket engines combustion chambers. Such instabilities which can cause highly destructive damages result from the complex interaction between propellants injection, §ame dynamics, and acoustic modes. In this study, it is proposed to focus on the primary atomization which remains a key point in the understanding and modeling of cryogenic rocket engines because the physical processes in the chamber are highly dependent on the characteristics of the spray produced by the injector. As a consequence, atomization is expected to be deeply linked to the high-frequency instabilities. Even if it has been extensively studied in the last decades, the parameters describing the spray, like the expansion angle, the penetration depth, and the droplet size distribution, are usually given as input data in the numerical simulation codes of such engines.
The aim of this paper is to compute the spray produced by a coaxial injector directly from the upstream global parameters (injection pressures, mass §ow rates, etc.) instead of imposing its characteristics as an input of the simulation. This goal can be reached if we realize a LES of the primary atomization. Taking into account primary atomization seems to be possible if we consider the current computational resources. However, this kind of simulation, even only applied to one single coaxial (liquid oxygen (LOx) gaseous hydrogen (GH 2 )), remains a tremendous challenge when considering atomization of the jet combined with the combustion of the spray. Indeed, the major issue is the large variety of multiscale phenomena such as combustion, evaporation of oxygen droplets, turbulence, and atomization. Nevertheless, we want to point out that the simulations presented in this paper are very challenging and new for our knowledge. We expect that CFD (Computational Fluid Dynamics) tools combined with experiments may thus help the development and enhancement of launcher propulsion systems.
As depicted in Fig. 1 of a coaxial cryogenic injector, the strong velocity di¨erence between the two phases (LOx and GH 2 ) generates §uctuating accelerations. Due to these §uctuations, Rayleigh Taylor instabilities destabilize the liquid to create ligaments. These instabilities then grow and eventually provoke the peeling of the main LOx jet, which is referred to as ¤primary atomization.¥ Large random-shaped liquid structures are thereby ejected towards the gas §ow, Figure 1 Con¦guration of a combustion chamber within liquid-propellant rocket engines subsequently undergoing ¤secondary breakup¥ when inertia forces exceed the liquid surface tension. This results in a spray of smaller LOx droplets, mainly, spherical, which are dispersed by the turbulent gas §ow and ¦nally vaporized to feed the combustion with hydrogen H 2 . Such a con¦guration, therefore, exhibits a two-phase §ow where the liquid phase is only composed of LOx, whereas the gas phase is made up with hydrogen, vaporized oxygen, and combustion products. Eventually, the resulting hot and high-pressure combustion products exhaust through a nozzle at supersonic speed, thereby providing the required thrust.
For the simulation, we propose an original fully Eulerian coupling strategy between a di¨use interface model for ¤separated LOxH 2 two-phase §ow¥ and a kinetic based model for the ¤LOx droplets dispersed in the §ow.¥ Among the di¨use interface models, we have selected the 4-equation locally homogeneous §ow (LHF) system which can be derived from the 7-equation general system after pressures, velocities, and temperatures relaxations. Concerning the Eulerian kinetic based model, the sectional approach instead of the sampling one has been selected because it is well adapted to fragmentation and atomization of the liquid jet. The coupling has been implemented in the ONERA£s CEDRE [1] platform for energetics which is adapted for multiphysics coupling solvers. The necessity of the coupling is due to the topology of the two-phase §ow. Following the classi¦cation proposed by Ishii [2], this kind of simulation exhibit a dense ¤separated¥ liquid phase near the injector as well as a ¤dispersed¥ liquid phase composed by oxygen droplets resulting from primary atomization of the jet and secondary atomization (fragmentation of breakup) of ligaments. Because mathematical models are di¨erent to deal with ¤separated¥ or ¤dis-persed¥ topologies, a coupling of solvers has been proposed. Note that other approaches [3,4] for propulsion and also for internal combustion engine have been proposed with, for example, Lagrangian description of the spray or Level Set and Volume of Fluid description of the dense liquid §ow near the injector.
In the ¦eld of modeling primary atomization, the two principal applications are cryogenic rocket engines and combustion in Direct Injection Diesel. There are a lot of works only dedicated to primary atomization and using tracking interface methods. Here, we focus on studies including together the simulation of the atomization of the jet and the combustion of the spray. Concerning cryogenic combustion, the works [3,5,6] initially started for Diesel engines are very interesting. A 4-equation di¨use interface model (LHF) is combined with a sophisticated transport equation for the interfacial area. This last equation contains source terms for creation and destruction of interface area which are closed in the context of turbulent Reynolds-averaged NavierStokes (RANS) approach. Then, the basis of the model has been used in [7,8] and the global strategy has been used for a coaxial injector in [9]. More recently, the authors of [10,11] proposed a new closure of the source terms for interfacial area. Concerning Diesel Injec-tion, a new model is proposed in [12] and combines the LHF-RANS approach with a Lagrangian description for the spray. Following this way, the authors of [4] make the source terms depending on ¤dispersed¥ or ¤separated¥ topology of the two-phase §ows. The strategy is continuing in the RANS context with [13] but it seems that the strong coupling between Eulerian and Lagrangian methods induces some di©culties [14]. Then, the authors of [1517] focus on the transport interfacial area equation in the LES context and use also accurate tracking interface methods.
In the ¦rst approach (Eulerian LHF model with interfacial area transport equation), the principal drawback is the simpli¦ed description of the spray which does not take into account the polydispersion. In the second approach, coupling EulerianLagrangian methods sometimes associated with a tracking interface method, several di©culties arise. The ¦rst one is due to the lack of robustness of the Lagrangian methods in the case of strong two-way coupling. The second one is related to the statistical interpretation of numerous numerical particles produced in the atomization area. Finally it seems also that taking into account for the compressibility of the §uid with the tracking interface method is not an easy task. In this context, the fully Eulerian coupling method between di¨use interface model for the jet and a kinetic approach for the spray seems to be a new approach and also promising. In e¨ect, the method is fully compressible and can be applied to general physical §ows. It provides the classical advantages of Eulerian methods (implicit time integration schemes, e©cient parallel computations, robust simulations in case of strong two-way coupling). It makes possible to perform fully coupled computations with atomization of the jet and a precise description of the spray.
The plan of this work is as follows. In the next section 2, we give a brief description of the derivation of the two di¨erent two-phase §ow models. We present the derivation of the 4-equation LHF model in the hierarchy of the di¨use two-phase §ow interface models starting from the 7-equation model of BaerNunziato [18]. Then, we derive the Eulerian kinetic based model starting from the Boltzmann-like equation [19]. For the description of the ¤dispersed two-phase §ow,¥ the sectional approach has been selected with a piecewise linear reconstruction.
In section 3, the principal numerical methods are described. The time splitting technique is presented as well as the new multislope MUSCL technique [20] suitable for general unstructured meshes. This method provides a good compromise between accuracy and robustness. This method has been developed to overcome the di©culties of this kind of two-phase §ow simulations relying on the presence of strong gradients and high density ratios due to the simultaneous presence of liquid and gas.
Then, in section 4, it is described how the 4-equation LHF model and the sectional approach are combined following and improving the strategy previously proposed in [2123]. Moreover, a formulation for the atomization model with new source terms between the dense LOx and the spray composed to oxygen droplets is proposed. The fully Eulerian coupling strategy between ¤separated¥ and ¤dispersed¥ two-phase §ows solvers has been implemented in the ONERA£s CEDRE code [1].
Finally, in section 5, this tool has been applied to the numerical simulation of the MASCOTTE test facility [24] on the 10-bar operating point corresponding to cryogenic rocket engines under transient operating conditions. With this new strategy, it was possible to carry out a LES of a single element rocket engine injector taking into account the primary atomization of the jet simultaneously with the combustion of the spray. Interactions occurring between neighbor elements in actual combustion chambers, like the Vulcain 2 engine of the Ariane 5 launcher where there are several hundreds of elements on the face plate, are, therefore, not taken into account, but they may be assumed negligible for the primary atomization phenomenon of interest in this paper. In parallel to this numerical work, a lot of experimental investigations have been realized in order to provide relevant experimental data [25,26].

MODELING TWO-PHASE FLOWS
In this section, the model used for the description of ¤separated and dispersed two-phase- §ow¥ is established.

Di¨use Interface Model
Here, one of the possible relaxation methods [27], which allows one to establish the hierarchy of di¨use interface models ranging from the 7-equation model to the 4-equation model precisely used in this paper, is brie §y described. The most general interface di¨use model to describe ¤separated two-phase §ow¥ is the 7-equation model of Baer and Nunziato [18] based on averaging procedures [18]. If one extends this system to the case of multispecies §uid, the convective part of the (5 + n 1 + n 2 )-equation system with relaxation pressures, velocities, and temperatures source terms can be written under the form: In the system, the set of conservative variables is: We have supposed that the two §uids are, respectively, composed to i = 1, . . . , n 1 and j = 1, . . . , n 2 chemical species with Y i 1 and Y j 2 standing, respectively, for the mass fractions of the two §uids. The notations are classical; α k are the volume fractions of each phase; ρ k are the phase densities; u k are the vector velocities; p k are the pressures; and e k = ε k + (1/2)u k · u k are the speci¦c total energies with ε k being the speci¦c internal energies. On the other hand, p I and u I stand for the pressure and the velocity at the interface. Then, the model contains relaxation parameters λ, µ, and θ which determine the rates at which pressures, velocities, and temperatures of the two phases reach equilibrium. Here, we are interested in situations where relaxation times are small compared with the other physical characteristic times relative to the problem presented in the previous section. Thus, let us set λ = λ ′ /ε, µ = µ ′ /ε, and θ = θ ′ /ε where ǫ tends to zero. In the following of this subsection, let us brie §y describe the relaxation method process [27] which corresponds to an asymptotic analysis in the case ε → 0. Let us consider the following hyperbolic system with sti¨source terms which writes in one-dimensional space: The objective is to put the system on equilibrium and we expect to ¦nd the solution to be close to the subset: Moreover, let us suppose that a parametrization of this space ζ is explicitly known which can be written as With these assumptions, let us dM u be the Jacobian matrix and P the projection matrix on kernel ker (R ′ (M (u))).
Then, one can prove that the reduced system obtained at the ¦rst order can be written as The CHARME solver for ¤the separated two-phase §ow¥ is obtained after pressures, velocities, and temperatures relaxations. The system obtained is composed to (2 + n g + 1) equations. In this solver, a §uid mixture composed by one gaseous phase of (n 1 = n g ) species and one liquid phase (n 2 = 1) standing for the dense LOx is considered. The classical NavierStokes system writes: In this system, conservative and primitive set of variables, respectively, writes: . . , n g , stand for the mass fractions of gaseous species while Y l stands for the mass fraction of the dense LOx. Then, ρ is the mixture density, e tot is the total energy, and P and T stand, respectively, for the pressure and the temperature. The convective and di¨usive §uxes can be written as We do not give more details on di¨usion and convection which are classical phenomena in §uid mechanics. Then, source term S(U) includes combustion modeling and coupling source terms between the CHARME solver and the SPIREE solver presented in the following subsection for the description of the ¤dispersed two-phase §ow.¥

Turbulent combustion
The H 2 O 2 combustion is then modeled using an in¦nitely fast chemistry assumption (high Damk ohler number). This means that kinetics e¨ects are not taken into account. The species production rates are related to the gap between the local and the equilibrium concentration. In other words, the reacting species are relaxed towards chemical equilibrium with a ¦nite relaxation time driven by a turbulent time scale. In the LES framework, such a time scale can be assumed from the resolved strain tensor. This approach is similar to the well-known ¤eddy breakup¥ (EBU) model as in both approaches, in¦nitely fast chemistry is assumed. But fortunately, taking into account a local equilibrium involving radical species, renders a much more accurate §ame temperature. The reaction rate writes: The constant has been chosen equal to 1 because this kind of model presents the same drawbacks than the EBU model of Magnussen and Hjertager [28] in term of prediction of the behavior of the §ame and there is no universal value for the constant. We know that the combustion model would have to be improved in the future but this can be done independently of the new strategy presented in this paper. So, we do not focus on turbulent combustion modeling at this step. Then, in the numerical approach, we have used upwind ¦nite volume schemes on generalized unstructured meshes. Classically, we have used the Monotically-Integrated-LES (MILES) approach [29,30] because the upwind schemes have intrinsic subgrid turbulence models coupled naturally to the resolved scales.

Kinetic Based Model
Here, we derive the sectional model used in this paper for the description of the dispersed phase [31] and the starting point is the WilliamsBoltzmann-like kinetic equation. Modeling of dispersed two-phase §ows is based on a mesoscopic description of the dispersed phase. Particles are supposed to be spherical and fully characterized by a small set of variables: position x, radius r (or, more generally, a size variable denoted by φ), velocity v, and temperature θ. In most applications, the particle number density function (n. In this balance equation, the left-hand side stands for the ¤transport¥ of the particles in the phase space (F, R, and H, respectively, correspond to the force acting on a particle, the evaporation rate, and the heat exchange rate) while • and Q in the right-hand side, respectively, stand for the e¨ect of fragmentation and collision phenomena. Note that F, R, and H depend on the local gas composition, velocity, and temperature. All §uid models for gasparticle §ows are based on conservation equations for some particular moments of the n.d.f. These models can be formally derived from the kinetic equation by assuming particular closure assumptions. Let us ¦rst introduce the following notations: where h stands for the particle enthalpy. Integrating the ¤kinetic¥ equation with respect to θ and v , one formally gets a system of balance equations of the following form: As we want to precisely describe the polydispersion of the spray, models of order 2 in size are derived. So, in addition, another equation for the mass conservation is introduced. Then, equations for the conservation of momentum vector velocity and energy are also derived. The term P stands for the second-order kinetic stress tensor (equivalent to a generalized pressure tensor). It is due to the dispersion of the droplet velocity distribution function. In this study, it is neglected. Now, we describe the option, often called ¤multi §uid¥ model or sectional model which is illustrated in Fig. 2. Information on the droplet size distribution is kept at the macroscopic level thanks to a ¦nite volume discretization with respect to the size variable. A set of equations is derived for each section and, in this type of model, sections are coupled thanks to mass, momentum, and heat Compared to other reconstruction, the principal advantages of the linear one are to conserve the positivity when computing the inversion between the mass and number (n k and m k ) and the coe©cients (α k , β k , s min , and s max ) and to reduce the computational cost. When using this reconstruction, we also keep the possibility to perform exact computation of the integrated source terms.
The SPIREE solver for the ¤dispersed two-phase §ow¥ writes: In this system of conservation laws, the conservative and primitive variables are: In this system, ρ d stands for the particle density; N d stands for the averaged number of particles by unit of volume; h d stands for the total energy; and variables D, v d T d , and α are, respectively, the particle diameter, the velocity, the temperature, and the volume fraction. The convective §ux is given by the following expression in which the pressure term is neglected: In the following, the classical source terms between the CHARME and SPIREE solvers are brie §y presented. The source term taking into account for the primary atomization will be speci¦cally described in details in section 4.

Modeling drag force
Because of the very high density of the particulate phase with respect to the §uid one (ρ/ρ d ≈ 10 −3 ), the only external force to be accounted for is the drag force with the following general expression: where C D is the drag coe©cient depending on the particle radius as well as on the characteristic of the §uid §ow around the particle. The dynamic relaxation time writes: This velocity response time represents the time required for a particle to respond to a change in the §uid velocity.

Modeling heat and mass transfer
A large variety of models can be found for the description of mass and heat transfer between the particles and the gas. The simplest one is based on the d 2 law [32] where the droplet temperature is supposed to be constant and the droplet heating is neglected. A linear regression in time for the droplet surface is then obtained. At the opposite, much more complex models which take into account the evolution of the temperature pro¦le inside the droplet give a very accurate description of the evolution of the droplet surface temperature. A classical intermediate model resolving droplet heating but still not resolving internal conduction is the in¦nite conductivity model. In the solver SPIREE, the classical AbramzonSirignano model [33] is implemented. The droplet temperature is supposed to be uniform but to depend on time and the variation of the droplet surface φ = S obeys: This expression of dS/dt is derived from the mass conservation equations under the assumption that the gas §ow around the droplet is stationary. The parameter B M stands for the Spalding dimensionless mass transfer number. To take into account the e¨ect of the convective transport on vaporization due to the droplet motion relative to the gas, the so-called ¤¦lm theory¥ is used and the convective Sherwood number Sh is introduced. This number has to be modi¦ed (Sh * = 2 + (Sh − 2)/F M (B M )) to take into account the presence of the Stefan §ow. Note that subscripts s and ∞ refer to the conditions at the droplet surface and at in¦nity from the droplet surface, e. g., in the external gas §ow and �·� represents an averaged value between conditions s and ∞. D and Y stand, respectively, for the binary di¨usion coe©cient and the vapor mass fraction in the gas. Then if we examine the energy conservation equation, we can deduce the evolution for H which corresponds to the heat §ux entering the droplet.

Modeling secondary breakup
The general form of the fragmentation reads: where v bup (v, r) denotes the breakup frequency. When We > We c , the frequency is given by the following model with σ denoting the surface tension coe©cient and τ bup corresponding to the average breakup time: Refer to [34] for h bup standing for the n.d.f. of the droplets produced by the fragmentation of a given droplet of radius r * and velocity v * . Then, the experimental correlation for the Sauter Mean Diameter after fragmentation has been chosen following the work of Wert [35].

NUMERICAL METHODS
The numerical methods used in this study are based on a ¦nite volume approach. For the coupling between CHARME and SPIREE solvers including a large variety of phenomena, here, a time splitting technique is proposed. Concerning spatial approximation, a new MUSCL multislope technique [20] has been developed for general unstructured meshes to ensure robustness of the simulation.

Time Splitting Technique
The method consists in splitting the integration into di¨erent operators based on di¨erent physical phenomena. The principal advantage is to reduce the computational cost when characteristic times are very di¨erent but also to improve the robustness when dealing with two-way coupling simulations in which interactions between particles and carrier phase are very strong. With this technique, one can choose an optimal time step for the computation and is able to use a dedicated numerical scheme for each operator: T g and T k stand, respectively, for transport of the §uid and the dispersed phase; S is the classical coupling operator between the §uid and the spray (drag force and heat and mass transfer); I is the operator dealing with the particles interaction such as the fragmentation or secondary breakup; and A stands for the primary atomization coupling the §uid and the spray.
In this study, the Lie formalism of order 1 has been applied:

New MUSCL Multislope Technique
When someone is interested in two-phase §ows, he/she has to deal with high density ratios, strong gradients, and, also, discontinuous solutions. So, the Figure 3 Forward and backward points for MUSCL multislope technique in two dimensions discretization has to be a good compromise between accuracy and robustness. Nowadays, the MUSCL technique remains an alternative to more recent technique and for this reason, it has been chosen to build a new multislope approach. The very di©cult task was to formulate the new approach for general unstructured meshes. A new reconstruction procedure is introduced in [20], operating on general unstructured meshes and computing the reconstructed values at the face centroids. All the most important properties to ensure precision and stability are proved in [20]. For sake of stability, the primitive variables (volume fraction, velocity, temperature, and diameter) are interpolated instead of the conservative ones. Figure 3 shows a scheme of principle for the new multislope method. The technique generalizes the multislope method introduced before for triangular or tetrahedral meshes. As in the original MUSCL method, both a backward and a forward scalar slopes are computed for each face of a given element. Figure 4 illustrates the coupling strategy between the two solvers CHARME and SPIREE. The solver CHARME (CH) for separated two-phase §ow stands for the liquid dense oxygen (O CH 2,l ), the gaseous oxygen (O CH 2,g ) resulting from droplets evaporation, the gaseous hydrogen (H CH 2 ), and the combustion products (PRDTS CH ). On the other hand, the solver SPIREE (SP) stands for the dispersed oxygen liquid droplets (O SP 2,l ). Figure 5 shows a preliminary computation to illustrate the coupling strategy. We represent instantaneous ¦elds of the following variables starting from the top left-hand side and moving to the bottom right-hand side following the Figure 4 Coupling strategy between CHARME and SPIREE solvers Figure 5 Preliminary results for the coupling between CHARME and SPIREE arrows. The LOx mass fraction is determined in the §uid solver. The atomization source term shows where the LOx mass is transfered from the §uid to the spray. The total volume fraction of the spray (comprising all sections) is only resulting from this mass transfer. The total vaporization rate accounts for the inverse transfer: oxygen mass is leaving the spray solver to return to the §uid but under gaseous form due to the evaporation process. This results in the gaseous oxygen ¦eld underneath. Finally, as a result of the gaseous oxygen and hydrogen ¦elds (bottom left-hand side), we get the water mass fraction and the temperature ¦elds representing the combustion process (bottom right-hand side). The classical source terms between CHARME-SPIREE are the drag force (velocity relaxation), the heat transfer (temperature relaxation), and the mass transfer: F D , φ C , ' m vap . These terms have been described in the previous section. The correlation of SchillerNaumann is used for the drag force and the AbramzonSirignano model [33] is used for the heat and mass transfer. Now, let us describe the source terms s(u ) for each section of particles of SPIREE and then focus on the primary atomization.

COUPLING STRATEGY
The ¦rst component traduces the growth of mass due to primary atomization and its decrease by vaporization. Now, let us look at the source term S(u) a¨ected to the carrier phase CHARME. The ¦rst n g components of gas species concern combustion and evaporation phenomena. The ¦rst term includes the evaporation of liquid oxygen droplets in the ¤dispersed¥ phase which creates gas oxygen in the carrier phase of CHARME. The component number n g + 1 standing for the transport equation of mass liquid fraction concerns primary atomization. It transforms the dense LOx of CHARME into ¤dispersed¥ LOx in the appropriate sections of SPIREE:

Modeling Primary Atomization
The model used to describe the mass transfer between solvers accounting for primary atomization reads: where ρY l is the liquid mass in a given volume of control; v ato is the characteristic frequency of the primary atomization process; and λ ato (Y l ) is the e©ciency function. Let us assume the atomization frequency to be directly connected to the strength of the velocity gradient, which is the only information locally available in the LHF framework (no velocity di¨erence is known). 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, it has been chosen to use the latter approach: The e©ciency function reads: It is designed to make sure that when some LOx mass is transferred from the §uid towards the spray in a given volume of control, the corresponding vanishing volume in the §uid 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 di¨usion which spreads the interface over several mesh elements in order to carry out the mass transfer in a smooth way.
So, the properties of the created droplets resulting from the primary atomization have to be assumed. They cannot be computed locally from resolved quantities as the LHF formalism provides too little information. Actually, these properties are estimated based on the instability analysis from Marmottant and Villermaux [36]. 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, etc.). Consequently, the knowledge of the steady operating conditions of the Mascotte con¦guration 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 §uid, which may be actually a rough approximation. Note that a very accurate and local value of the created droplet diameter is not of the ¦rst-order importance, as the use of a secondary breakup model is expected to rapidly modify and somehow correct the local droplet diameters. In fact, the zone of secondary atomization is expected to be correctly computed. Concerning the zone of primary atomization, the computation is limited by the LHF model in which only one velocity is available. Finally, the temperature of the created droplet is just set to the constant value which was used to describe the liquid phase in the §uid, namely, 85 K, corresponding to the LOx injection temperature.

NUMERICAL RESULTS
In this section, some numerical results obtained with the coupling strategy are presented.
The three-dimensional geometry is depicted in Fig. 6. The overall device is approximately 50 cm long, with a 50-millimeter wide section. The LOx post has a 5-millimeter diameter, whereas the total diameter of the injector (axial LOx + coaxial H 2 ) is 12 mm. A tetrahedral unstructured mesh made up with approximately 9.8 M of elements has been used. The mesh has been built so that the ¦nest re¦nement is located near the injector exit, where atomization takes place. The smallest cell size is the order of 20 µm (in the zone marked by white square) whereas the maximum cell size is the order of 3 mm at the end of the chamber. Figure 6 also represents the mesh computation with a zoom near the injector. For the sake of simulation, the computational geometry has been split into 1920 subdomains, and then dispatched into 480 processors to allow parallel computing. The physical time step of the computation is of the order of 2 · 10 −8 s induced by a maximum Courant number imposed to 1.5. The total physical time computed is 17 ms which corresponds to 130 jobs of 15 h on 480 processors. The CPU time required for a simulation of 1 s of physi-  cal time is about 50 · 10 6 h. The numerical results obtained are presented in the rest of the section. Comparisons with experimental results provided in [26] are only qualitative because the results are not yet converged as illustrated in Fig. 7. Figure 7 stands for the time evolution of the pressure and temperature obtained by the resolution of the LHF solver as well as the volume fractions and penetration depth obtained by the resolution of the spray solver. Results appear clearly not converged at this point of the computation because they include the end of the transient regime. The mean pressure decreases in the interval 04 ms and reaches a minimum value equal to 7 bar. Then, the ¦rst droplets appear and feed the combustion which induces an increasing of the pressure. Between 4 and 13 ms, the pressure goes up to a maximum value equal to 10.7 bar and then decreases to reach 10.1 bar. In the experiments, the nominal pressure is equal to 11 bar. The maximum temperature can be related to the formation of the stable §ame. The maximum value between 3500 and 3600 K is obtained at 4 ms. Figure 7 also represents the evolution of the volume fraction for the three sections of the spray which are not empty. The volume fraction of the three sections globally increases with a total volume fraction which tends to 0.02. In Also, the penetration depth has been plotted which is evaluated in the simulation with the position of the isoline of the liquid volume fraction equal to 0.99. This length increases during the simulation such as expected and seems to stabilize between 11 and 13 mm. This length is approximately stable at 8 ms. Comparison with theoretical value of 8 mm or experimental values between 8 and 41 mm is not realistic at this level of convergence. In addition, the interpretation of the penetration depth must be carried carefully because of numerical di¨usion.
In Figs. 814, di¨erent mean ¦elds for the §uid solver are presented. The averaged ¦eld is computed between 13 and 17 ms. Each variable is represented in the plane (XZ) and the temperature and velocity (see Fig. 8) isovalues have been plotted as well as liquid and gaseous oxygen mass fractions (see Fig. 9) and the H 2 gas (see Fig. 10) and the H 2 O product of combustion (see Fig. 11). The norm of velocity is represented with logarithmic scaling. One can observe the recirculation of the coaxial H 2 with high velocity around the liquid oxygen which has a low velocity and which is atomized. In the lateral position, one can also observe the helium ¦lm. A recirculation of H 2 is also observed. The mean value of the mass fraction clearly shows the region of transition between separated and dispersed two phase §ow.   In addition, a qualitative comparison between experiments [26] and results concerning the visualization of the radical OH has been realized which is very representative of the position of the front §ame in the plane (XY ) (see Fig. 12).
For the spray solver, only some instantaneous ¦elds are presented. Figures 13  and 14 stand for the isovalues of volume fraction and the evaporation rate in each section. Figure 15 stands for the frequency of fragmentation at ¦nal time 17 ms and for each section and in the plane (XZ). Figure 16 represents the total volume fraction as well as the primary atomization source term in the plane (XY ). Figure 17 stands for the time evolution of the liquid dense oxygen surface as well as the temperature (isovalues). This ¦eld gives one a good approximation of the location of the front §ame. Finally, in Fig. 18a, qualitative comparison between experiments and numerical results for the instantaneous ¦eld is given.

CONCLUDING REMARKS
In this paper, a LES of the primary atomization in cryogenic combustion chamber has been proposed. To reach this crucial challenge, a fully Eulerian coupling strategy between a di¨use interface 4-equation (1-velocity) model and a kinetic based model has been suggested. Speci¦c numerical methods have been developed to ensure accuracy and robustness of the computation. The ¦rst results seem to be very promising but need to be converged. For this reason, the comparisons with experiments are only qualitative at this moment. In the future, we intend to use a 7-equation (2-velocity) model in order to improve the physical modeling of the primary atomization.