Numerical simulation of supercritical CH4/O2 combustion

The modeling and simulation of CH4-O2 combustion at high pressure requires a dedicated kinetic mechanism, such as RAMEC, to obtain a precise description of the fuel decomposition. Nevertheless, to reduce its CPU cost, the optimized and reduced chemistry method ORCh is applied and validated for a set of canonical test-cases. A very good agreement is obtained by comparison with the original RAMEC detailed mechanism. To further evaluate the efﬁciency of the reduced mechanism, a premixed ﬂame expanding in a homogeneous isotropic turbulence (HIT) at high pressure along with the combustion of pocket of methane in an environment of pure oxygen, with or without HIT, are presented. The new reduced chemical mechanism performs well with a relative error of 4% for the maximum of temperature in the premixed case and 6% for the non-premixed one.


Introduction
The design of the future European rocket will require the use of methane as fuel [8]. To achieve this objective the use of numerical simulations will be necessary, and the choice of a validated kinetic scheme usable on supercomputers becomes a key point of the whole validation process. Very few mechanism are available in the literature [20,14,22] for such extreme thermodynamic conditions, i.e. high-pressure and cryogenic inlet temperatures. For example, the mechanism GRI3.0 [20] is only valid for a temperature T ∈ [1000, 2500] K, a pressure P ∈ [0.001, 10] bar and an equivalence ratio Φ ∈ [0. 1,5]. In addition, the chemical reactions involving the species CH3O2, CH3O, HO2 and H2O2 play a major role in the decomposition of methane at high pressure [22], and must therefore be preserved. Such important species for high-pressure combustion appears in the RAMEC mechanism [14] or in the one developed by Zhukov and Kong [22]. However, simulating the combustion with a detailed or skeletal kinetic mechanism is a recurring issue [2], because the computational cost becomes prohibitive due to the required high-pressure thermodynamics [13,17], and the number of chemical species involved. This specific point is illustrated in [3] where numerical simulations are carried out with three different kinetic schemes for the combustion of methane, decane and iso-octane. These three chemical schemes have 29, 91 and 857 reactive species (N s), respectively. All other things being equal, the switching from the combustion of methane to the combustion of decane (multiplication by 3 of N s) increases the calculation time of one iteration by a factor of 6, and by 614 if the methane is replaced by iso-octane (multiplication by 30 of N s). A reduction in the kinetics must then be achieved to decrease the CPU cost that is particularly demanding when studying high pressure flames [5] with a fully compressible numerical code. Obtaining a reduced CH4/O2 combustion mechanism is therefore the main objective of this study. Starting from the RAMEC mechanism that contains 38 species and 190 reversible reactions [14], a chemical reduction is performed using the ORCh procedure [7]. The performances of the reduced mechanism are examined through the computation of auto-ignition delays, as well as the structure of laminar premixed and non-premixed flames, in order to delimit its domain of validity. This reduced mechanism is then compared to the detailed mechanism on a turbulence-flame interaction configuration, at high pressure. For this specific case, the choice of pressure, i.e. 56 bar, was guided by the experience of Singla et al. [19], who studied LOx-CH4 combustion with a coaxial injector, because this configuration is actually representative of the conditions encountered in rocket engines. Finally, a single isolated pocket of methane surrounded by oxygen at 56 bar is fired using the reduced and detailed mechanisms. These two simulations are performed using the numerical solver SiTCom-B [10,11,15,16] that solves the fully compressible Navier-Stokes equations.

Reduced Chemistry
The Optimized and Reduced Chemistry (ORCh) method [7] is used to reduce the detailed chemistry of RAMEC. Briefly, species and reactions that contribute poorly to the production of defined important (target) species are suppressed by directed relation graph methods introduced by Pepiot et al. [9] once stochastic trajectories are computed [7]. Following this procedure, a reduced mechanism having 17 species (17S) and 44 reversible reactions (44R) is obtained. Further reducing this kinetic scheme could still lead to a good description of methane oxydation, but over a more restricted range of pressures and equivalence ratios, which is not the objective of the present work.The validation of this new mechanism is now addressed by comparison with the original detailed mechanism on auto-ignition cases, freely propagating laminar premixed 1D flames (not shown here) and counterflow diffusion 1D flames (not shown here). The Cantera software [1] is used for computations. A constant pressure reactor is simulated for an equivalence ratio varying from 0.4 to 14 and a pressure and temperature set to 100 bar and 1200 K, respectively. Defining the ignition delay time, τ ign , as the time required for a temperature increase of 300 K, i.e. until reaching here 1500 K, we can see in Fig. 1 that the ignition delay time is very well captured by the reduced mechanism. Thus a good agreement is found between the reduced kinetic and the detailed one, on canonical flame configurations.

Theoretical formulation
Direct numerical simulations are carried out with the finite volume code SiTCom-B (cartesian mesh) and the fully coupled conservation equations of momentum, species and energy for reactive flows are given by: T is the vector of conservative variables. u i are the velocity components, E is the total non chemical energy decomposed as the sum of sensible energy and kinetic energy. Y k is the k th chemical species in a pool of N sp species. F j c represents the conservative variables flux vectors and is given by The source term vector S = (0, 0,ω E ,ω k ) T witḣ ω k the chemical rate of species k andω E , the heat release due to combustion. τ ij is modeled with the Newtonian viscous stress tensor, and q j , the j th component of the heat flux vector, is given for a multi-component flow in [5]. For supercritical simulations, the Soave-Redlich-Kwong (SRK) equation of state (EoS) links pressure, p, temperature, T and density, ρ, as instead of the ideal EoS, which is unable to predict correctly dense fluids. R u is the universal gas constant and W the molecular weight of the fluid mixture, W = Nsp k=1 X k W k with X k the mole fraction of species k (among N sp species). The Van der Waals mixing rules are used for mixtures: κ ij is the binary interaction coefficient. The constants a i (attractive forces), b i (co-volume of particles) and α i are determined from universal relationships [12], involving the critical temperature and partial pressure of species k, T c k and p c k , respectively, as well as the acentric factor. The evaluation of such coefficients at high pressure is based on kinetic theory considerations for dense gas mixtures. They are expressed as D bin kl = F c × D bin,• kl , with D bin,• kl , the binary diffusion coefficients at standard pressure given by the kinetic theory of dilute gases [6], i.e. using polynomial fits, or coming from the empirical correlation of Fuller et al. [12,4]. F c is a correction factor that takes into account for volume effects. Using the Takahashi approach [21], the scaling factor F c is tabulated as a function of the reduced temperature and pressure.

Turbulence-Premixed Flame Interaction
To further evaluate the efficiency of the reduced mechanism to describe the combustion of CH4/O2 at high-pressure, a flame evolving in a turbulent field is simulated. The configuration under study corresponds to the superposition of a homogeneous stoichiometric CH4/O2 mixture at 300 K with a homogeneous isotropic turbulence (HIT) generated by the Passot-Pouquet spectrum. HIT reference values are for characteristic parameters: speed of sound (c 0 = 356.3 m/s), kinematic viscosity (ν = 2.983 10 −7 m 2 /s), turbulent velocity (U p = 3.5 m/s) and the length of the most energetic scale (Le = 1.0 10 −5 m). This setting leads to approximately seven integral length scales in a box of 0.1 mm length. A uniform structured mesh is used and each direction is discretized with 400 mesh cells. With a mesh resolution of 0.25 µm, the flame thickness (≈ 2×10 −6 m) is described with 8-10 mesh cells. The mixture is ignited with a hot spot defined by hyperbolic tangents for the temperature and species mass fractions and located in the center of the domain. The maximum temperature is 2500 K and the composition of burnt gases is close to stoichiometric conditions. The pressure is set to 56 bar. The radius of the ignition sphere is 1 × 10 −5 m. The time evolution of the initial hot spot is shown in Fig. 2 by means of the heat release rate (HRR) using the reduced kinetic scheme (RAMEC-17S-44R) and the ideal gas state equation. This hot gas pocket makes it possible to start combustion, i.e. activate Arrhenius' laws, and thus generate a premix flame. This flame then interacts with the isotropic homogeneous turbulence present in the reactive mixture. The initially circular hot spot is deformed and the flame front progresses in all directions. For t = 1.4 × 10 −6 s, the species composition and temperature at r = 0 are close to values found in the calculations of one-dimensional premixed laminar flames. The behavior of this flame simulated in two-or three-dimensions is similar, but the wrinkling of the flame is less pronounced in 2D. The choice of models changes only very slightly the average behavior of the flame evolution as shown in figure 3, in which we have compared a detailed and reduced chemistry, as well as a thermodynamics of the real gas and ideal gas type. This is confirmed by comparing the flame lengths evolving over time (Fig. 4), for the different cases studied: 2D or 3D, real gas (RG) or ideal gas (IG), reduced or detailed chemistry. In this case, two flame perimeters, one external and the other internal, are calculated to define the flame length (L f ) which corresponds to the average of these two perimeters. For the 3D simulations, the median plane of the normal z is used to carry out these measurements, and the set of results, 2D and 3D, is plotted in figure 4. Comparing the results coming from 2D and 3D simulations is of limited interest because the cut along z can only give partial information for L f , and it would not be less wrong to analyze another cut with median plane of the normal x or y. In the present case, L f is always shorter for 3D simulations than for 2D simulations. In Fig. 4, the flame length values from calculations using reduced chemistry always underestimate the values of profiles using detailed chemistry. This is due to the slight difference observed in the auto-ignition delay between the two kinetic schemes. We also note that the impact of thermodynamic modeling, i.e. ideal gas or real gas, is here negligible. Finally, we note that the difference between the two kinetic schemes tends to be reduced over time. This means that a simulation of flames stabilized by means of a plate splitter for example, or with a coaxial injector [18], can be very well captured by the reduced kinetic, the errors being maximum during the transient phenomenon.

Partially Premixed Combustion
In the previous case, a premixed flame was tracked over time, with no chance of observing a non-premixed combustion. However, the latter regime is preferred for the combustion occurring in rocket engines. One way of approaching this reality, but keeping an achievable configuration in order to validate the new reduced scheme, is to consider a pocket of methane surrounded by oxygen, and to analyze the mode of propagation of the flame once the mixture is fired. A triple flame should appear, the lean branch on the oxygen side and the rich branch on the methane side. Between the two, a diffusion zone must also exist. This benchmark will thus be decisive in the validation of the reduced kinetic scheme. In Fig. 5, a pocket of methane having a diameter of 3.0 × 10 −5 m is fired using a hot spot of diameter 1.5 × 10 −5 m. A two-dimensional uniform structured mesh is used (box 0.15 mm long) and each direction is discretized with 600 mesh cells, i.e. the mesh resolution is 0.25 µm. The pressure is set to 56 bar but only the ideal gas equation of state is considered. In Fig. 5, a twin flame spreads all around the pocket of methane from right to left before merging again. Three reaction zones are clearly visible on the field of heat release rate : one is on the methane side, the other on the oxygen side, and the latter between the two first and located on the stoichiometric isoline and corresponding to a diffusion flame. The field of temperature develops around the pocket of methane and reach a value close to the one observed in 1D counterflow flames. The mass fraction of methane decreases with time and diffuse to the right, towards the hot burnt gases. Fig. 6 shows that the reduced chemical mechanism is able to reproduce the behavior of the detailed chemical mechanism, either with or without adding any HIT as initial condition, similarly to the premixed case. The flame front is a little bit delayed in time, this behavior being already observed in previous premix simulations. It is then consistent to see the same effect on the flame propagation through the domain. In the case of an addition of HIT to the initial condition (Fig. 6, bottom), a triple flame structure is still observed. This flame propagates on either side of the methane pocket following a turbulent interface. The methane pocket moves globally towards the lower boundary condition, and the code is stopped as soon as the limit is reached in order to avoid non-physical behaviors. Moreover, as shown in Fig. 7, the maximum temperature evolution is well predicted by the reduced chemical mechanism, with a maximum relative error about 6% with or without HIT.

Conclusion
To reduce the computational time of numerical simulations, the detailed kinetic schemes, required to fully describe the oxidation of hydrocarbon fuels (methane here), can be reduced by removing elementary reactions and chemical species. The RAMEC chemical mechanism was  thus reduced with the ORCh method from 38 species and 190 reactions to 17 species and 44 reactions. In order to verify the good performance of this new reduced scheme, since it must keep a high level of precision over a wide range of pressure and equivalence ratio, several one-dimensional test cases have been simulated and a very good agreement between the detailed reference chemistry and the reduced one was found. A premixed flame interacting with homogeneous isotropic turbulence has been successfully simulated using the new reduced mechanism. This is because the flame topology is similar to that found when using the detailed mechanism, but time shifted which results from a slight difference in the auto-ignition delay. A good agreement is found for the flame length between the two mechanisms, and the few observed differences disappear over time, which means that the reduced kinetic can be used to simulate an established flame. Finally, the combustion of a pocket of methane into an environment of pure oxygen is simulated with success with the reduced mechanism. The resulting triple flame, observed either with or without HIT, is recovered meaning that the combustion of