A Ruiz

L Selle

B Cuenot

T Poinsot

STABILIZATION OF A SUPERCRITICAL HYDROGEN / OXYGEN FLAME BEHIND A SPLITTER PLATE

The numerical simulation of fluid dynamics and combustion in cryogenic rocket engines is addressed in this paper, with the intent to elucidate flame stabilization mechanisms. A model configuration is devised to allow a fully resolved simulation, both for the dynamics and the flame structure: a twodimensional splitter plate represents the lip of an injector and the operating point is typical of a real engine. The non-reacting flow field is first scrutinized to evaluate the impact of the large density gradients between the fuel (hydrogen) and oxidizer (oxygen) streams. It is found that the turbulence generated by the splitter is very intense and strongly distorts the high-density-gradient front at both small and large scales. Under reacting conditions, the flame stabilizes right at the lip of the injector, which is a common feature of hydrogen / oxygen flames under these conditions. A particularly complex flame structure is evidenced at the anchoring point, with turbulent transport playing an important role.

Introduction

The development of liquid-propellant rocket engines is a long and costly trial-and-error process. One of the major challenges in the design of the combustion chamber is to ensure a stable combustion. Indeed, the occurrence of combustion instabilities can cause structural damage but also raises critical issues for mission success, as well as crew safety. Another critical point is to ensure a steady combustion regime, avoiding flashback and blowoff events, for example. A key parameter controlling these phenomena is the stabilization of the flame downstream the injector, which is the subject of the present study. The literature on flame stabilization is extremely large, with detailed investigations using both experiments and numerical simulations. The main parameter describing flame stabilization is the distance from the injector to the flame front and a wide variety of phenomenon have been reported with chemistry, transport and various sources of heat losses potentially playing an important role.

Cryogenic rocket engines operate at a very large pressure, at which reactants are usually under a supercritical state. In such conditions, the phase-change phenomena no longer occurs and fluids transition from a liquid-like state to a gas-like state in a continuous way. However, because the equation of state under such conditions is highly non-linear, very large density gradients due to temperature or composition variations are common. It has been shown that these gradients affect the overall dynamics of jets as well as the development of turbulence, mixing processes and eventually combustion [START_REF] Yang | Modeling of supercritical vaporization, mixing, and combustion processes in liquidfueled propulsion systems[END_REF][START_REF] Bellan | Theory, modeling and analysis of turbulent supercritical mixing[END_REF]. This particular behavior is a major issue for the development of accurate and quantitative numerical simulations. However, over the past 10 years, a number of research groups were able to conduct detailed simulations under fairly complex conditions, with developed turbulence [START_REF] Okong'o | Direct numerical simulations of O2/H2 temporal mixing layers under supercritical conditions[END_REF][START_REF] Selle | Modelling of subgrid-scale phenomena in supercritical transitional mixing layers: An a priori study[END_REF]Zong & Yang, 2008) and sometimes combustion [START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF][START_REF] Palle | Analysis of high-pressure hydrogen, methane, and heptane laminar diffusion flames: Thermal diffusion factor modeling[END_REF].

The objective of this paper is to present a fully-resolved numerical simulation of a model injector to scrutinize turbulence, mixing and combustion processes downstream the injector. The numerics and models are first presented (Sec. 2) followed by the description of the configuration (Sec. 3). Then the validation of the chemistry is presented in Sec. 4, followed the analysis of the cold and reacting flows (Sec. 5).

Numerics and Models

The compressible Navier-Stokes equations are solved using the AVBP code (Gourdain et al., 2009b,a). In the present study, a two step Taylor-Galerkin scheme called TTG4A is used, which is third order in space and time [START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF]. Real-gas thermodynamics is accounted for through the Peng-Robinson equation of state [START_REF] Peng | A new two-constant equation of state[END_REF] while transport coefficients are modeled based on the theory of corresponding states for the dynamic viscosity and the thermal conductivity [START_REF] Chung | Applications of kinetic gas theories and multiparameter correlation for prediction of dilute gas viscosity and thermal conductivity[END_REF][START_REF] Chung | Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties[END_REF] and constant Schmidt numbers (c.f. Tab. 1). This overall numerical and modeling methodology has already been validated for non reacting flows [START_REF] Schmitt | Large-eddy simulation of supercritical-pressure round jets[END_REF]. The critical-point coordinates of the intermediate species OH, O, H, H 2 O 2 , HO 2 (for which no experimental data is available) is estimated using the Lennard-Jones potential-well depth σ i , and the molecular diameter i of the i-th species from the CHEMKIN transport coefficients of the San Diego Mechanism, according to the following expression [START_REF] Giovangigli | Detailed modeling of planar transcritical H2-O2-N2 flames[END_REF]:

V c,i = 3.29N σ 3 i (1) T c,i = 1.316 i k (2) (3)
where N is the Avogadro number and k is the Boltzmann constant. The acentric factor ω ac is set to zero for radical species. The numerical values for all species are summarized in Tab. 1. The configuration is a two-dimensional H 2 /O 2 flame stabilized behind a splitter plate at a mean pressure of 100 bar. The computational domain, presented in Fig. 1, is 11h long in the x-direction and 10h in the y-direction, where h = 500 10 -6 m is the splitter height. Hydrogen is injected above the splitter at a temperature T in H 2 = 150 K and a bulk velocity u in H 2 = 125 m/s. Below the splitter, oxygen is fed at T in O 2 = 100 K and u in O 2 = 30 m/s. These conditions were chosen to mimic a typical liquid rocket engine, at the nominal operating point. The shape of the inlet velocity profiles follows a 1/7th power law. Although developed turbulence is generally present in the feeding lines of rocket engines, here, no velocity perturbation is added through the inflow boundary condition. And yet, strong turbulence levels caused by vortex shedding are observed downstream the splitter allowing for a developed turbulent mixing layer and strong flame / turbulence interactions (c.f. Sec. 5). The outlet boundary condition is derived from the NSCBC technique [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF][START_REF] Baum | Accurate boundary conditions for multicomponent reactive flows[END_REF] and accounts for both real-gas effects (Okong'o & Bellan, 2002) and transverse terms [START_REF] Granet | Comparison of nonreflecting outlet boundary conditions for compressible solvers on unstructured grids[END_REF]. The upper and lower boundaries are treated as symmetries while the splitter is an adiabatic no-slip wall. The mesh resolution is ∆ = h/500 in a layer containing the splitter, with a 3h vertical extent. Outside this zone, a transverse stretching factor of approximately 1.02 is employed. The mesh contains approximately 13.5 million nodes.

Parameters

Chemistry

The combustion of hydrogen and oxygen is modeled using a detailed scheme accounting for 8 species and 12 reactions [START_REF] Boivin | An explicit reduced mechanism for H2-air combustion[END_REF], which is derived from the San Diego mechanism [START_REF] Petrova | A small detailed chemical-kinetic mechanism for hydrocarbon combustion[END_REF]. The forward rate coefficients are given in Tab. 2. The backward reaction rates are classically computed using low-pressure entropy and enthalpy NIST/JANAF tables (Jan, 1998). A more general treatment of the reaction rates in the real-gas framework, using chemical potentials directly a Units are mol, s, cm 3 , kJ, and K. b Chaperon efficiencies are 2.5 for H 2 , 16.0 for H 2 O and 1.0 for all other species. c Chaperon efficiencies are 2.5 for H 2 , 12.0 for H 2 O and 1.0 for all other species. d Chaperon efficiencies are 2.0 for H 2 , 6.0 for H 2 O and 1.0 for all other species. computed from the equation of state has been investigated by [START_REF] Giovangigli | Detailed modeling of planar transcritical H2-O2-N2 flames[END_REF], however, the impact on the flame structure appears to be rather small, showing that the perfectgas treatment of the reaction rates is a good approximation: this could be due to the high temperature (and hence perfect-gas behavior) at the flame location. Note that since high pressure is considered in the present study, the high-pressure limit of the falloff reactions have been taken.

The validation of the implementation in AVBP is achieved by comparing the flame structure using CANTERA [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] and AVBP in a counterflow flame configuration. The numerical setup is is that of an opposed-jet flame [START_REF] Pons | The structure of multidimensional strained flames under transcritical conditions[END_REF] computed, in AVBP, on a simple square mesh (size h) with a constant grid size identical to that of the splitter case. The left boundary condition is a symmetry, the right side an outlet, while hydrogen comes from the top and oxygen from the bottom. The boundary velocity are determined to impose a constant strain, a, on the flame:

u O 2 (x) = ax (4) v O 2 (y) = -ay (5) u H 2 (x) = ax ρ O 2 ρ H 2 1 2 (6) v H 2 (y) = -ay ρ O 2 ρ H 2 1 2 (7) (8)
The hydrogen velocity is chosen so that the momentum flux ρu 2 is the same between the oxygen and hydrogen stream, to place the stagnation point in the middle of the computational domain. As for the value of the strain, a, the validation is conducted at a value based on the splitter height h and the mean velocity difference between the two streams at the injection: this yields a = 3800 s -1 . A typical result of such strained diffusion flame is presented in Fig. 2 showing streamlines superimposed on the temperature field. The validation procedure comprises three computations:

1. CANT PG: a computation is performed using CANTERA using the perfect-gas equation of state.

The temperature of the fresh gases is 300 K to limit real-gas effects.

2. AVBP PG: this is the same computation as CANT PG, performed in AVBP. Because CANTERA solves for the full transport matrices while AVBP assumes constant Schmidt numbers and mixtureaveraged transport coefficients, this computation is both a validation of the implementation of the chemistry in AVBP as well as a validation of the simplified transport.

3. AVBP RG: for this computation, the temperature of the fresh gases is lowered to match those of the splitter case and evaluate real-gas effects on the flame structure.

Then, a cut through the flame at x = h/2 in the computational domain of AVBP is compared to the flame structure from CANTERA. Figure 3 shows the temperature and mass fraction of HO 2 (the main initiator of the combustion) versus the mixture fraction for the three computations. First comparing simulations CANT PG and AVBP PG, the agreement between the two codes is excellent: the maximum discrepancy for the temperature is of the order of 70 K, which gives a relative error of 2 %. The differences for the mass fraction Y HO 2 are even smaller. This agreement validates the implementation in AVBP as well as the assumptions on the transport, for this temperature. The laminar flame in the thermodynamic conditions of the splitter case (AVBP RG) is very similar to the higher-temperature computation (AVBP PG) in the flame region. This conclusion is similar to that of Ribert et al. [START_REF] Ribert | Counterflow diffusion flames of general fluids: Oxygen/hydrogen mixtures[END_REF], which they justify by the fact that the combustion processes take place in hot regions where real-gas effects are negligible.

Results and analysis

This section is devoted to the computation of the splitter configuration presented in Sec. 3. First, the non-reacting flow is analyzed with specific attention on the variety of scales caused by the interaction between turbulence and the large density gradients. Combustion is then initiated using the kinetics validated in Sec. 4.

Cold Flow

The spatial resolution in the present study is such that 500 grid points are used to discretize the splitter height h. As a comparison with previous efforts to compute such configurations [START_REF] Oefelein | Thermophysical characteristics of shear-coaxial LOX-H2 flames at supercritical pressure[END_REF], the resolution is ten times greater. As the resolution is increased, the range of flow structures widens considerably, especially in the density field as shown in Fig. 4. The first conclusion from this instantaneous visualization of the flow is that strong turbulence levels develop right downstream the splitter, which is an a posteriori validation for the lack of turbulence injection. Small vortices are regularly shed from to the edges of the lip, due to the shear-induced instability. Then many pairing events occur thatresult in large vortices with a length scale comparable to that of the splitter height. While impacting the LOX stream, these vortices roll-up the LOX interface, resulting in the apparition of large waves-like structures at the boundary of the LOX stream. These density waves interact with each other and merging can occur, creating pockets where intense turbulent mixing between O 2 and H 2 is present. It is speculated that the wave structures identified in the density-gradient surface could correspond to the "finger-like" structures observed in many experimental studies of transcritical flows [START_REF] Chehroudi | Visual characteristics and initial growth rate of round cryogenic jets at subcritical and supercritical pressures[END_REF][START_REF] Oschwald | Injection of fluids into supercritical environments[END_REF], using shadowgraphs. The density field is also clearly distorted at much smaller scales with pockets of dense fluid eventually extracted from the dense stream. In the absence of surface tension, these small pockets do not form droplets and are therefore further broken down by turbulence.

The scalar dissipation rate χ = 2D ∇Z.∇Z (where Z is the mixture fraction) is a measure of the molecular mixing between hydrogen and oxygen. This mixing is favorable for combustion and is greatly enhanced by turbulence: the wrinkled exchange surface enables a larger total mass flux and the compressive strain locally steepens the composition gradients, which is also favorable for species diffusion. Figure 5 shows a snapshot of χ where, again, the intense turbulent activity downstream the splitter is evident.

The comparison of Figs. 4 and 5 also suggests that the majority of the scalar dissipation occurs on the side of the light hydrogen stream. This would indicate that this supercritical jet behaves like a liquid jet in the sense that first, some oxygen must leave the dense stream (through evaporation in the case of a liquid jet) and is then mixed with the lighter species.

Reacting Flow

Combustion is initiated in the computational domain with the following procedure: first the field of mixture fraction is computed from a non-reacting established flow, then, the flame structure from a laminar diffusion flame is applied onto this field with subsequent replacement of the composition and temperature. During this procedure, the cold-flow pressure and velocity fields are not altered. As expected, significant acoustic perturbations are generated when the flow adapts from this approximate solution based on an educated guess, however, the procedure is successful because the pressure waves eventually leave the computational domain and a stable combustion regime is reached.

Flame stabilization. An instantaneous temperature field is presented in Fig. 6 with superimposed iso-contours of density gradient (green) and heat release (black and grey). The first observation is that the flame is not lifted and the hot gases touch the splitter plate. This is certainly partially due to the adiabatic boundary condition imposed at the solid boundary. However, hydrogen / oxygen flames are known to be extremely reactive and experimental evidence under comparable thermodynamic conditions showed that the flame is stabilized very close to the splitter [START_REF] Singla | Flame stabilization in high pressure LOx/GH2 and GCH4 combustion[END_REF]. The temperature field is very convoluted on both oxygen and hydrogen sides and similarly to what was observed in the cold flow, pockets of oxygen are sometimes captured in the hot gases, especially in the second half of the computational domain. This observation allows to speculate that these pockets of warmed-up oxidizer could eventually enter the hydrogen stream and generate a corresponding unsteady heat release. Further analysis of the different combustion regimes is required to substantiate this claim. The green iso-contour highlights the region of high density gradient and allows to distinguish the dense stream from the lighter regions. Contrary to the observations on the cold flow, the large-scale wave-like structures are not present in the reacting regime; only smaller convolutions are observed, at least within the computational domain. Such reduction in the vertical extent of the mixing layer under reacting conditions has already been observed at low pressure [START_REF] Renard | Dynamics of flame/vortex interactions[END_REF]. Finally one can note that, as expected from the cold-flow scrutiny, mixing and combustion occur in the light stream.

A close-up view of the vicinity of the splitter is presented in Fig. 7. Again, the burnt gases touch the solid boundary, which is consistent with the adiabatic treatment of the splitter, but not realistic. Regions of strong reaction rate are visible at the edges of the high-temperature pocket, both on the hydrogen and oxygen sides. This is an indication that turbulent mixing is very intense in this region allowing each species to be transported across the splitter height. At this point in the analysis, this transport mechanism is not clear but animations of the splitter region show that large-scale vortices allow for fresh hydrogen to be transported towards the oxygen stream while at other times, burnt gases with excess oxygen are transported up to the hydrogen stream. It is therefore concluded that even though chemistry plays a key role in the stabilization of such flame, turbulent transport is also active and influences the combustion regime. It is also foreseen that for fuels with less chemical activity, turbulence could become a key stabilization process.

Flame structure. In an attempt to elucidate the combustion regimes behind the splitter, a scatter- Curvature of the density field. In order to quantify the visual differences observed in the reacting and non-reacting cases, the probability density function (PDF) of the curvature of a density iso-surface is presented in Fig. 9. The value chosen for the curvature analysis is the iso-surface at the intermediate density between the hydrogen and oxygen streams:

ρ 1/2 = (ρ in H 2 + ρ in O 2 )/2 ∼ 600kg.m -3
. While the PDF for the reacting flow is roughly symmetric and centered around a zero value, the curvature PDF for the cold flow is double peaked: the peak at -10 4 m -1 corresponds to a concave radius of curvature of h/5 in the density isocontour. This indicates that although visually, the large density waves appears to have a characteristic radius of curvature of h, they actually consist in a collection of highly curved pieces. The difference between the hot and cold flow curvature pdf also seems to confirm the visual observation that the density waves disappear in the hot flow.

Conclusions

Simulations of a splitter-plate configuration representative of a cryogenic rocket engine injector have been conducted. The dynamic and thermodynamic conditions are chosen to be representative of a real engine under steady operation: both fuel (hydrogen) and oxidizer (oxygen) are above their critical pressure. Both non-reacting and reacting configurations were computed, on a two-dimensionnal domain with full resolution of turbulence and chemistry. The cold flow shows that intense turbulence develops downstream the splitter resulting in a highly convoluted region of large density gradient. This interface exhibits both large and small-scale structures similar to experimental direct visualizations. It is also reported that mixing occurs mainly in the light hydrogen stream, which is confirmed under reacting conditions as the flame develops outside the dense oxygen stream. Under reacting condition, the flame is anchored on the splitter, which is due to both the adiabatic boundary condition and the high reactivity of hydrogen. A complex stabilization process involving both chemistry and turbulent transport is suggested but further investigations are required to confirm and quantify this claim. Finally, the larger scales in the high density-gradient magnitude regions seem to disappear under reacting conditions. This work is to be extended to boundary conditions accounting for heat losses at the splitter in order to investigate stabilization mechanisms in a more accurate way.

Figure 1 :

 1 Figure 1: Computational domain for the two-dimensional splitter-plate configuration.

Figure 2 :

 2 Figure 2: Strained diffusion flame computed in AVBP: streamlines superimposed on the temperature field. Thermodynamic conditions correspond to the splitter case: hydrogen at 150K from the top, oxygen at 100K from the bottom and ambient pressure is 100 bar.

Figure 3 :

 3 Figure 3: Comparison of flame structure versus mixture fraction between AVBP and CANTERA. (a) temperature and (b) HO 2 mass fraction. The vertical bar indicates the stoichiometric mixture fraction.

Figure 4 :

 4 Figure 4: Density field, ranging from 15 kg.m -3 (blue) to 1200 kg.m -3 (black) using a logarithmic scale.

Figure 5 :

 5 Figure 5: Scalar dissipation rate, χ, in the cold flow (same physical time as in Fig. 4).

Figure 6 :

 6 Figure 6: Temperature field with superimposed density gradient (green) and heat-release (black: max heat release of case AVBP RG; grey: 10 % of case AVBP RG).

Figure 7 :Figure 8 :

 78 Figure 7: Close-up view of the flame stabilization zone behind the splitter. Temperature field (yellow is 3800 K)with superimposed density gradient (green) and heat-release (black: max heat release of case AVBP RG; grey: 10 % of case AVBP RG).

Figure 9 :

 9 Figure 9: Curvature PDF of the median density (ρ 1/2= (ρ in H 2 + ρ in O 2 )/2 ∼ 600kg.m -3) isocontour for the cold and reacting flows.

Table 1 :

 1 Species critical-point properties (temperature T , pressure P , molar volume V and acentric factor ω) and Schmidt numbers.

		H2	O2	H2O	O	H	OH	H2O2	HO2
	T c,i (K)	33	154.581 647.096 105.28 190.82 105.28 141.34 141.34
	P c,i (MPa)	1.2838 5.0430 22.064 7.0882 31.013 7.0883 4.7861 4.7861
	V c,i (cm 3 /mol)	64.284 73.368 55.948 41.205 17.069 41.205 81.926 81.926
	ω ac	-0.216 0.0222 0.3443	0.0	0.0	0.0	0.0	0.0
	Schmidt Number 0.28	0.99	0.77	0.64	0.17	0.65	0.65	0.65

Table 2 :

 2 +H 2 → H 2 O 2 +H 1.62 10 11 0.61 100.14 12 H 2 O 2 +M → 2OH+M d 2.62 10 19 -1.39 214.74 Forward rate coefficients in Arrhenius form k = AT n exp (-E/RT ) for the skeletal mechanism.

	Reaction		A a	n	E a
	1 H+O 2	OH+O	3.52 10 16 -0.7	71.42
	2 H 2 +O	OH+H	5.06 10 4	2.67	26.32
	3 H 2 +OH	H 2 O+H	1.17 10 9	1.3	15.21
	4 H+O 2 +M → HO 2 +M b 4.65 10 12 0.44	0.0
	5 HO 2 +H → 2OH	7.08 10 13 0.0	1.23
	6 HO 2 +H	H 2 +O 2	1.66 10 13 0.0	3.44
	7 HO 2 +OH → H 2 O+O 2	2.89 10 13 0.0	-2.08
	8 H+OH+M H 2 O+M c 4.00 10 22 -2.0	0.0
	9 2H+M H 2 +M c	1.30 10 18 -1.0	0.0
	10 2HO 2 → H 2 O 2 +O 2	3.02 10 12 0.0	5.8
	11 HO 2

Acknowledgments

The funding for this research is provided by Snecma, which is the prime contractor for the European launcher Ariane 5 cryogenic propulsion systems and CNES (Centre National d'Etudes Spatiales), which is the government agency responsible for shaping and implementing France's space policy in Europe. Their support is gratefully acknowledged. This work was granted access to the high-performance computing resources of CINES under the allocation 2011-c2011025082 made by Grand Equipement National de Calcul Intensif.