Combustion Regimes in Turbulent Non-Premixed Flames for Space Propulsion

Abstract

Direct numerical simulations of non-premixed fuel-rich methane–oxygen flames at 20 bar are conducted to investigate the turbulent mixing burning of gaseous propellants in rocket engines. The reacting flow is simulated by using an EBI-DNS solver within an OpenFOAM frame. The transport of species is resolved with finite-rate chemistry by using a complex skeletal mechanism that entails 21 species. Two different flames at low and high Reynolds numbers are considered to study the sensitivity of the flame dynamics to turbulence. Regime markers are used to measure the probability of the flow to burn in premixed and non-premixed conditions at different regions. The local heat release statistics are studied in order to understand the drivers in the development of the turbulent diffusion flame. Despite the eminent non-premixed configuration, a significant amount of combustion takes place in premixed conditions. Premixed combustion is viable in both lean and fuel-rich regions, relatively far from the stoichiometric line. It has been found that a growing turbulent kinetic energy is detrimental to combustion in fuel-rich premixed conditions. This is motivated by the disruption of the local premixed flame front, which promotes fuel transport into the diffusion flame. In addition, at downstream positions, higher turbulence enables the advection of methane into the lean core of the flame, enhancing the burning rates in these regions. Therefore, the primary effect of turbulence is to increase the fraction of propellants burnt in oxygen-rich and near-stoichiometric conditions. Consequently, the mixture fraction of the products shifts towards lean conditions, influencing combustion completion at downstream positions.

1. Introduction

Methane is a promising propellant for the new generation of space propulsion systems [1,2]. Compared to kerosene, this fuel has superior cooling capabilities and a specific impulse [3]. In addition, it produces a negligible amount of soot [4] and displays low toxicity [2]. These characteristics confer a competitive edge to methane over the established propellants in the space industry, such as hydrazine or the previously mentioned kerosene [2,5]. Consequently, substantial research activities have been fostered over the last decade to investigate the usage of green propellants [6,7,8], particularly methane for space propulsion applications. Nevertheless, critical aspects related to turbulent ignition, mixing, and combustion challenge the development of new thrust chambers propelled by methane [9,10]. Computational fluid dynamics (CFD) software provides essential support in this process. The conventional methods applied therein, such as Reynolds-averaged Navier–Stokes (RANS) solvers, require turbulent combustion models to predict unclosed turbulent terms. This question has challenged the combustion community during the last century. Reacting flows escalate the inherent complexity of turbulence with the introduction of additional mechanisms of vorticity generation and dissipation [11,12]. Turbulence and combustion can promote or inhibit each other, depending on the overlap between their characteristic time and length scales [13,14,15]. Combustion in rocket engines is a particularly challenging case. Due to the high pressures and the absence of nitrogen, the flow’s reactivity is strongly accentuated, dwindling the chemical length and timescales [16]. As a consequence, the Damköhler number is very high, enabling the occurrence of exotic phenomena such as flame-generated turbulence [17,18] and intermittency between combustion in lean and fuel-rich conditions [19]. Knowledge of these phenomena is essential in the formulation of effective turbulent combustion models. Therefore, it is necessary to investigate the conditions in which the combustion of injected propellants develops in a rocket combustor. The present text addresses this question by using direct numerical simulations, representing turbulent mixing and combustion near the injection region of a methane–oxygen engine at high pressure. The resulting databases are post-processed to investigate the occurrence of combustion in different mixing configurations within this environment.

The rest of the paper is organized as follows: The next section provides a theoretical background, which provides necessary concepts concerning turbulence and non-premixed combustion. The third section describes the numerical simulation strategy. The numerical results are analyzed in the fourth section. Finally, the fifth section summarizes the concluding remarks of this research.

2. Theoretical Background

The present section introduces the basic notions that are required to discuss the results of the conducted numerical simulations. The text is structured in three parts. First, the basic concepts of turbulence for variable-density flows are presented. Second, the basic parameters for characterizing turbulent diffusion flames are presented. Third, the classification techniques of local combustion regimes are described, along with their phenomenological interpretation.

2.1. Turbulence

Turbulent flows are characterized by an unstable nature, which precludes reaching steady-state conditions. Due to their chaotic features [20,21], a statistical approach is deemed convenient. For a turbulent quantity $\mathrm{q}$, the most immediate indicator is its Reynolds number, i.e., time average, which is defined as follows:

$$\overline{\mathrm{q}}=\frac{1}{\mathsf{\tau }}{\int }_0^{\mathsf{\tau }}\mathrm{q}\left(\mathrm{t}\right)\mathrm{d}\mathrm{t}$$

(1)

In variable-density flows, such as turbulent flames, it is convenient to introduce density weighting to account for the specific volume variations. The Favre average [22] arises from the following consideration:

$$\tilde{\mathrm{q}}=\frac{1}{\mathsf{\tau }\overline{\mathsf{\rho }}}{\int }_0^{\mathsf{\tau }}\mathrm{q}\left(\mathrm{t}\right)\mathsf{\rho }\left(\mathrm{t}\right)\mathrm{d}\mathrm{t}$$

(2)

The averaged values can be taken as a reference to address the instantaneous values as fluctuations. Therefore, the Reynolds fluctuations are defined as ${\mathrm{q}}^{\prime }\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right)-\overline{\mathrm{q}}$ and the Favre fluctuations as ${\mathrm{q}}^{″}=\mathrm{q}\left(\mathrm{t}\right)-\tilde{\mathrm{q}}$. Higher-order statistics are used to quantify the statistics of these fluctuations. The most relevant indicator is the turbulent kinetic energy, which is essentially the variance of the velocity fluctuations, as follows:

$$\tilde{\mathrm{k}}=\frac{\overline{\mathsf{\rho }{\mathrm{u}}_{\mathrm{i}}^{″}{\mathrm{u}}_{\mathrm{i}}^{″}}}{\overline{\mathsf{\rho }}}$$

(3)

The turbulent kinetic energy is a fundamental manifestation of turbulence, as it indicates the magnitude of the flow unsteadiness and its influence on the momentum conservation equations. Therefore, it denotes the degree of forced convection, which is essential in turbulent combustion and mixing processes. Indeed, turbulent combustion models, such as the eddy dissipation concept (EDC) [23,24] and the flamelet model [25], contain implicit or explicit dependencies on the turbulent kinetic energy [26].

In turbulent premixed flames, turbulent kinetic energy enhances the burning rate, although the flame speed sensitivity to $\tilde{\mathrm{k}}$ decreases with growing turbulence following the so-called bending effect [27]. This is caused by a change in the overlap between chemical and turbulent scales [28,29,30]. If the vortexes are significantly larger than the flame front, their effect on flame development is mainly an increase in the effective surface, due to wrinkling. If the turbulent structures are smaller, the enhancing effects of turbulence primarily occur in terms of diffusion increase. In both cases, it is expected that greater turbulence yields higher flame speeds. Nonetheless, previous research [31] has suggested the possibility of quenching in premixed flames with high turbulent intensity; however, there is no scientific consensus concerning this claim [32]. In diffusion flames, if the chemistry is very fast compared to the turbulence, an enhancement of the burning rate with respect to chemistry takes place [33]. This is primarily motivated by the fact that the mixing intensity grows with turbulent kinetic energy. However, if turbulence becomes significantly faster than the mixing process, it can disrupt the local diffusion flame structure producing local quenching [15].

2.2. Turbulent Diffusion Flames

In a diffusion flame, the oxidizer and the fuel are simultaneously mixing and reacting. The development of these processes is measured with the mixture fraction and the progress variable, which are briefly described in this sub-section.

The mixture fraction $\mathrm{Z}$ is used to assess the local mixing conditions, regardless of the combustion progress. For a chemical mixture consisting of $\mathrm{N}$ species, this parameter is determined as follows:

$$\mathrm{Z}=\sum _{\mathrm{k}=1}^{\mathrm{N}}{\mathrm{Y}}_{\mathrm{k}}{\mathrm{Z}}_{\mathrm{k}}$$

(4)

where ${\mathrm{Y}}_{\mathrm{k}}$ and ${\mathrm{Z}}_{\mathrm{k}}$ stand for the mass fraction and the mixture fraction of the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ chemical species, respectively. The mixture fraction of a molecule is its hydrocarbon content mass fraction. For a generic chemical species with the formulation ${\mathrm{C}}_{{\mathrm{n}}_{\mathrm{C}}}{\mathrm{H}}_{{\mathrm{n}}_{\mathrm{H}}}{\mathrm{O}}_{{\mathrm{n}}_{\mathrm{O}}},$ this parameter is determined as follows:

$${\mathrm{Z}}_{\mathrm{k}}=\frac{{\mathrm{A}}_{\mathrm{C}}{\mathrm{n}}_{\mathrm{C}}+{\mathrm{A}}_{\mathrm{H}}{\mathrm{n}}_{\mathrm{H}}}{{\mathrm{A}}_{\mathrm{C}}{\mathrm{n}}_{\mathrm{C}}+{\mathrm{A}}_{\mathrm{H}}{\mathrm{n}}_{\mathrm{H}}+{\mathrm{A}}_{\mathrm{O}}{\mathrm{n}}_{\mathrm{O}}}\approx \frac{12{\mathrm{n}}_{\mathrm{C}}+{\mathrm{n}}_{\mathrm{H}}}{12{\mathrm{n}}_{\mathrm{C}}+{\mathrm{n}}_{\mathrm{H}}+16{\mathrm{n}}_{\mathrm{O}}}$$

(5)

where ${\mathrm{A}}_{\mathrm{C}}$, ${\mathrm{A}}_{\mathrm{H}}$, and ${\mathrm{A}}_{\mathrm{O}}$ stand for the atomic mass weights of carbon, hydrogen, and oxygen, respectivley. As can be seen, the mixture fraction is zero when only the oxidizer is present, and it tends to unify when the mixture consists only of the fuel. Near the injection site, the mixture fraction profile presents sharp gradients, indicating the mixing of clusters with $\mathrm{Z}\approx 1$ and $\mathrm{Z}\approx 0$. At downstream positions, the $\mathrm{Z}$ field is smoothed, tending to the global mixture fraction ${\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}={\dot{\mathrm{m}}}_{\mathrm{C}\mathrm{H}4}/({\dot{\mathrm{m}}}_{\mathrm{C}{\mathrm{H}}_4}+{\dot{\mathrm{m}}}_{{\mathrm{O}}_2})$.

Besides the mixing of propellants, it is also necessary to evaluate the local development of chemical reactions. The progress variable is used for this purpose. For non-premixed combustion, the progress variable can be determined as follows, as suggested by Bray et al. [34]:

$$\mathrm{c}=\frac{{\mathrm{Y}}_{\mathrm{c}}}{{\mathrm{Y}}_{\mathrm{c},\mathrm{e}\mathrm{q}}\left(\mathrm{Z}\right)}$$

(6)

where ${\mathrm{Y}}_{\mathrm{c}}$ is the reaction progress variable and ${\mathrm{Y}}_{\mathrm{c},\mathrm{e}\mathrm{q}}\left(\mathrm{Z}\right)$ is the equilibrium composition for the present mixture fraction. The reaction progress variable is a scalar that indicates the development of the combustion process. For methane combustion, the definitions devised by Pierce and Moin [35] are often used. However, in space propulsion applications, due to the absence of inert gases, high temperatures take place, with a significant proliferation of free radicals and hydroxyl molecules. For this reason, the mixture fraction definition used is the same in recent works [17,36], as follows:

$${\mathrm{Y}}_{\mathrm{c}}=\frac{{\mathrm{Y}}_{{\mathrm{C}\mathrm{O}}_2}}{{\mathrm{M}}_{\mathrm{C}{\mathrm{O}}_2}}+\frac{{\mathrm{Y}}_{\mathrm{C}\mathrm{O}}}{{\mathrm{M}}_{\mathrm{C}\mathrm{O}}}+\frac{{\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}}{{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}}+\frac{{\mathrm{Y}}_{\mathrm{O}}}{{\mathrm{M}}_{\mathrm{O}}}+\frac{{\mathrm{Y}}_{\mathrm{H}}}{{\mathrm{M}}_{\mathrm{H}}}+\frac{{\mathrm{Y}}_{\mathrm{O}\mathrm{H}}}{{\mathrm{M}}_{\mathrm{O}\mathrm{H}}}$$

(7)

where ${\mathrm{M}}_{\mathrm{l}}$ denotes the molecular weight of the ${\mathrm{l}}^{\mathrm{t}\mathrm{h}}$ chemical species. If the chemical reactions are very fast compared to the turbulent processes, then the progress variable and the mixture fraction are algebraically related [37]. Therefore, the progress variable dynamics are driven by the mixing intensity. This physical phenomenon is evaluated with the scalar dissipation of the mixture fraction, which is defined as follows:

$$\mathsf{\chi }=2{\mathrm{D}}_{\mathrm{m}}{\left(\frac{\partial \mathrm{Z}}{\partial {\mathrm{x}}_{\mathrm{i}}}\right)}^2$$

(8)

where ${\mathrm{D}}_{\mathrm{m}}$ is the mass diffusivity of the reacting mixture. The scalar dissipation rate is analogous to the dissipation rate of the turbulent kinetic energy. Indeed, closure models [25,38,39] relate them as follows:

$$\tilde{\mathsf{\chi }}\propto \frac{\tilde{\mathsf{\epsilon }}}{\tilde{\mathrm{k}}}\widetilde{{\mathrm{Z}}^{″2}}=\frac{\sqrt{\tilde{\mathrm{k}}}}{\tilde{\mathrm{L}}}\widetilde{{\mathrm{Z}}^{″2}}$$

(9)

where $\tilde{\mathrm{L}}={\tilde{\mathrm{k}}}^{\frac{3}{2}}/\tilde{\mathsf{\epsilon }}$ denotes the length scale of the large eddies. In stable non-premixed flamelets with fast chemistry, the heat release rate scales with the scalar dissipation rate at stoichiometric conditions [33] are $\dot{\mathrm{q}}\propto {\mathsf{\chi }}_{st}$. Following (9), it is clear that the local mixing intensity and the burning rate are primarily driven by the turbulent kinetic energy.

In anchored flames with a stable recirculation zone, a standing flame can be maintained, regardless of the injection velocity, since continuous mixing is enabled near the injection posts. Therefore, as the velocity increases, the heat release rate grows as well, due to the greater rate of incoming fresh reactants. This fact can hinder the comparison between flames burning with different injection velocities. Such is the case for the numerical simulations performed in the present work. For this reason, the heat release rate should be normalized with the mean flow rate at the inlet to account for this fact. One possible normalization is as follows:

$$\frac{1}{{\mathsf{\delta }}_{\mathrm{q}}}=\frac{\dot{\mathrm{q}}}{\dot{\mathrm{V}}{\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}\mathrm{E}}$$

(10)

where $\mathrm{E}$ is the specific energy content of methane ($\mathrm{E}\approx 55.5 \mathrm{M}\mathrm{J}/\mathrm{k}\mathrm{g}$) and $\dot{\mathrm{V}}=\mathsf{\rho }{\mathrm{u}}_{\mathrm{i}\mathrm{n}\mathrm{j}}$ is the volumetric flow rate at injection. This value has dimensions of 1/length, indicating the longitude required for the mixture to obtain all of the available energy. Therefore, one suitable reference value to normalize it is the laminar flame thickness at stoichiometric conditions, as follows:

$$\frac{{\mathsf{\delta }}_{\mathrm{L}0}}{{\mathsf{\delta }}_{\mathrm{q}}}=\frac{\dot{\mathrm{q}}{\mathsf{\delta }}_{\mathrm{L}0}}{\dot{\mathrm{V}}{\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}\mathrm{E}}$$

(11)

The most important point of this normalization is the correction with respect to the incoming velocity, which allows the comparison of flows with different injection mass flow rates. More sophisticated normalizations might be possible, but the one presented in this paper suffices within the present work’s scope.

2.3. Flame Classification

For stability and performance purposes, liquid rocket engines typically operate in non-premixed conditions with fuel excess. Despite this global configuration, at local and instantaneous levels, combustion can take place in premixed and lean conditions as well. The occurrence of these processes is dependent on the turbulent mixing of the propellants. Direct numerical simulations allow the detailed observation of these local phenomena. The present sub-section discusses the objective criteria to quantify and classify the local burning regime based on the instantaneous physical fields. The presented methods will be applied to the output of the performed numerical simulations to investigate the occurrence of the different burning conditions.

2.3.1. Fuel-Rich vs. Lean

The classical procedure to calculate the stoichiometric ratio of a hydrocarbon is based on the assumption that the reactants are oxidized into a mixture that is entirely composed of water and carbon dioxide. For the case of methane, this corresponds to the following global reaction scheme:

$$\mathrm{C}{\mathrm{H}}_4+2{\mathrm{O}}_2\to 2{\mathrm{H}}_2\mathrm{O}+\mathrm{C}{\mathrm{O}}_2$$

(12)

which yields a stoichiometric mass O/F ratio of approximately 4, corresponding to $Z_{st}\approx 0.2$. In atmospheric combustion applications that operate in lean conditions, this approximation is quite accurate, since the carbon monoxide content is very low compared to that of carbon dioxide. However, for combustion in space propulsion applications, fuel-rich mixtures are desired in order to maximize the specific impulse. In addition, the high temperatures significantly promote carbon monoxide generation. As a consequence, the equilibrium mass fraction of carbon monoxide significantly outweighs that of carbon dioxide. In an extreme case, the corresponding adjusted global reaction can be written as follows:

$$\mathrm{C}{\mathrm{H}}_4+\frac{3}{2}{\mathrm{O}}_2\to 2{\mathrm{H}}_2\mathrm{O}+\mathrm{C}\mathrm{O}$$

(13)

which implies a stoichiometric O/F of roughly 3, corresponding to ${\mathrm{Z}}_{\mathrm{s}\mathrm{t}}\approx 0.25$. Besides the significant portion of carbon monoxide, additional particularities of combustion in the frame of space propulsion further entangle the equilibrium composition. Due to the high pressures and the absence of nitrogen, very high temperatures are reached in the flame (${\mathrm{T}}_{\mathrm{a}\mathrm{d}}>3400 \mathrm{K}$). Consequently, significant mass fractions of free radicals and hydroxyl are present at equilibrium. All of these substantial deviations from the assumption in (13) demand the redefinition of the stoichiometric condition. A reasonable criterion is to choose the O/F that maximizes the volumetric heat release. The integrated heat release rate is calculated as a function of the mixture fraction using Cantera [40] and a complex chemical mechanism [41]. The result is displayed in Figure 1. As can be seen, the mixture fraction that maximizes heat release takes place at values higher than the classical assumption of ${\mathrm{Z}}_{\mathrm{s}\mathrm{t}}\approx 0.2$. More specifically, the maximum value is found at $\mathrm{Z}\approx 0.25$. This is in good agreement with the previous reasoning, since ${\mathrm{Y}}_{\mathrm{C}\mathrm{O},\mathrm{e}\mathrm{q}}/{\mathrm{Y}}_{{\mathrm{C}\mathrm{O}}_2,\mathrm{e}\mathrm{q}}\approx 4$. Therefore, it can be expected that the redefined stoichiometric value is closer to the one in the assumption in Figure 1. Besides the graph in Figure 1, additional evidence for the better suitability of an $\mathrm{O}/\mathrm{F}\approx 3$ as a stoichiometric value is provided in the discussion of the results.

The shift in the maximum heat release is relevant, since most models for non-premixed combustion require knowledge of the properties at stoichiometric conditions. Within the present work, it is assumed that the stoichiometric mixture fraction in high-pressure methane–oxygen combustion corresponds to ${\mathrm{Z}}_{\mathrm{s}\mathrm{t}}\approx 0.25$. It is important to note that this specific definition is assumed to address the decoupling between the local maximum heat release and the standard stoichiometric oxidizer-to-fuel ratio. From a global standpoint, the classical definition should be considered.

2.3.2. Flame Index

In burning devices that operate in non-premixed conditions, the reacting flow simultaneously displays both premixed and non-premixed configurations at a local level. These variations are driven by the local mixing conditions, which can enable premixed combustion despite the predominant diffusive nature of the flame. Regime markers can be used to identify the local burning regime. The most commonly used marker is the flame index $\mathrm{F}\mathrm{I},$ as defined by Yamashita et al. [42], as follows:

$$\mathrm{F}\mathrm{I}=\frac{1}{2}+\frac{\nabla {\mathrm{Y}}_{{\mathrm{C}\mathrm{H}}_4}\cdot \nabla {\mathrm{Y}}_{{\mathrm{O}}_2}}{2\left|\nabla {\mathrm{Y}}_{{\mathrm{C}\mathrm{H}}_4}\right|\left|\nabla {\mathrm{Y}}_{{\mathrm{O}}_2}\right|}$$

(14)

The original formulation has been re-normalized so that zero corresponds to non-premixed conditions and unity corresponds to premixed conditions. The relationship between this numerical value and the local combustion regime can be easily deduced by the physical meaning of the mathematical expression. If the gradients of the fuel and the oxidizer are aligned, then the local mixing conditions can be deemed premixed. Contrarily, the opposed gradients of the reactants correspond to a counterflow non-premixed configuration. There are several available definitions of a flame index that have been proposed in recent works [43,44,45,46,47,48]. These markers have been extensively reviewed in a recent study [49], where it was found that most definitions yield similar classifications overall. In the frame of the present work, the definition of Yamashita et al. was the one considered, due to its simplicity and extended usage.

An example of the different classifiers presented in this section can be visualized in Figure 2. In this figure, the instantaneous detail of the shear layer for one of the simulated flames is displayed. This graph shows the flame index, mixture fraction, and temperature in order to illustrate the different regimes in which combustion can take place in space propulsion applications. As can be seen, the central region of the shear layer presents the classical counterflow structure of a diffusion flame. Heat release is concentrated in the nearby stoichiometric conditions, where the local O/F maximizes the reacting capabilities of the flow. However, in certain positions far away from the stoichiometric line, the gradients of methane and oxygen are aligned, denoting local premixed combustion. This circumstance is very specific to space propulsion applications. Due to the absence of nitrogen, the flammability limits are significantly widened, enabling combustion in extremely fuel-rich and/or oxygen-rich conditions. Direct numerical simulations enable accessibility to every relevant physical field. Therefore, they offer the possibility to study the occurrence of these events and how they influence the combustion development process.

3. Numerical Simulation Case Setup

The numerical simulations aim to represent the turbulent combustion process that occurs in modern rocket engines in locations near to injection and far away from walls. The turbulent reacting flow is resolved with the compressible solver EBI-DNS [50,51,52,53], which is embedded in the open-source software OpenFOAM [54,55]. This code uses an implicit scheme to resolve the transient conservation equations for mass, momentum, energy, and species with the finite volume method (FVM) [56,57]. EBI-DNS has been applied and validated in several combustion-related problems in recent years [58,59,60,61,62]. A recent study by Zirwes et al. [63] comprehensively assessed the numerical accuracy of this solver for a wide variety of setups, with different reactive and compressibility features. This research shows that EBI-DNS is capable of providing results that are comparable with high-fidelity explicit solvers under a wide range of scenarios. Direct validation of DNS for space propulsion applications is unfeasible, due to the prohibitive computational cost that would be required for a complete system simulation. This difficulty can be addressed by applying convenient normalizations over the results of a reduced engine instance. Following such an approach, a good agreement between experimental and simulated combustion completion points was reported in a previous research paper [64].

The Cantera [40] routines are used for the computation of detailed chemistry and transport properties, using the mixture-averaged transport model described by Kee et al. [65]. The reaction mechanism developed by Slavinskaya et al. [41] is used to determine the chemical source terms of methane combustion using the finite-rate method. This skeletal mechanism consists of 21 species and 97 reactions, and it was devised for combustion in high-pressure space propulsion applications. This mechanism has been applied in several numerical studies of combustion for scale rocket combustors [66,67]. Despite the relatively high complexity of the used reaction scheme, the formation of large aromatic compounds and soot has not been considered. Such a simplification is appropriate for methane combustion, where soot formation is negligible. Nonetheless, it is important to keep in mind the possibility of soot-related phenomena being ignored in the fuel-rich side of the flame. The relevant data concerning combustion are provided in

Table 1. Main turbulence and combustion parameters of the simulated flame.

Symbol	Description	Value
${\mathsf{\delta }}_{\mathrm{L}0}$	Laminar flame thickness at stoichiometric conditions	$5.505 \mathsf{\mu }\mathrm{m}$
${\mathrm{s}}_{\mathrm{L}0}$	Laminar flame speed at stoichiometric conditions	$2.574 \mathrm{m}/\mathrm{s}$
${\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}$	Global mixture fraction	$0.33$
${\mathrm{T}}_{\mathrm{i}}$	Injection temperature	$300 \mathrm{K}$
$p$	Combustion pressure	$20 \mathrm{b}\mathrm{a}\mathrm{r}$

.

The inlet boundary conditions at the injection site in the main simulation are patched using the results of a precursor simulation. In this prior simulation, periodic synthetic turbulence is generated at the inlet via the procedure described in [68]. This scheme is derived from the original formulation proposed by Kraichnan and the subsequent improvements suggested by Shur et al. [69] and Morsbach et al. [70]. These methods are based on discretizing a reference turbulent kinetic energy spectrum in harmonics with random amplitude and phases. These waves are thereafter superposed to provide a synthetic velocity field that complies with the prescribed 2nd order turbulent statistics. The statistics were selected to resemble the output of RANS simulations for a scale methane rocket combustor discussed in a previous work [71]. The described simulation strategy is summarized in Figure 3.

Due to the high computational cost associated with DNS, the hydraulic diameter of the injector ${\mathsf{\varphi }}_{\mathrm{H}}$ is one order of magnitude smaller than the actual setup. As a consequence, the global Reynolds number $\mathrm{R}\mathrm{e}=\mathrm{U}{\mathsf{\varphi }}_{\mathrm{H}}/\mathsf{\nu }$ and the turbulent Reynolds number $\mathrm{R}{\mathrm{e}}_{\mathrm{T}}={\mathrm{k}}^2/\mathsf{\nu }\mathsf{\epsilon }$ are smaller than the reference device. Despite this shift, most of the relevant turbulent and chemical features are preserved. Therefore, the simulations can be assumed to capture the relevant aspects of turbulent mixing and combustion. Relevant turbulent flow parameters at the injector plane are provided in

Table 2. Main turbulence and combustion parameters of the simulated flame.

Symbol	Description	$\mathbf{S}\mathsf{\alpha }$	$\mathbf{S}\mathsf{\delta }$
${\tilde{\mathsf{\eta }}}_{{\mathrm{C}\mathrm{H}}_4}-{\tilde{\mathsf{\eta }}}_{{\mathrm{O}}_2}$	Average Kolmogorov scale	$2.33–2.51 \mathsf{\mu }\mathrm{m}$	$0.517–0.558 \mathsf{\mu }\mathrm{m}$
$\mathrm{R}{\mathrm{e}}_{\mathrm{T},{\mathrm{C}\mathrm{H}}_4}-\mathrm{R}{\mathrm{e}}_{\mathrm{T},{\mathrm{O}}_2}$	Average turbulent Reynolds number	$100–115$	$452–516$
${\mathrm{R}\mathrm{e}}_{{\mathrm{C}\mathrm{H}}_4}-{\mathrm{R}\mathrm{e}}_{{\mathrm{O}}_2}$	Global Reynolds number	$1291–1425$	$5810–6411$
${\mathrm{U}}_{\mathrm{i}}$	Injection velocity	8.3 m/s	37.5 m/s

.

The main simulation comprises a volume of 0.3 × 0.2 × 3 mm, resolved with 144 × 216 × 2160 cubic cells with identical geometry. Periodicity is imposed in the directions perpendicular to the injection. The grid size is small enough to resolve the smallest Kolmogorov eddies in both of the simulations, following the standard requirements [72]. The flame front is resolved with varying resolutions, since the flame front thickness varies. The chosen grid size ensures a resolution of the instantaneous flame front realizations in the order of 10 cells. This is generally enough, although a higher resolution would be desirable. The convergence of the turbulent statistics for the used setup has been addressed in a prior work [36]. The effect of increasing the resolution is marginal. More specifically, a grid size 2.5 times larger yields maximum errors in the order of 5% and axially averaged properties below 1%. A finer grid (halving the cell volume) leads to discrepancies in the axially averaged properties in the order of 0.5%. Therefore, it can be concluded that the resolution of the used mesh is high enough to allow the usage of the simulation results for turbulence research purposes. The timestep size is chosen in an adaptive fashion to ensure that the maximum acoustic CFL number is below 0.9. This corresponds to maximum convective CFL numbers in the order of 0.05 for $\mathrm{S}\mathsf{\alpha }$ and 0.25 for $\mathrm{S}\mathsf{\delta }$. Moreover, due to the implicit nature of the used solver, the fulfillment of the general condition $\mathrm{C}\mathrm{F}\mathrm{L}<1$ is not required to achieve numerical stability. Such low CFL numbers, coupled with the cubic spatial differentiation scheme, ensure the minimization of numerical dissipation, as discussed elsewhere [62].

4. Results Analysis

This section is devoted to the discussion and analysis of the results obtained in the performed numerical simulations. The discussion is structured in three parts. The first and second sub-sections i.e., Section 4.1 and Section 4.2, discuss the overall results of simulations $\mathrm{S}\mathsf{\alpha }$ and $\mathrm{S}\mathsf{\delta },$ respectively. The flame development is discussed from a phenomenological standpoint to identify the main events and the driving physical mechanisms. In the third sub-section (i.e., Section 4.3), the different interactions between turbulent mixing and combustion are discussed in order to investigate the sensitivity of the burning process, with respect to the turbulent kinetic energy.

4.1. Simulation $S\alpha$

The Simulation $\mathrm{S}\mathsf{\alpha }$ is characterized by low inlet turbulence, with a flow behavior close to laminarity. An example of instantaneous flame development can be observed in Figure 4. In this illustration, the flame classification is displayed at the top, and both a positive and a negative heat release rate can be visualized at the center and bottom, respectively.

The flow structure near injection resembles a triple flame, with a non-premixed branch surrounded by lean and fuel-rich premixed branches. Negative heat release can be observed in the transition from fuel-rich premixed to fuel-rich non-premixed conditions. This process is on account of the cracking reactions, which can lead to endothermic results due to the formation of larger hydrocarbons (especially ethylene). The fuel-rich premixed branch displays a discontinuous axial development, with local breakups. This is caused by the forced advection of methane into the flame core. Due to its higher diffusivity, methane benefits more than oxygen from the moderate turbulence. This different sensitivity to turbulence shifts the local O/F in the thick premixed regions outside of the flammability limits. Contrarily, the lean premixed branch displays a more stable configuration. As expected, most of the heat release is obtained from the non-premixed flamelets. Nevertheless, a significant amount of combustion occurs in the lean premixed side, while both of the shear layers are still present. This trend is interrupted by $\mathrm{z}\approx 250 {\mathsf{\delta }}_{\mathrm{L}0}$, when both of the shear layers merge into a single structure. From this point, the combustion at the lean core significantly decreases, and the majority of chemical reactions develop on the non-premixed side.

Figure 5 was elaborated to statistically summarize the observed processes. The right manifold can be used to investigate the probability of occurrence of the different classifications. Most of the observations take place in non-premixed conditions, ranging from moderately lean to significantly fuel rich. This is in good agreement with the example in Figure 4, where it can be observed that the non-premixed flamelet is predominant. There is a significant portion of data in the premixed conditions with extremely lean and fuel-rich mixtures. However, only the lean premixed regions deliver heat release rates in the same order of magnitude as the non-premixed flamelets. Intermediate flame indexes are rather unusual. This is due to the quasi-bimodal distribution of the flame index, which can be observed in the examples in Figure 2 and by other researchers [49]. These graphs provide valuable insights concerning the overall combustion process. However, it is important to remark on their strong dependency on geometry. For example, if the simulation domain were significantly longer, the statistics would shift toward equilibrium conditions. This would imply a greater concentration of points towards $\mathrm{Z}={\mathrm{Z}}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}$, with a lower averaged heat release rate.

4.2. Simulation $S\beta$

An analogous analysis was conducted with the database of Simulation $\mathrm{S}\mathsf{\delta }$. Unlike $\mathrm{S}\mathsf{\alpha }$, this simulation presents a Reynolds number high enough to display fully turbulent features. An example of an instantaneous flame classification and heat release rate is shown in Figure 6.

As can be seen, there is a greater disruption in the fuel-rich premixed region, as this combustion regime only occurs in isolated regions. This is caused by the greater turbulent kinetic energy, which, furthermore, promotes the transport of methane into the flame core. Another remarkable fact is that the negative heat release caused by cracking reactions near to the injection site takes place intermittently. These variations are caused by the flame stretching, which can locally promote or inhibit the heat transfer from the warm products to the cold methane. If heat transfer is promoted, the initial transitory ethylene formation is dampened, due to the local enhancement of the burning rate. The opposite occurs if the flame stretching defocuses the local heat fluxes from the products to the reactants. This effect is discussed in greater detail in the following sub-section. Negative heat release can be observed at downstream positions as well; however, this result is mainly motivated by recombination reactions of the combustion products.

Overall, the absolute heat release rate is substantially larger in $\mathrm{S}\mathsf{\delta }$ compared to $\mathrm{S}\mathsf{\alpha }$. This is mainly a consequence of the greater mass flow rate of the propellants, which enhances the incoming rate of fresh reactants and the turbulent kinetic energy $\tilde{\mathrm{k}}$, which is proportional to the scalar dissipation rate $\tilde{\mathsf{\chi }}$. Nevertheless, the flow in simulation $\mathrm{S}\mathsf{\delta }$ is able to keep up a significant burning rate in the flame core, with lean premixed conditions following the merging of the shear layers. To complete this section, the 2D joint functions are displayed in Figure 7. Few differences can be observed in the 2D joint probability density function. The most remarkable one is the greater prevalence of lean premixed configurations. This is caused by the greater presence of this regime at downstream positions, as can be seen by comparing Figure 4 and Figure 6. Concerning the normalized heat release, it can be seen that it is higher overall, regardless of the considered configuration. However, the enhancements are more remarkable in the lean regions. The differences in the heat release trends are discussed in greater detail in the following sub-section.

4.3. Comparison

Despite the remarkably different flow characteristics, for both of the flames, the chemistry remains very fast compared to the turbulence. To prove this claim, the progress variable, as a function of the mixture fraction, is represented in Figure 8. As can be seen, both of the flames follow the trend expected for fast chemistry [33,37]. Flame $S\delta$ displays a slight deviation, due to its higher turbulent kinetic energy, but the differences are negligible overall. The maximum value of the progress variable is obtained at $\mathrm{Z}\approx 0.25$, which further supports the discussed necessity to consider a shift in the stoichiometric value. The maximum variance of the progress variable conditional to $\mathrm{Z}$ is in the order of $2\times {10}^{-3}$, corresponding to $\mathrm{Z}\approx {\mathrm{Z}}_{\mathrm{s}\mathrm{t}}$, while, for the rest of the domain, $\widetilde{{\mathrm{c}}^{″2}}<{10}^{-4}$ holds. Therefore, it can be assumed that combustion is limited by mixing, since, for a given mixture fraction, the corresponding progress variable is reached almost immediately. The turbulent timescales near injection can be used to examine the validity of the flamelet assumption as well. The mean turbulent timescale in this region is in the order of $50 \mathsf{\mu }\mathrm{s}$ and $10 \mathsf{\mu }\mathrm{s}$ for simulations $S\alpha$ and $S_{\delta }$, respectively. Since the chemical timescale is in the order of 2 $\mathsf{\mu }\mathrm{s}$, it is evident that chemistry is significantly faster in both cases. These points are crucial for the interpretation of the results discussed in this section. Since $\mathrm{c}=\mathrm{f}\left(\mathrm{Z}\right),$ the local heat release rate $\dot{\mathrm{q}}$, which scales with $\dot{\mathrm{c}}$, is essentially dependent on the mixture fraction dynamics. Therefore, the influence that turbulence exerts on the flame development occurs in terms of heat and mass transfer enhancements caused by forced convection.

To further study the influence of turbulence in the mixing process, the probability distributions in Figure 5 and Figure 7 are compared. To assess the sensitivity of these joint distributions, with respect to inlet turbulence, the relative increase is displayed in Figure 9. This graph provides a rough indication of the effects that a turbulence enhancement has with respect to mixing. As can be seen, there is a greater amount of data toward the horizontal center. This is due to improved mixing, which promotes the concentration of data, with respect to the global mixture fraction, i.e., ${\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}\approx 0.33$. However, the probability of presenting extreme conditions, i.e., $(\mathrm{Z}=0,\mathrm{Z}=1$), is enhanced as well. The motivation for this difference is the smaller residence time of the molecules in simulation $\mathrm{S}\mathsf{\delta }$. Since turbulence and mean velocity increase simultaneously, the enhanced mixing is accompanied by a shorter stay time, which increases the possibility of an unmixed cluster reaching the outlet. These effects illustrate that, although these comparative analyses provide valuable information, they should be approached carefully and critically.

To study the variations in heat release, this value was averaged conditionally to the mixture fraction for both of the simulations. The result is displayed in Figure 10. As can be seen, $\mathrm{S}\mathsf{\delta }$ provides a higher heat release overall, with the only exception being the extremely lean regions.

Thanks to the fact that $\mathrm{c}$ and $\mathrm{Z}$ are algebraically related, most of the dependencies with the geometry and residence time are eliminated. Therefore, the differences in Figure 10 are linked to the mixing enhancement driven by turbulence. These effects are strongly dependent on which part of the flame is considered. To simplify the phenomenological description, the different trends are commented on in two different positions: near injection and downstream (where both of the shear layers have merged). The physical driving mechanisms in both of these cases are commented on in the text below.

4.3.1. Near Injection Mixing

To approach the turbulent mixing dynamics near injection, details of the simulated flame are presented in Figure 11. In this figure, the instantaneous stoichiometric line is presented along with the line marking the instantaneous border between the premixed and non-premixed classifications. Therefore, the region between these two isolines corresponds to the lean side of the local diffusion flame. Overall, flame $S\alpha$ exhibits a greater amount of normalized heat release. $S\delta$ presents similar values in the stoichiometric line, but it only keeps up in the stretched regions due to the enhanced convection. Moreover, in these positions, the negative heat release caused by cracking reactions can be suppressed thanks to the high local curvature, which fosters the heat transfer from the stoichiometric line to the fuel-rich zone. Indeed, it can be seen that the regions with non-negative heat release only appear where there is significant stretching.

4.3.2. Downstream Mixing

The flame development dynamics are dramatically altered once the merging of the shear layers takes place. To study the driving mechanisms, details of the local heat release for both simulations are presented in Figure 12. In this illustration, the transparency is varied in order to identify the local burning regime. As can be seen, the lean side of the diffusion flame generates similar heat release in both of the simulations, although, in Simulation $\mathrm{S}\mathsf{\delta }$, the observed values are slightly larger.

The main differences between the simulations take place at the flame core, where the mixture fraction is significantly below its stoichiometric value. In Simulation $\mathrm{S}\mathsf{\alpha }$, intermittency between premixed and non-premixed combustion can be observed. This result points to unstable burnability, due to insufficient mixing intensity. This behavior contrasts with what is observed in $\mathrm{S}\mathsf{\delta }$, where the flame core burns mostly in the premixed conditions. There, the heat release greatly varies between the two flames, since it depends on the local strain and availability of reactants. In both cases, a heat release rate drop can be observed from the lean edge of the diffusion flame to the oxygen core. However, heat release near the diffusive region is significantly larger in $\mathrm{S}\mathsf{\alpha }$, and the values decrease faster. This decreasing trend is motivated by the lower availability of burnable reactants, as the mixture fraction tends to ${\mathrm{Z}}_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{e}}$. Although heat release at the core is smaller in Simulation $\mathrm{S}\mathsf{\alpha }$, the significantly higher burning rate close to the diffusive side enables a greater heat release in lean conditions overall. This process is behind the trend observed in Figure 10 towards the lean side.

The differences in the burning rate between the simulations are primarily driven by the different evolution of the scalar dissipation. To discuss the influence of this variable, instantaneous 2D cuts of its value are displayed in Figure 13.

As expected, the scalar dissipation rate is greater in $\mathrm{S}\mathsf{\delta }$ overall. This is caused by the larger turbulent kinetic energy at the inlet. The most remarkable difference between the simulations corresponds to the behavior following the merging of the shear layers. In Simulation $\mathrm{S}\mathsf{\alpha }$, the mixing intensity decreases after the shear layers collapse. Contrarily, in $\mathrm{S}\mathsf{\delta }$, an increase in scalar dissipation follows the merging of the shear layers. This is motivated by the larger wrinkling of the shear layer, due to the integrated effects of the turbulence departing from the injection site. In $\mathrm{S}\mathsf{\alpha }$, the shear layers are mostly aligned with the injection direction and, consequently, between themselves. Therefore, their merging produces little additional shear, and global dissipation dominates. Nevertheless, in $\mathrm{S}\mathsf{\delta }$, turbulence has enough energy to create substantial oscillations of the diffusive flame structures prior to their merging. In such a scenario, the collapse becomes more abrupt, enhancing the local mixing intensity. This phenomenon enables larger heat release rates, as premixed combustion becomes viable at the core, as displayed in Figure 12 for $\mathrm{S}\mathsf{\delta }$. If the scalar dissipation rate becomes too small, the possibility of combusting in a premixed condition in a stable way is disabled. This is evident in $S\alpha$, but it can occur in localized regions in $\mathrm{S}\mathsf{\delta }$ as well, as can be observed in Figure 6 and Figure 12.

5. Conclusions

Non-premixed methane–oxygen turbulent diffusion flames were studied using direct numerical simulations with finite-rate chemistry. The main purpose of this research was to study the flame development in the frame of modern space propulsion systems. A novel definition of the stoichiometric oxidizer-to-fuel ratio was proposed to capture the effects of significant carbon monoxide formation. The proposed definition is better suited to reaction engines, since the amount of carbon dioxide in equilibrium is small compared to carbon dioxide, altering the main assumption in which the classical definition of stoichiometric conditions is grounded.

The redefined stoichiometric value was used in conjunction with flame regime markers to classify the local combustion regime. Therefore, two main degrees of freedom were considered: fuel-rich vs. oxygen-rich and premixed vs. non-premixed conditions. The statistics of the heat release rate were analyzed conditional to these conditions, with varying Reynolds numbers, to investigate the influence of turbulence in the flame development process. This investigation was performed by a comparative analysis of one mildly turbulent flame ($\mathrm{S}\mathsf{\alpha }$) and one fully turbulent one ($\mathrm{S}\mathsf{\delta }$). For both of the flames, the chemistry was very fast compared to the turbulence. Therefore, combustion progress is limited by mixing, and differences in the heat release dynamics can be primarily attributed to the forced convection variations. It has been found that, despite the eminent non-premixed configuration of the burning device, a significant amount of heat release occurs in lean premixed conditions. This process occurs primarily at the flame core, following the merging of the shear layers. With growing turbulence, this merging becomes more abrupt, due to the increased decorrelation between the two interacting diffusive flame structures.

The results presented here point out several challenges concerning combustion in systems that resort to using oxygen as an oxidizer. Due to the wider flammability limits, combustion in lean premixed conditions can become significant. This effect indicates the necessity to consider a combination of premixed and non-premixed combustion models in turbulent combustion simulations for space propulsion systems. In addition, the definition of the stoichiometric conditions should be redefined to account for the large amount of carbon monoxide in equilibrium.