Analysis of turbulent mixing in a methane–oxygen recessed injector for space propulsion

Martinez-Sanchis, Daniel; Sternin, Andrej; Jocher, Agnes; Haidn, Oskar

doi:10.1063/5.0160320

Turbulent mixing in a methane–oxygen recessed injector is studied using direct numerical simulations. The operating point is chosen to be fuel-rich and at high pressure to recreate a representative environment for space propulsion applications. The results are used to investigate the transport of the turbulent mixture fraction statistics and the validity of conventional transport models. It is observed that molecular diffusion is only relevant near the boundary layer of the injection recess cavity and at the recirculation zone. Moreover, turbulent mixing in the axial direction is negligible as radial turbulent diffusion dominates. Radial turbulent diffusion near injection is driven by Kelvin–Helmholtz Instabilities (KHI) manifesting at large scales in the order of the injector geometry. The dominance of this process over microscale mixing originates negative turbulent diffusion, which produces a mixture resegregation and the appearance of lean pockets far from the oxidizer injection plane. Gradient models display poor capabilities for the prediction of this sort of phenomena. Closure models for the turbulent mixing transport terms are proposed and evaluated. An anisotropic gradient model is devised, providing performance improvements within the recess cavity and the recirculation region. In addition, a novel filtered Reynolds-averaged Navier Stokes approach based on the mixing state is proposed. This new methodology shows excellent prediction capabilities in the regions dominated by KHI, accurately predicting negative turbulent diffusivity. The challenges associated with this model are commented on, and strategies to enable its application are proposed.

I. INTRODUCTION

Methane will play a key role in the future generation of liquid rocket engines.^1–3 This propellant has excellent storability and handling properties, coupled with a high specific impulse.⁴ These qualities make it an ideal substitute for hydrazine in upper-stage applications.⁵ In addition, when compared to kerosene, methane displays a higher effective exhaust velocity, simpler processing, and more stable combustion.^4,6 Hence, methane has additionally established itself as an attractive option to replace kerosene in first-stage engines.⁷ Despite its numerous advantages, fundamental questions concerning ignition, mixing, and stable combustion require further investigation.⁶ This necessity has motivated a substantial number of research works focused on methane–oxygen combustion in recent years.^8–12

Combustion chambers of modern rocket engines burn propellants in a non-premixed regime. Fuel and oxidizer are injected separately, with their burning limited by the mixing rate. Hence, fast mixing is essential to complete chemical reactions within the minimum volume to reduce the thrust chamber size. The coupling between mixing and combustion development is especially evident in the frame of the flamelet model,^13–15 where infinitely fast chemistry is assumed. In such a case, the development of chemistry is strictly limited by the mixing of the reactants. This is usually the case for liquid rocket engines, where high pressures and the absence of inert gases yield very small chemical timescales.^16–18

Mixing of the propellants is also relevant from a mechanical standpoint. Cyclic oxidation and reduction of copper at the combustion chamber wall generates a particular material degradation process named blanching.^19–21 Blanched surfaces display a sponge-like appearance with a significant loss of thermal conductivity and increased surface roughness.²² These physical changes yield higher temperatures at the combustion chamber wall, with a subsequent decrease in the liner life.²³ This effect poses a significant challenge for the design of new engines, especially those conceived to be reusable. Rocket engines burning light hydrocarbons or hydrogen typically operate in fuel-rich configurations, with oxygen injected in the inner side of the coaxial injector. Although such a configuration should preclude oxygen from reaching the walls, blanching has been observed in several engines where oxygen is the limiting reactant. In hydrogen–oxygen engines, the Space Shuttle Main Engine (SSME)^20,22 with an equivalence ratio $ϕ \approx 1.31$ ⁠, HM60^24,25 (⁠ $ϕ \approx 1.58$ ⁠), and Vulcain²⁶ (⁠ $ϕ \approx 1.49$ ⁠) are prominent examples of combustion devices that exhibited blanching at some point. Blanching has also been reported in fuel-rich methane–oxygen rocket combustors with $ϕ \approx 1.15 .$ ²⁷ The most common strategies to eliminate blanching are injection trimming and improved cooling.²⁶ Although such measures effectively prevent this undesired phenomenon, they do not address the fundamental problem, which is the lack of a complete understanding of the turbulent mixing process. Moreover, blanching countermeasures can be detrimental to the overall system performance. Injection trimming increases the propellant supply system complexity and creates a more heterogeneous mixture with the subsequent combustion performance downgrade. Cooling improvements typically involve either increasing the turbomachinery power requirements or thinning the combustion chamber walls. The former option reduces the amount of energy available for thrust generation, whereas the latter drives up manufacturing costs.

For all the stated reasons, an improved understanding of the turbulent mixing process is fundamental to supporting future combustion chamber design for space propulsion applications. Recent research is mainly focused on the mixture fraction probability density function (PDF) modeling.^28–30 However, modeling the mixture fraction transport in the frame of Reynolds Average Navier–Stokes (RANS) has received little attention, although this is the most common approach followed industrially. In most extended solvers, several transport sources are often neglected or modeled with simple gradient assumptions.^31,32 Nonetheless, counter-gradient turbulent diffusion has been observed both experimentally³³ and in numerical simulations,³⁴ directly challenging the most fundamental assumption of commonly used models. Understanding the physical mechanisms behind these processes is needed to devise improved numerical models and understand the performance drivers in the turbulent mixing process.

The present paper addresses the problem of RANS turbulent closure for mixture fraction statistics in the frame of non-premixed flames for space propulsion applications. Turbulent mixing and combustion in a shear coaxial methane–oxygen injector is studied using direct numerical simulations (DNS) with a complex chemical reaction mechanism. The transport budget of the mixture fraction statistics is evaluated to understand the main drivers in the mixing process. The performance of conventional models is assessed, and new strategies are proposed to circumvent the deficiencies of standard approaches.

The rest of the paper is structured as follows. A fundamental theoretical background is laid in Sec. II, where basic concepts concerning turbulence and turbulent mixing are introduced. The simulation setup is detailed in Sec. III. The results of the DNS database are analyzed in Sec. IV to investigate the drivers influencing. Finally, the concluding remarks of this research are provided in Sec. V.

II. THEORETICAL BACKGROUND

In laminar flows, mixing is primarily driven by molecular diffusion, which accounts for the tendency of heterogeneous flows to homogenize due to Brownian molecular motion. This mixing mechanism is relevant to a wide number of applications but of limited magnitude. To achieve high thrust-to-weight ratios, it is necessary to devise small combustion chambers which completely mix and burn the propellants within minimum volumes. Turbulence is necessary to accelerate these processes. Thanks to the chaotic flow motion, heat, and mass transfer phenomena are enhanced, enabling significant improvements in terms of vehicle performance. However, turbulent flows present a remarkable analytical complexity due to their unsteady nature. The present section discusses the fundamental aspects of turbulent flows with a special focus on the mixing processes. First, basic notions of turbulent flows with variable density are introduced. Second, the transport equations of the turbulent mixture fraction statistics are discussed. In Sec. II C, the gradient assumption for closure is discussed along with its physical implications. Finally, a mathematical framework to filter flow statistics based on the mixing configuration is introduced.

A. Turbulence characterization

Turbulent flows display a chaotic behavior.^35,36 Hence, the relevant physical quantities are typically approached from a statistical standpoint instead of deterministically. Let us consider a turbulent quantity

q

⁠, whose value fluctuates over time. The most immediate statistic is its time average, which is generally described as

\bar{q} = \frac{1}{τ} \int_{0}^{τ} ρ (t) q (t) d t .

(1)

This value is typically referred to as the Reynolds average. Reynolds averages are useful to study any turbulent quantity in constant-density flows. However, if significant density variations are foreseen, it is convenient to weigh the non-inertial quantities with density. This motivation is behind the Favre-average,³⁷ which is mathematically expressed as

\tilde{q} = \frac{1}{\bar{ρ} τ} \int_{0}^{τ} ρ (t) q (t) d t .

(2)

Taking average values as reference, it is possible to define Reynolds and Favre fluctuations as

q^{'} (t) = q (t) - \bar{q}

and

q^{″} (t) = q (t) - \tilde{q}

⁠, respectively. The statistical behavior of the fluctuations is relevant due to its coupling with the transport of average properties. Since the average value of fluctuations is zero, higher-order statistics are typically used to characterize the behavior of fluctuations. The most general second-order statistic is the covariance, which indicates the correlation between the fluctuations of two turbulent scalars

q

and

r

⁠,

\tilde{q ″ r ″} = \frac{1}{\bar{ρ} τ} \int_{0}^{τ} ρ (t) q ″ (t) r (t) ″ dt .

(3)

In turbulent flows, the Reynolds stress tensor (RST)

\tilde{u_{i}^{″} u_{j}^{″}}

containing the covariances of velocity fluctuations plays a fundamental role since it is directly related to the transport of momentum. Hence, an accurate prediction of the velocity field and all its related parameters depends on the RST. The magnitude of this tensor is usually summarized with the turbulent kinetic energy, which represents the extent of the velocity variance: This tensor is often simplified with the turbulent kinetic energy, which represents the magnitude of velocity variances,

\tilde{k} = \frac{1}{2} \tilde{u_{k} u_{k}} .

(4)

The turbulent kinetic energy constitutes an integral quantity, which arises from the superposed effect of multiple vortices with varying sizes and energy. The specific length scales at which turbulence manifest are crucial to characterize this phenomenon. A characteristic large eddy size may be defined as³⁸

\tilde{L} = \sqrt{{\tilde{k}}^{3}} / \tilde{ε},

(5)

where

\tilde{ε}

denotes the Favre-averaged dissipation rate of the turbulent kinetic energy, which represents the dissipation of turbulent fluctuations due to molecular viscosity.

This effect causes a successive breakup of vortical structures. The fragments of originally large eddies reach a size and are fragmented into smaller entities; they are small enough to be so that viscosity dominates. The size of the smallest eddies is

\tilde{η} = {(\frac{{(\tilde{μ} / \bar{ρ})}^{3}}{\tilde{ε}})}^{1 / 4} = {(\frac{{\tilde{ν}}^{3}}{\tilde{ε}})}^{1 / 4} .

(6)

The discussed eddy sizes are particularly relevant in the frame of scale-resolving simulations since they are tightly coupled with the resolution requirements. The Kolmogorov scale determines the maximum cell grid size that should be chosen for a simulation to be labeled as DNS, with

Δ x_{max} \approx 2.1 η .

³⁹ Larger cells might require sub-grid scale models to account for the non-resolved turbulent kinetic energy. Large and smallest eddies are related through the turbulent Reynolds number,

R e_{T} = \frac{{\tilde{k}}^{2}}{\tilde{ε} ν} = {(\frac{\tilde{L}}{\tilde{η}})}^{4 / 3},

(7)

which is an indicator of the turbulence bandwidth. Hence, as the turbulent Reynolds number grows, the range of existing eddy sizes widens. The turbulent Reynolds number assesses the magnitude of turbulence at a specific flow location. In canonical flows, it is convenient to define a global Reynolds number that characterizes the overall flow turbulence. For channel flows, the following definition is commonly adopted:

Re = \frac{ϕ_{H} U}{ν},

(8)

where

ϕ_{H}

denotes the hydraulic diameter and

U

the bulk flow velocity. For a specific setup, turbulent and global Reynolds numbers are proportional,

R e_{T} \propto Re,

due to the similarity between the scales of the large eddies and the channel geometry, i.e.,

\tilde{L} \propto ϕ_{H}

⁠, and

\sqrt{\tilde{k}} \propto U

⁠.

B. Turbulent mixing

In non-premixed flames, fuel and oxidizer are injected separately, and mixing is required for combustion to develop. The mixing state is evaluated with the mixture fraction, which essentially denotes the relative weight of hydrocarbon content. The mixture fraction of the

k^{th}

molecular species, with a composition

C_{n_{C}} H_{n_{H}} O_{n_{O}}

is expressed as

Z_{k} = \frac{n_{C} M_{C} + n_{H} M_{H}}{n_{C} M_{C} + n_{H} M_{H} + n_{O} M_{O}},

(9)

where

M_{C}

⁠,

M_{H}

⁠, and

M_{O}

denote the atomic weights of carbon, hydrogen, and oxygen, respectively. More refined approaches⁴⁰ are convenient if the fuel contains oxygen (e.g., alcohols). However, for the case of light hydrocarbons the definition in suffices. Within a given mixture containing

N

different chemical species, each with their corresponding mass fraction

Y_{k}

⁠, the overall mixture fraction is determined as

Z = \sum_{k = 1}^{N} Y_{k} Z_{k} .

(10)

For methane–oxygen flames, the stoichiometric mixture fraction corresponds to

Z_{st} \approx 0.2

⁠, with higher values denoting fuel-rich concentrations and smaller lean configurations. The mixture fraction is a conserved scalar in the sense that it remains unaffected by chemical reactions since these alter molecules but not atoms. In turbulent flows, this parameter displays a fluctuating behavior due to its coupling with the momentum conservation equations that are characterized by the mentioned chaotic behavior. Hence, the mixture fraction shall be analyzed through its statistics as any turbulent quantity. We first consider the Favre-average mixture fraction

\tilde{Z}

⁠, whose distribution over space and time is governed by the following transport equation:

\frac{\partial \bar{ρ} \tilde{Z}}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} \tilde{Z})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (\bar{ρ D \frac{\partial Z}{\partial x_{i}}} - \bar{ρ u_{i}^{″} Z ″}) = D_{Z, M_{i}} + D_{Z, T_{i}},

(11)

where

D_{Z, M_{i}}

and

D_{Z, T_{i}}

stand for the molecular and turbulent diffusive transport. These terms are unclosed, and modeling is required for their determination. The molecular diffusion term is usually neglected in flows with sufficiently high Reynolds numbers, since turbulent diffusion is the dominant mixing mechanism in such a context. The transport term associated with turbulent diffusion is usually modeled following the gradient assumption, which is detailed in Sec. II C. In addition to average values, a more detailed statistical description of the mixture fraction is usually needed to close certain combustion models. For this reason, modeling the transport of the mixture fraction variance

\tilde{Z^{″ 2}}

is also of great interest in several models, such as the flamelet approach. The formalisms involving this closure problem are not detailed in this text for the sake of brevity. However, it is worth noting that it presents important similarities with the challenge of closing .

C. Gradient assumption

The gradient assumption plays a pivotal role in most RANS solvers since it is the most extended approach for closing unresolved turbulent fluxes. In its general form, an unresolved turbulent flux

\bar{ρ u_{i} ″ q ″}

is closed as

\bar{ρ u_{i} ″ q ″} = - \bar{ρ} D_{T} \frac{\partial \tilde{q}}{\partial x_{i}},

(12)

where

D_{T}

is the turbulent diffusivity. Using this expression, the unclosed turbulent fluxes in and can be modeled as

\bar{ρ u_{i} ″ Z ″} = - \bar{ρ} D_{T} \frac{\partial \tilde{Z}}{\partial x_{i}},

(13)

and

\bar{ρ u_{i}^{″} Z^{″ 2}} = - \bar{ρ} D_{T} \frac{\partial \tilde{Z^{″ 2}}}{\partial x_{i}} .

(14)

The turbulent diffusivity

D_{T}

is usually calculated from the turbulent viscosity, which is modeled with the following expression:⁴¹^,

ν_{T} = C_{μ} \frac{\tilde{k}}{\tilde{ε}},

(15)

where

C_{μ}

is a model constant, with a recommended value

C_{μ} \approx 0.09 .

^42,43 The turbulent diffusivity may be retrieved as

D_{T} = ν_{T} / {Sc}_{T}

⁠, where the turbulent Schmidt number

{Sc}_{T}

is typically assumed to be a constant on the order of unity. If the turbulent viscosity is inserted into the turbulent Reynolds number definition, it is possible to rewrite it as

{Re}_{T} = \frac{ν_{T}}{C_{μ} \tilde{ν}},

(16)

which illustrates that the Reynolds number essentially indicates the ratio between turbulent and molecular viscosity. We can take advantage of this definition to relate turbulent and molecular diffusivities,

D_{T} = \tilde{D} \frac{{Re}_{T} S c C_{μ}}{{Sc}_{T}} .

(17)

If we insert this definition in the transport equation of the average mixture fraction, neglecting the correlation between molecular diffusivity and mixture fraction gradients, the following expression is derived:

\frac{\partial {\tilde{u}}_{i} \tilde{Z}}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (\frac{\partial \tilde{Z}}{\partial x_{i}} \tilde{D} (1 + {Re}_{T} C_{μ} \frac{Sc}{{Sc}_{T}})) \approx \frac{\partial}{\partial x_{i}} (\frac{\partial \tilde{Z}}{\partial x_{i}} \tilde{D} (1 + {Re}_{T} C_{μ})),

(18)

where the time derivative term was erased since it tends to zero in statistically stationary flows. This simplified formula serves to illustrate the role of turbulence in the mixing process. High diffusivity (molecular or turbulent) enables the completion of the mixing process using minimum length. In turbulent flows

R e_{T} C_{μ} ≫ 1

⁠, and turbulence is the predominant mixing mechanism. This simplified expression can also be used to deduce the consequences of exotic phenomena such as backscattering. Reversed energy transport has been observed in turbulent reacting flows in several previous studies for both premixed^44,45 and non-premixed flames.⁴⁶ In practice, the effect of backscatter is the counter-diffusion transport of turbulent statistics, which manifests as negative turbulent viscosity and diffusivity. From a formal standpoint, such a scenario can be interpreted as

C_{μ}

not being constant and exhibiting negative values, effectively leading to negative turbulent diffusion. Such circumstances lead to counterintuitive processes in which the average concentrations of fuel and oxidizer tend to resegregate instead of further mixing. These effects play a fundamental role in the overall turbulent mixing dynamics, critically influencing the performance of injection systems. Due to the high-pressures and subsequent fast chemistry, such a scenario is likely to occur in turbulent non-premixed flames for space propulsion applications. However, backscattering is typically neglected in most RANS solvers, where

C_{μ}

is assumed to have an invariant value.

D. Mixing-based filtering

The discussion of new models for the unclosed transport terms requires the introduction of a novel filtering methodology in the averaged turbulent quantities. The main idea behind this approach is to address the conservation equations at lean and fuel-rich conditions to decouple the fluctuations between these modes from ones within the diffusion flame front. The necessary mathematical framework is introduced in this subsection.

We first begin defining the lean

(λ)

and fuel-rich

(φ)

filters as

λ (Z) = \{\begin{array}{l} 1, & i f Z \leq Z_{st} \\ 0, & i f Z > Z_{st}, \end{array}

(19)

and

φ (Z) = \{\begin{array}{l} 1, & i f Z \leq Z_{st} \\ 0, & i f Z > Z_{st} . \end{array}

(20)

With the following properties holding at any time,

\{\begin{array}{l} λ + φ = 1 \\ λ φ = 0 \end{array} \forall t, \overset{⇀}{x} .

(21)

These features are very advantageous from an operational standpoint. They enable a seamless decomposition of the conservation equations of average properties without the introduction of cross correlation terms. The Reynolds average of a turbulent quantity

q

conditioned to a given mixing conditions can be expressed as

{\bar{q}}_{θ} = \frac{\bar{q θ}}{\bar{θ}},

(22)

where

θ

denotes

λ

or

φ

indistictively. Closure between global averages and the filtered ones requires the determination of the Favre-averages of each filter, which can be generally expressed as

\tilde{θ} = \frac{\bar{ρ θ}}{\bar{ρ}} = \int_{0}^{1} θ (Z) \tilde{f} (Z) dZ .

(23)

The Favre-averages of the filters

\tilde{λ}

and

\tilde{φ}

are unclosed, and modeling is required for their determination. In the Annex of the present text, a methodology to derive these quantities is proposed, applied and validated. Once the Favre-average values of the filters are known, their Reynolds average can be derived as

\{\begin{array}{l} \bar{λ} = \frac{\tilde{λ} {\bar{ρ}}_{φ}}{\tilde{λ} {\bar{ρ}}_{φ} + \tilde{φ} {\bar{ρ}}_{λ}} \\ \bar{φ} = \frac{\tilde{φ} {\bar{ρ}}_{λ}}{\tilde{λ} {\bar{ρ}}_{φ} + \tilde{φ} {\bar{ρ}}_{λ}} . \end{array}

(24)

A detailed demonstration of this result is provided in the annex. If the average filter values are known, global statistics can be explicitly obtained from filtered ones. For example, Reynolds averages can be obtained as

\bar{q} = \bar{q (λ + φ)} = \bar{q λ} + \bar{q φ} = \bar{λ} {\bar{q}}_{λ} + \bar{φ} {\bar{q}}_{φ},

(25)

and Favre-averages

\tilde{q} = \frac{\bar{ρ q}}{\bar{ρ}} = \frac{\bar{ρ q λ} + \bar{ρ q φ}}{\bar{ρ}} = \frac{\bar{ρ λ} {\tilde{q}}_{λ} + \bar{ρ φ} {\tilde{q}}_{φ}}{\bar{ρ}} = \tilde{λ} {\tilde{q}}_{λ} + \tilde{φ} {\tilde{q}}_{φ} .

(26)

Second-order statistics can be easily retrieved as well. A global covariance between two scalars may be expressed as

\tilde{q ″ r ″} = \tilde{λ q ″ r ″} + \tilde{φ q ″ r ″} .

(27)

Global fluctuations can be expressed as a function of the filtered ones,

q^{″} = q - \tilde{q} = q_{θ}^{″} + {\tilde{q}}_{θ} - \tilde{q} = q_{θ}^{″} + Δ {\tilde{q}}_{θ},

(28)

where

Δ {\tilde{q}}_{θ}

denotes the offset between global and filtered averages, i.e.,

Δ {\tilde{q}}_{θ} = {\tilde{q}}_{θ} - \tilde{q}

⁠. Resorting to this notation and applying some fundamental relationships, a global covariance can be expressed,

\tilde{q ″ r ″} = \frac{\bar{λ} {\bar{ρ}}_{λ} (\tilde{q_{λ}^{″} r_{λ}^{″}} + Δ {\tilde{q}}_{λ} Δ {\tilde{r}}_{λ}) + \bar{φ} {\bar{ρ}}_{φ} (\tilde{q_{φ}^{″} r_{φ}^{″}} + Δ {\tilde{q}}_{φ} Δ {\tilde{r}}_{φ})}{\bar{λ} {\bar{ρ}}_{λ} + \bar{φ} {\bar{ρ}}_{φ}},

(29)

where only the filtered covariances

\tilde{q_{λ}^{″} r_{λ}^{″}}

require modeling. Every other quantity within is closed in the frame of the proposed approach. The unclosed covariances may be retrieved using conventional models such as the ones derived from the gradient assumption. Under these considerations, the closure of non-resolved turbulent statistics does not present any additional formal challenge compared to the traditional RANS approach. Third- and higher-order statistics can be derived following analogous procedures.

III. SIMULATION SETUP

The numerical simulation setup is described in this section. The discussion is structured in two main parts. First, the simulation strategy is presented, discussing the boundary conditions along with the characteristics of the operating regime. In Sec. III B, the DNS environment is detailed, and the compliance with resolution requirements is thoroughly assessed.

A. Simulation strategy

The performed numerical simulations can be observed as two different subsets. First, the unmixed methane and oxygen flows develop into physically meaningful turbulence in a precursor simulation. This simulation can be considered as the supply manifolds transporting the propellants to the combustion chamber. Hence, the resulting velocity fields are subsequently used to feed the turbulent inlet of the main simulation, which is where turbulent mixing and combustion take place. At the inlet of the precursor simulation, periodic synthetic turbulence is enforced, following the procedure detailed in Ref. 47, which derives from the works of Morsbach and Franke⁴⁸ and Shur et al.⁴⁹ These different methodologies share the common principle of superimposing a set of randomized harmonics derived from the discretization of a reference spectrum, as originally proposed by Kraichnan.⁵⁰ The synthetically generated velocity field develops with enough length to ensure convergence into a mature turbulent flow as proved in a previous study.⁴⁷ The synthetically generated turbulent field at the flow core is enforced with $k = 10 m^{2} / s^{2}$ for both propellants, and the dissipation rate was chosen to fulfill the condition $L = ϕ_{inj} / 6$ ⁠. This allows to define the spectral kinetic energy function as detailed bay Pope.³⁸

The main simulation is conceived with a splitter-plate configuration, with periodic boundary conditions in the spanwise direction. The remaining two directions are denoted as radial (x) and axial (z), as illustrated in Fig. 1. Although such naming is not entirely accurate, it simplifies the discussion and facilitates analogies with shear coaxial injectors. The hydraulic diameter of the entire injector element $ϕ_{H}$ is taken as a reference length scale as displayed in Fig. 1. Gaseous methane and oxygen are initially injected into the recess cavity, which spans axially over one hydraulic diameter preceding the combustion chamber. Within this region, the first stages of turbulent mixing and combustion take place, with the formation of two symmetrical diffusion flames anchored at the two injection posts. The partially burnt flow is second injected into the last simulation, which is referred to as the combustion chamber, and it spans axially over ten hydraulic diameters, wherein radial periodic boundary conditions are enforced. These boundary conditions aim to recreate a representative scenario in which a large number of injector elements are installed, and the influence of the combustion chamber wall on the central ones is negligible. Such is the case in most modern injection systems, which can feature hundreds of individual injection elements. The described simulation setup is summarized in the schematic displayed in Fig. 1.

FIG. 1.

View large Download slide

Overview of the simulation strategy in a recessed methane–oxygen injector.

Both the recess cavity and the combustion chamber section constitute prismatic volumes. The recess region is characterized by dimensions of 0.3 × 0. $2$ × $0.3$ mm³, and it is resolved with 240 × 160 × 240 uniform cubic cells. The combustion chamber segment comprises a volume of $0.6$ × $0.2$ × $3 mm$ ³, and it is resolved with 480 × 160 × 2400 cubic cells whose geometry is identical to those in the recess segment. The spanwise length is limited by the computational cost. It would be desirable to use a larger configuration to ensure the absence of pollution due to the periodic boundary condition. However, large eddy simulations of this configuration show marginal changes in first and second-order averages if substantially larger spanwise lengths are used. In addition, the analysis of the autocorrelation function in the spanwise direction shows that correlation. The dimensions and operating regime were chosen to resemble those of a family of scale rocket combustors, which has been intensively studied over the last decade.^8,27,51,52 However, relevant differences motivated by various constraints should be remarked. Compared to the reference system, the dimensions of the numerical simulation are one order of magnitude smaller, and the mean bulk injection velocity is halved. Consequently, the global Reynolds number in the numerical simulations is roughly 20 times smaller than the reference one. More specifically, methane and oxygen are injected with identical bulk velocities corresponding to 37.3 m/s, which provides global Reynolds numbers of 2113 and 4672 for methane and oxygen, respectively. Within the recess segment, the Reynolds number decreases from 6943 to 3612 due to the viscosity increase as combustion progresses. The mentioned reductions in size and velocity arise from the existing limitations in computational cost. However, chemical timescales are preserved since the propellant combination, and the combustion chamber pressure are identical to the ones in space propulsion applications. The geometrical changes primarily influence the large vortical scales. Since the most relevant turbulence–chemistry interactions occur at small scales,^53,54 it can be concluded that the used geometry recreates the most consequential phenomena. Detailed parameters concerning the chemical flow characteristics are provided in Table I.

TABLE I.

Main combustion parameters of the simulated flame.

P

Pressure

20 bar

T_{u}

Temperature of the unburnt reactants

300 K

(O / F)

Global oxidizer-to-fuel ratio

1.99

{(O / F)}_{st}

Stoichiometric oxidizer-to-fuel ratio

3.9898

ϕ_{CH 4}

Hydraulic diameter of methane injector

0.05 mm

ϕ_{O_{2}}

Hydraulic diameter of oxygen injector

0.1 mm

ϕ_{H}

Hydraulic diameter of the entire injection element

0.3 mm

δ_{L 0}

Laminar flame thickness at stoichiometric conditions

5.505 μ m

s_{L 0}

Laminar flame speed at stoichiometric conditions

2.5735 m / s

An example of the transient fields obtained in the performed simulations is displayed in Fig. 2 (multimedia available online). In this illustration, an instantaneous cross section of the mixture fraction and the carbon monoxide mass fraction is represented. As can be seen, the domain size enables the observation of several relevant phenomena, such as flame anchoring at the injection post, initial mixing, or recirculation near the faceplate. However, it is important to remark that the domain size is not large enough to enable combustion completion, i.e., reaching a quasi-steady state chemical equilibrium state. The authors estimate that the domain axial length should be $z \in [20 - 40] mm$ to achieve combustion completeness with the current injection conditions.

FIG. 2.

$Example of instantaneous fields in a 2D cross section: Mixture fraction Z (a) and carbon monoxide mass fraction Y C O (b). Multimedia available online.$

View large Download slide

Example of instantaneous fields in a 2D cross section: Mixture fraction $Z$ (a) and carbon monoxide mass fraction $Y_{C O}$ (b). Multimedia available online.

B. Numerical solver and resolution requirements

The main simulations were conducted with the reactive compressible solver EBI-DNS,^55,56 which is coded within the open-source environment OpenFOAM.^57,58 The transient equations for the conservation of mass, momentum, energy, and species are resolved implicitly using the Finite Volume Method (FVM).^59,60 This solver has been applied and validated over a wide variety of combustion-related problems.^61–63 In a recent publication,⁶⁴ it was found that the fidelity of EBI-DNS is paired with well-established reactive DNS codes, while displaying excellent scalability properties.

Detailed chemistry and transport properties are computed with Cantera,⁶⁵ using the mixture-averaged transport model as described by Kee et al.⁶⁶ The skeletal mechanism developed by Slavinskaya et al.¹⁶ is applied to determine the reaction rates of the methane combustion process using the finite rate approach. This chemical mechanism was devised for high-pressure combustion in space propulsion applications, and its uncertainty has been thoroughly examined for this sort of environment.

The main challenge associated with DNS is their high computational cost due to the requirement of resolving all relevant physical processes. In turbulent reacting flows, the cell size should be small enough to ensure an appropriate resolution of the fluid dynamic and chemical phenomena. The resolution quality is objectively assessed with the number of cells used to represent the Kolmogorov eddies and the flame front.

The reference recommendation for the Kolmogorov eddies is to use a minimum of 1/2.1 cells.³⁹ The Kolmogorov scale is a turbulent parameter that varies over the simulation domain due to its dependency on viscosity and the dissipation rate of turbulent kinetic energy. Therefore, the resolution quality varies over space accordingly. The most comprehensive way to assess the resolution of the smallest vortexes is to analyze the spatial field of the Kolmogorov length scale and compare it with the local cell size. This analysis is contained in Fig. 3. As can be seen, the most critical regions are the recess cavity and the injection faceplate. This is motivated by the small degree of combustion development. In these regions, the flow density remains high, and the viscosity is low, yielding greater Reynolds numbers and smaller Kolmogorov eddies. The resolution in specific critical axial locations can be observed in Fig. 3(b). As can be seen, the resolution does not comply with the recommended requirements at certain locations at the inlet plane and at the faceplate. This marginally insufficient resolution near the wall may yield a slight underprediction of turbulent production. However, the cell size remains very close to the recommended value even in locations with the poorest resolution, and it may be safely concluded that such effects should remain small. Moreover, near-wall turbulence is not the primary goal of this study. The Kolmogorov eddies are very well resolved at the locations where turbulent mixing and combustion take place, which are the focus of the present research. In a recent work with the same computational environment, Zirwes et al.⁶⁷ achieved an excellent agreement with experimental results even though the smallest eddies were slightly under resolved near injection. For all the stated arguments, it can be concluded that the chosen grid size complies with the resolution requirements from the standpoint of fluid mechanics.

FIG. 3.

View large Download slide

Resolution of the Kolmogorov eddies: Number of cells used to resolve the Kolmogorov eddies through the entire domain (a) and radial-wise resolution of the Kolmogorov eddies at specific axial positions (b).

The flame front resolution constitutes the main requirement from the standpoint of chemistry. In the context of premixed flames, 10–20 cells per laminar thermal flame front thickness are usually recommended.^68,69 However, in diffusion flames the realized flame front thickness varies over space and time due to the fluctuating mixing conditions. In methane–oxygen flames at the chosen operation conditions, it has been observed that the flame front thickness is typically 2–4 times larger than the one at stoichiometric conditions.¹² A local indication of the flame front thickness may be derived from the thermal gradients,

δ_{th} = \frac{Δ T}{‖ \nabla T ‖},

(30)

where

Δ T

indicates the offset between the temperature of the reactants and the products, which corresponds to

Δ T \approx 3150 K

for the chosen operation regime. With this definition, it is possible to evaluate the instantaneous realized flame thickness. Since

δ_{t h}

varies over space and time, this parameter should be approached statistically. Figure 4(a) was elaborated to summarize the realized flame front resolution. In this illustration, the probability density functions (PDFs) of the instantaneous flame front resolution are displayed. As can be seen, the overwhelming majority of the data are resolved with more than

10

cells, and the lowest recorded resolutions are in the order of

δ_{t h} / Δ \approx 7

⁠. The flame front is remarkably thinner within the recess segment compared to the combustion chamber. This is because the early mixing and combustion take place at conditions nearer to stoichiometric configurations enabling a higher heat release rate with greater temperature gradients. This is further illustrated in Fig. 4(b), where an instantaneous plot of the flame front resolution is presented. As can be seen, the two symmetrical diffusion flames propagate from their anchoring point at the injection post. The highest resolution is found toward the center, where the mixture is locally stoichiometric and combustion is mostly completed, with the subsequent high viscosity. Therefore, the largest gradients are found at slightly lean and slightly rich mixing conditions. Another relevant aspect to consider is flame stretching. The wrinkles of the flame-front structure can improve or degrade the transport of heat and mass transfer, significantly influencing the local flame dynamics. This can be observed in some zones at the region cavity. The effect of stretching in the flame resolution can be considered to be captured in the probability density functions displayed in Fig. 4(a).

FIG. 4.

View large Download slide

Flame front resolution: Probability density function of the flame front resolution within the recess segment and the combustion chamber for different temperature ranges (a), and example of instantaneous flame front resolution (b).

In transient simulations, the used time step is a crucial parameter to be considered for the simulation accuracy. Generally, there are three potential limiting factors:

Numerical flow stability, given by the convective (momentum), and acoustic (pressure) CFL (Courant–Friedrichs–Levy) numbers.
Heat conduction, given by the Fourier number $F O$ ⁠.
Chemical timescale, given by the highest chemical reaction rates.

In the present simulation, the strictest requirement corresponds to the acoustic CFL number. The simulation time steps are adaptive, ensuring that the maximum acoustic CFL is below 0.95 through all the simulation domain, corresponding to maximum convective CFL numbers on the order of 0.1–0.15. In a previous study⁶⁷ it was proved that this arrangement coupled with a global fourth order discretization scheme ensures negligible numerical dissipation. In the present work, identical differentiation schemes (i.e., cubit 3rd order in space and 1st order in time) were used, allowing to claim a low risk of numerical diffusion.

With the described numerical setup, the solver obtains instantaneous flow fields for every relevant turbulent quantity, such as the ones displayed in Fig. 2 (multimedia available online). These results are averaged in time and space to obtain first and higher-order statistics. Space averaging is done in two different ways:

Due to the splitter-plate configuration, the boundary conditions are constant spanwise, thus precluding the transport of statistics in this direction (y). Therefore, results for all points sharing axial and radial coordinates can be averaged since they have identical stochastics.
Symmetry allows using the data at positions with the same distance to the symmetry to determine any statistic, e.g., $\bar{q} (x) = \bar{q} (- x)$ ⁠. This condition allows doubling the available data.

These spatial averages are applied together with time averaging to determine the turbulent flow statistics. In total, the statistics configure symmetrical 2-dimensional fields. The results for all relevant physical quantities (velocity, density, pressure, temperature, and species concentration) are stored every 0.5

μ s

⁠, during 316

μ s

⁠. Merging time and spatial averages, a total of 333 (time) × 160 (spanwise) × 2 (symmetry) = 202 560 datapoints are used to compute the statistics. Such a large amount of data are generally enough to ensure ergodicity. However, due to the coherent structure of vortical structures over space and time, independency between all the measurements cannot be assumed. The pseudo-number of equivalent measurements at a specific location can be estimated following the procedure discussed in Ref. 12:

n_{ps} = 2 \frac{t_{sim}}{\tilde{k}} \tilde{ε} \frac{L_{y}}{{\tilde{k}}^{3 / 2}} \tilde{ε} = 2 t_{sim} L_{y} \frac{{\tilde{ε}}^{2}}{{\tilde{k}}^{5 / 2}} .

(31)

Using this equivalent number of measurements, it is possible to estimate the error for a fluctuating variable with a given confidence interval

C I

as

ϵ_{CI} = \frac{t_{ps} (CI, n_{ps}) σ}{\sqrt{n_{ps}}} \leftrightarrow P (ϵ < ϵ_{CI}) = CI,

(32)

where

σ

is the Reynolds or Favre-average of the considered variable. Assuming a confidence interval

CI = 95 %

⁠, it is possible to obtain the normalized error value

ϵ_{C I} / σ

⁠. The spatial distribution of this indicator is displayed in Fig. 5(a). As can be seen, the error normalized with the standard deviation is below 0.05 through most of the domain. Since the standard deviation is typically one or two orders of magnitude below the characteristic scale of the considered quantity, it is easy to see that the statistical error will be below 1% in most cases. An example of the estimated error with a confidence interval of 95% is presented in Fig. 5(b). As can be seen, the error is mostly below 0.01, which is sufficiently small considering that

Z

is on the order of unity.

FIG. 5.

$Error in first-order statistics with a confidence interval of 95%: Error normalized with the standard deviation (a) and example of estimated error field for the Favre-average mixture fraction Z ̃ (b).$

View large Download slide

Error in first-order statistics with a confidence interval of 95%: Error normalized with the standard deviation (a) and example of estimated error field for the Favre-average mixture fraction $\tilde{Z}$ (b).

IV. RESULTS

The present section is devoted to the analysis of the simulation results. First, the averaged fields are commented on to discuss the overall turbulent flame development and mixing process. Second, turbulent diffusion is examined in detail to understand the physical motivations behind the observed statistical results. Finally, different modeling approaches for the turbulent diffusion transport term are studied, including new proposed modeling methodologies. The performance of all the different strategies is evaluated, pinpointing their advantages and shortcomings.

A. Turbulent mixing and flame development

The Favre-averaged mixture fraction and its variance are displayed in Fig. 6. During the recess cavity and right-after injection, two symmetrical flames anchored at the prime injection post-tip are formed. These systems display a canonical shear configuration therein, with high concentration gradients localized near the mean stoichiometric line. The reacting mixture flows through the cavity before it is second injected in the combustion chamber. It is here where the flame structure starts to exhibit a remarkable oscillatory behavior. These fluctuations are primarily caused by the interaction between stagnated flow at the injection faceplate and the jet flow leaving the recess cavity.

FIG. 6.

$Statistics of the mixture fraction: average mixture fraction (a) and mixture fraction variance (b).$

View large Download slide

Statistics of the mixture fraction: average mixture fraction (a) and mixture fraction variance (b).

The shear between these regions with different velocities and densities promotes Kelvin–Helmholtz instabilities sort (KHI),^70–73 causing large-scale fluctuations of the mixing layer. The vortex shedding develops with an amplitude in the order of the injector diameter $ϕ_{H}$ ⁠. It is important to mention that this sort of instability takes place within the recess region as well. However, this phenomenon is significantly weaker in this segment due to two principal reasons:

The shear is less pronounced since the recirculation zone at the injection post is significantly smaller compared to that in the injection faceplate.
Instead of a single stream, there are three different jet flows (fuel–oxidizer–fuel) whose individual oscillating behavior produces a mutual dampening.

As a consequence, KHI become the driving mixing mechanism after injection in the combustion chamber. This leads to an increased dominance of turbulent diffusion over its molecular counterpart. This process is manifested through the high variance increase in Fig. 6(b). By $z \approx 3 ϕ_{H}$ ⁠, these fluctuations are high enough to collapse the mean shear layer structure, and turbulent diffusion becomes the dominant mixing mechanism. Turbulent diffusion of fuel to the outer radial positions continues after the average mixture fraction fields reach a homogeneous configuration by $z \approx 4 ϕ_{H}$ ⁠. The continued fuel advection generates lean pockets far from the lean flame core at $z \in [5 ϕ_{H}, 7 ϕ_{H}]$ ⁠. This result evinces negative turbulent diffusivity and the possibility of significant amount of oxygen reaching radial positions far from the flame lean core.

To further investigate the dynamics of the turbulent mixing process, the budget of the average mixture fraction transport is displayed in Fig. 7. In this figure, the total transport is displayed along with the contribution owing to molecular and turbulent diffusion, as formalized in (11). Positive (red) values can be interpreted as a positive net flux of fuel, whereas negative (blue) values indicate that, on average, more oxidizer than fuel is received. As can be seen, molecular diffusivity is relevant near the wall within the recess cavity and in the faceplate region. This is a consequence of the smaller turbulent Reynolds numbers of these regions due to the proximity to the wall, which limits the magnitude of turbulent diffusion.

FIG. 7.

$Transport of the average mixture fraction Z ̃: Total transport (a), molecular diffusivity component (b), and turbulent diffusivity component (c).$

View large Download slide

Transport of the average mixture fraction $\tilde{Z}$ ⁠: Total transport (a), molecular diffusivity component (b), and turbulent diffusivity component (c).

Once the flow enters the combustion chamber, turbulent diffusion becomes the driving mixing mechanism due to the magnitude increase in KHI. Initially, turbulence amplifies mixing, enhancing the transport of fuel from the outer regions to the flame core and vice versa for oxygen. Hence, the advection of fuel and oxidizer is aligned with the average concentration gradients, with a positive flux a mixture fraction from the rich outer regions to the injection lean core. Positive diffusivity ensures the stability of this process, since diffusive advection decreases as the flow homogenizes, and concentration gradients tend to zero. However, turbulent diffusivity can be negative, producing a resegregation of the mixing components. This process can be observed by $z \approx 5 ϕ_{H}$ ⁠, where the fuel and oxidizer separate after reaching an almost homogeneous configuration.

B. Turbulent diffusion dynamics

To understand the observed trends of turbulent diffusion, the spatial variability of this phenomenon and the physical motivations behind it are investigated in this section. For the sake of simplicity, we shall limit the scope of the analysis to radial turbulent diffusion, which is the dominating source. This is claim is validated with the decomposition of the turbulent diffusion transport into radial and axial components represented in Fig. 8. As can be seen, the axial component is negligible compared to the radial one throughout most of the domain. This is an expectable result, due to the null radial velocity component of the injected propellants, which primarily restricts mixing to the radial direction. It shall be noted that in other injection configurations such as impinging, axial turbulent diffusion may play a non-negligible role.

FIG. 8.

View large Download slide

Individual components of the turbulent diffusivity term: radial component (a), axial component (b).

The spatial variability of the turbulent diffusivity may be approached through the analysis of the departure from the turbulent diffusivity model summarized in . The reasoning behind this expression attributes the covariance between mixture fraction and momentum fluctuations to turbulent diffusion originated from mixture fraction gradients. This model has the degree of freedom

C_{μ} \approx 0.09

⁠, which is constant if the model is ideal, with turbulent diffusion behaving analogously to its molecular counterpart. The actual in the frame of DNS, this model parameter can be treated instead as a coefficient that displays different values at each position, since the actual value of the unclosed statistic

\bar{ρ u_{x}^{″} Z ″}

is known,

\frac{C_{μ}}{{Sc}_{T}} = \frac{\tilde{ε}}{\tilde{k}} \frac{\bar{ρ u_{x}^{″} Z ″}}{\frac{\partial \tilde{Z}}{\partial x}} .

(33)

The conventional gradient model is rooted in the flow isotropy assumption since diffusion of turbulent fluctuations from other directions is neglected. For the turbulent transport term

\bar{ρ u_{k}}

⁠, Daly and Harlow⁷⁴ proposed an anisotropic extension of the gradient model. In the current work, a similar adaptation of this anisotropic approach to the present context has been devised. The following formulation for the turbulent flux

\bar{ρ u_{x}^{″} Z ″}

was deduced:

\bar{ρ u_{k}^{″} Z ″} = \frac{C_{s}}{{Sc}_{T}} \frac{\bar{ρ}}{\tilde{ε}} (\frac{\partial \tilde{Z}}{\partial x_{k}} {\tilde{u_{k}^{″} u_{k}^{″}}}^{2} (1 - δ_{ik}) + \frac{\partial \tilde{Z}}{\partial x_{i}} \tilde{u_{k}^{″} u_{k}^{″}} \tilde{u_{i}^{″} u_{i}^{″}} δ_{ik}),

(34)

where

δ_{i k}

is the Kronecker delta, and

C_{s}

is a modeling parameter analogous to

C_{μ}

⁠. Within the present work, the radial component

\bar{ρ u_{x}^{″} Z ″}

reduces to 2D case,

\bar{ρ u_{x}^{″} Z ″} = \frac{C_{s}}{{Sc}_{T}} \frac{\bar{ρ}}{\tilde{ε}} (\frac{\partial \tilde{Z}}{\partial x} {\tilde{u_{x}^{″} u_{x}^{″}}}^{2} + \frac{\partial \tilde{Z}}{\partial z} \tilde{u_{x}^{″} u_{x}^{″}} \tilde{u_{x}^{″} u_{z}^{″}}) .

(35)

The local value of

C_{s} / {Sc}_{T}

can be determined using a similar formulation as . The spatial variability of the introduced model coefficients is displayed in Fig. 9. If gradient models would hold consistently, the coefficients

C_{μ} / S c_{T}

and

C_{s} / S c_{T}

would exhibit positive values throughout all the domain with small or no variations. This is obviously not the case.

FIG. 9.

View large Download slide

Local variability of the turbulent viscosity coefficients $C_{μ} / {S c}_{T}$ (a) and $C_{s} / {S c}_{T}$ (b).

The computed coefficients exhibit a strong variability through all the domain, and negative values can be observed at several locations. The negative values near $z \approx 4 ϕ_{H}$ denote the region where negative turbulent diffusivity is present since it is implied that the average mixture fraction is advected in a sense opposed to the local gradients. There are a few remarkable differences between the behavior of the conventional isotropic model and the one devised in the present work. Near the injection faceplate, the coefficient $C_{μ} / {S c}_{T}$ displays an irregular behavior, including negative values, whereas $C_{s} / {S c}_{T}$ exhibits a more consistent (i.e., regular) distribution. This result indicates the better performance of the anisotropic model, which is an expectable outcome due to the highly anisotropic flow nature in this zone. There are some noticeable performance differences within the recess cavity as well. In this segment, the anisotropic model constant displays a more consistent distribution, which illustrates its capability to capture the mixing layer anisotropy. The improved performance of the anisotropic model compared to the conventional isotropic gradient model can be observed with greater detail in Fig. 15, where the predicted fields of $\bar{ρ u_{x}^{″} Z ″}$ using the best-fitting values for the different models are displayed. For the sake of completeness, the best-fitting values found for the coefficients at the most relevant regions are provided in Table II.

TABLE II.

Best fitting parameters of the model coefficients for the isotropic and anisotropic gradient models.

Region	$C_{μ} / {Sc}_{T}$	$C_{s} / {Sc}_{T}$
Recess region	0.23	0.8
Injection faceplate	0.05	0.17
Combustion camber	0.28	0.2

Region	$C_{μ} / {Sc}_{T}$	$C_{s} / {Sc}_{T}$
Recess region	0.23	0.8
Injection faceplate	0.05	0.17
Combustion camber	0.28	0.2

Once the regions with negative turbulent diffusivity have been clearly identified, we shall analyze their statistical behavior in detail. Its comparison with that of regions where the gradient assumption holds shall shed light on the physical motivations behind this phenomenon. To investigate the invalidity of the gradient model for closure of the turbulent flux $\bar{ρ u_{x}^{″} Z ″}$ the joint probability density function of velocity and mixture gradient fluctuations was computed at specific paradigmatic locations. This result is presented in Fig. 10. Two different points are chosen for the discussion:

A, right after the flow enters the combustion chamber, where the shear layer structure is consistent, and the gradient model holds reasonably well.
C, where the mean shear layer structure has collapsed due to the large magnitude of fluctuations motivated by the KHI, and negative turbulent diffusion is observed.

FIG. 10.

$Joint probability density function of the radial velocity fluctuations and the mixture fraction at two reference points (a). Position of the reference points within the simulation domain (b). Velocity is normalized with the density-weighted average velocity at the combustion chamber injection plane ( ϕ H = 1), which was determined to be u i = 44.015 m/s.$

View large Download slide

Joint probability density function of the radial velocity fluctuations and the mixture fraction at two reference points (a). Position of the reference points within the simulation domain (b). Velocity is normalized with the density-weighted average velocity at the combustion chamber injection plane (⁠ $ϕ_{H} = 1$ ⁠), which was determined to be $u_{i}$ = 44.015 m/s.

At A, the average shear layer is rather consistent and most of the oscillations do not deviate significantly from the mean reference configuration. In such a case, there is an evident positive correlation between radial velocity fluctuations and mixture fraction fluctuations. However, an incipient non-linearity for very high radial velocities can be observed. If the radial velocity displays very high positive values, mixture fraction reductions occur. This is caused by the interaction with the contiguous flame from the opposite direction. If velocity anomalies are sufficiently high, the flame structure can display significant displacements, and the line of maximal mixture fraction can overtake the reference point. In such cases, the mixture is on a configuration driven by the amount of oxygen advected from the contiguous injector element instead of the mixing with the lean side of the nominal flame. This effect is further illustrated with instantaneous examples displayed in Fig. 11(a). As can be seen, if the fluctuations in velocity are moderate, the flame structure oscillates within its rich side. In this reference state, a positive correlation between the fluctuations in velocity and mixture fraction can be observed. For very large anomalies in velocity, the line of maximum mixture fraction is shifted toward the injection core. In this situation, the mixture founds itself at the opposite size of the local maximum, and the radial mixture fraction gradient sign is consequently reversed.

FIG. 11.

View large Download slide

Nominal gradient transport with negligible non-linear counter-gradient features [point A in Fig. 10(b)]: Representation of joint-PDF examples (a) and schematic of the different flow configurations and their relation to the gradient model (b).

The discussed non-linear effect has a reduced magnitude in A, and the overall gradient assumption remains mostly valid. However, at position C, the magnitude of the flame structure oscillations is significantly larger due to the strong onset of KHI. Hence, the mentioned non-linearities are dominant, and the validity of the gradient model is severely challenged. In this position, the mixture fraction oscillates between extreme lean and fuel-rich modes without much statistical relevance left for intermediate mixing configurations. This process is illustrated in Fig. 12. As can be seen, there is an overall positive value of the statistic $\bar{ρ u_{x}^{″} Z ″}$ ⁠, which arises from the global correlation between lean and rich modes. However, the gradient observed when the flame is at one of these specific displays the opposite sign.

FIG. 12.

View large Download slide

Counter-gradient transport due to Kelvin–Helmholtz oscillations [point C in Fig. 10(b)]: Representation of joint-PDF examples (a) and schematic of counter-gradient process (b).

Therefore, the observed negative turbulent diffusion arises from the prevalence of large-scale oscillations driven by the KHI originated at the combustion chamber injection. To further investigate the fluctuating behavior of the mixture fraction, this signal was analyzed in time and frequency domain as represented in Fig. 13. In this figure, the mixture signal time and frequency domain corresponding to the points A, B, C [displayed in Fig. 10(b)] can be observed.

FIG. 13.

$Analysis of mixture fraction fluctuations: Behavior in time domain (a) and power spectral density plot in frequency domain (b). The sub-figures 1–3 correspond to the points A–C, respectively [see Fig. 10(b)].$

View large Download slide

Analysis of mixture fraction fluctuations: Behavior in time domain (a) and power spectral density plot in frequency domain (b). The sub-figures 1–3 correspond to the points A–C, respectively [see Fig. 10(b)].

The spectral analysis reveals three relevant modes, which shall be commented:

KHI following the combustion chamber injection $f_{cc} \approx 32 kHz$ ⁠, controlled by the characteristic timescale of the combustion chamber injection.
KHI following prime injection at the recess cavity $f_{rl} \approx 250 kHz$ ⁠, controlled by the characteristic timescale of the prime injection.
Turbulent mixing dynamics $f_{T} \approx 100 kHz$ ⁠, controlled by the scalar dissipation rate at the flame $χ$ ⁠.

In Fig. 13(b.1), the flame has just entered the combustion chamber, and the modes due to the two different KHI can be spotted. The mode corresponding to the prime injection at the recess cavity is a relic which tends to dissipate since its source is no longer active. By Fig. 13(b.2), the mode $f_{rl}$ remains observable, but its fast vanishing is evident. By this point, the mode due to the turbulent diffusion flame mixing intensity can be observed. By Fig. 13(c), the mode due to the prime injection at the recess cavity is completely dissipated, and the only outstanding modes are those due to the injection in the combustion chamber $f_{cc}$ and turbulent mixing $f_{T}$ ⁠. Similar interactions have been found in previous works.^75,76 In addition, high frequency combustion instabilities in methane propulsion engines have been associated with Kelvin–Helmholtz instabilities in previous works.⁷⁷

The Strouhal number is an important characteristic parameter of flows subject to the KHI. This dimensionless number is usually defined as

Str = \frac{f_{KHI}}{u_{c}} l_{c},

(36)

where

f_{KHI}

denotes the vortex-shedding frequency motivated by the KHI with

u_{c}

and

l_{c}

symbolizing characteristic velocity and length scales. The frequency of the KHI can be easily determined performing a spectral analysis such as the one displayed in Fig. 13(b). However, the definition of characteristic velocity and length scales is not trivial due to the injection of a reacting flow. At the recess cavity, there are three different jet streams with the potential to originate KHI. Since fuel and oxidizer are injected with identical mean bulk velocities, this velocity can be taken as the characteristic one. However, fuel and oxidizer are injected with different injection diameters. An equivalent diameter is defined weighting by the ones of fuel and oxidizer with their momentum

ϕ_{H, eq} = Z_{global} ϕ_{H, f} + (1 - Z_{global}) ϕ_{H, ox},

(37)

where

Z_{global} = 1 / 3

is the global mixture fraction. Using this equivalent hydraulic diameter, a Strouhal number for the recess cavity

Str \approx 0.208

is obtained. For the injection in the combustion chamber, the density-weighed average axial velocity is taken as a reference velocity scale, and the injection hydraulic diameter

ϕ_{H}

is used as reference length scale. The choice of these references provides a Strouhal number

Str \approx 0.216

for the vortex shedding within the combustion chamber. The obtained values are in excellent agreement with the results obtained in vortex shedding studies due to blunt bodies, which consistently reveal Strouhal numbers in the order

Str \approx 0.2

for moderate Reynolds numbers.^78–82

C. Numerical modeling

In the frame of RANS simulations, determining the turbulent diffusion of the average mixture fraction requires from closure models to predict the non-resolved turbulent flux $\bar{ρ u_{x}^{″} Z ″}$ ⁠. In this section, different available modeling methodologies are commented.

The most commonly extended strategy in commercial solvers is to resort to the gradient assumption to predict the unclosed statistic $\bar{ρ u_{x}^{″} Z^{″}}$ ⁠. The obtained fields using such models are displayed in Figs. 14(a), 14(b), 15(a), and 15(b). The best-fitting values for the coefficients $C_{μ} / {Sc}_{T}$ and $C_{s} / {Sc}_{T}$ ⁠, were used to produce the fields $T_{μ}$ and $T_{s}$ ⁠, which denote the predicted value for $\bar{ρ u_{x}^{″} Z ″}$ using the isotropic and anisotropic gradient models, respectively. As can be seen, the models based on the gradient assumption display a poor performance through the combustion chamber region, and they fail to predict the negative turbulent diffusivity toward $z \approx 5 ϕ_{H}$ ⁠. Nonetheless, the gradient models, especially the anisotropic one, display a very good performance within the recess region, as observed in Figs. 15(a) and 15(b). The reported performance differences are tied to the dominance of large-scale fluctuations associated with the vortex-shedding frequency.

FIG. 14.

Performance of different numerical models for the turbulent scalar flux ρ u x ″ Z ″ ¯ over the entire simulation domain: Isotropic gradient model (a), anisotropic gradient model (b), mixing-based filtered RANS model (c), closed mixing-based filtered RANS model (d), and observed turbulent scalar flux (e).

View large Download slide

Performance of different numerical models for the turbulent scalar flux $\bar{ρ u_{x} ″ Z ″}$ over the entire simulation domain: Isotropic gradient model (a), anisotropic gradient model (b), mixing-based filtered RANS model (c), closed mixing-based filtered RANS model (d), and observed turbulent scalar flux (e).

FIG. 15.

Performance of different numerical models for the turbulent scalar flux ρ u x ″ Z ″ ¯ within the recess segment: Isotropic gradient model (a), anisotropic gradient model (b), mixing-based filtered RANS model (c), closed mixing-based filtered RANS model (d), and observed turbulent scalar flux (e).

View large Download slide

Performance of different numerical models for the turbulent scalar flux $\bar{ρ u_{x} ″ Z ″}$ within the recess segment: Isotropic gradient model (a), anisotropic gradient model (b), mixing-based filtered RANS model (c), closed mixing-based filtered RANS model (d), and observed turbulent scalar flux (e).

The difficulties of standard models to capture the statistical behavior of the KHI in shear coaxial injectors for space propulsion applications poses a significant challenge. The relevance of this phenomenon for initial mixing and combustion has been noted by several previous researchers.^83–87 In addition, high-frequency combustion instabilities have been associated with KHI after injection.⁷⁷ Moreover, the blanching issue mentioned in the introduction is likely related to this phenomenon, since the advection of oxygen to the combustion chamber walls must involve some sort of negative turbulent diffusion. Hence, exploring modeling improvements which can capture this process is of vital importance.

After the recess cavity, the statistic

\bar{ρ u_{x}^{″} Z ″}

is primarily driven by large-scale fluctuations between lean and rich modes. Therefore, it is more suitable to rewrite this term using the mixing-based filtering strategy presented in the introduction of this paper

\bar{ρ u_{x}^{″} Z ″} = \bar{ρ} (\tilde{λ} (\tilde{Z_{λ}^{″} u_{x, λ}^{″}} + Δ {\tilde{u}}_{x, λ} Δ {\tilde{Z}}_{λ}) + \tilde{φ} (\tilde{Z_{φ}^{″} u_{x, φ}^{″}} + Δ {\tilde{u}}_{x, φ} Δ {\tilde{Z}}_{φ})) .

(38)

Since fluctuations between lean and rich modes are dominant, this expression can be simplified as

\bar{ρ u_{x}^{″} Z ″} \approx \bar{ρ} (\tilde{λ} Δ {\tilde{u}}_{x, λ} Δ {\tilde{Z}}_{λ} + \tilde{φ} Δ {\tilde{u}}_{x, φ} Δ {\tilde{Z}}_{φ}) .

(39)

The prediction obtained with this expression is denoted as

T_{θ}

⁠, and it is displayed in Figs. 14(c) and 15(c). As can be seen, this model exhibits excellent prediction capabilities outside the recess segment, where large-scale fluctuations are the dominant mixing mechanism. This model only struggles to deliver accurate predictions within the recess cavity and right after injection, before the onset of the large-scale vortex shedding. In these locations, microscale turbulent mixing is non-negligible, and the filtered covariance

\tilde{Z_{θ}^{″} u_{x, θ}^{″}}

owing to microscale turbulent mixing is relatively significant. The application of the proposed model requires to split the transport equations for lean and fuel-rich conditions, to obtain the average fields of turbulent quantities conditional to the mixing state. The most immediate consequence of such an approach is the doubling of the computational cost. This is probably an affordable price given the evident performance improvements. However, the implementation of such a strategy may require a significant programing effort for the available solvers. Moreover, it might be necessary to investigate the performance of extended turbulence models when applied to lean and fuel-rich time averages instead of global averages. For these reasons, it is convenient to approach the closure of

\tilde{u_{x} ″ Z^{″}}

using global averages while capturing the dominance of lean-rich fluctuations manifested in . The main idea is to model the right-hand side of this expression solely resorting to closed global averages. In the frame of conventional RANS solvers, the statistical quantities

\tilde{λ}, \tilde{φ}, {\tilde{Z}}_{λ}, {\tilde{Z}}_{φ}, {\tilde{u}}_{x, λ}

and

{\tilde{u}}_{x, φ}

are not accessible. However, the first four quantities (i.e.,

\tilde{λ}, \tilde{φ}, {\tilde{Z}}_{λ}, {\tilde{Z}}_{φ}

⁠) can be modeled quite accurately following the procedure described in the annex of the present text. Moreover, the mean velocities at lean and fuel-rich conditions are linked through

\tilde{u} = \tilde{λ} {\tilde{u}}_{x, λ} + \tilde{φ} {\tilde{u}}_{x, φ}

⁠. Hence, only one additional equation is required to close the quantity

\tilde{Z^{″} u_{x}^{″}}

⁠. If the fluctuations in velocity and mixture fraction are perfectly correlated, this covariance becomes the product of the standard deviations of the involved scalars

σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}} = \sqrt{\tilde{u_{x}^{″} u_{x}^{″}} \tilde{Z^{″ 2}}},

(40)

which is a closed quantity. The overall correlation coefficient between fluctuations of radial momentum and mixture fraction may be written as

r = \frac{\tilde{u_{x}^{″} Z ″}}{σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}}} \approx \frac{\tilde{λ} Δ {\tilde{u}}_{x, λ} Δ {\tilde{Z}}_{λ} + \tilde{φ} Δ {\tilde{u}}_{x, φ} Δ {\tilde{Z}}_{φ}}{σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}}},

(41)

which is evidently unclosed. Since

r

denotes the correlation between radial momentum fluctuations and mixture fraction fluctuations, it is reasonable to assume that it scales with the average flux in fuel-rich conditions. Hence, the following closure model can be formulated:

r \approx κ \frac{{\tilde{u}}_{x, φ} \tilde{φ}}{{\tilde{u}}_{inj}},

(42)

where

κ \approx 4.48

is a modeling parameter. Applying fundamental relationships, it is possible to completely close the unresolved covariance

\tilde{u_{x}^{″} Z ″}

with the following formula:

\tilde{u_{x}^{″} Z ″} = σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}} r = σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}} κ \frac{{\tilde{u}}_{x} ({\tilde{Z}}_{λ} - \tilde{Z})}{{\tilde{Z}}_{λ} - {\tilde{Z}}_{φ} + κ σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}}} .

(43)

For numerical stability reasons, a floor value of unity is set to the denominator' absolute value, i.e.,

| {\tilde{Z}}_{λ} - {\tilde{Z}}_{φ} + κ σ_{{\tilde{u}}_{x}} σ_{\tilde{Z}} | \geq 1

⁠. This model can be further refined by considering the role of the mixture fraction degree of bimodality. For a given average mixture fraction

\tilde{Z},

the maximum possible variance corresponds to the case in which the mixture fraction behaves as a bimodal fluctuating variable, which can be determined as

{\tilde{Z^{″ 2}}}_{max} = \tilde{Z} (1 - \tilde{Z}) .

(44)

The mixture fraction variance is hence ceiled by this value. Therefore, the proximity to a bimodal distribution can be evaluated as the ratio between the observed variance and the one associated with the bimodal case

b = \frac{\tilde{Z^{″ 2}}}{\tilde{Z} (1 - \tilde{Z})},

(45)

which displays a value

b \approx 1

if the mixture fraction oscillates between

\tilde{Z} \approx 1

and

\tilde{Z} = 0

and intermediate values are stochastically negligible. This parameter can be assumed to scale with the correlation coefficient

r

⁠, since intermediate values contribute to the degradation of the overall correlation. Hence, a more accurate model may be expressed as

\tilde{u_{x}^{″} Z ″} = f {\tilde{u}}_{i} \frac{{\tilde{Z}}_{λ} - \tilde{Z}}{\tilde{φ} (\frac{κ}{u_{inj}} - \frac{{\tilde{Z}}_{φ} - {\tilde{Z}}_{λ}}{σ})},

(46)

where

f \approx 2.41 b

is a correction factor accounting for the correlation degradation given by the presence of intermediate values. The results of this model can be compared against the observed statistic in Figs. 14(d) and 15(d), where the predicted turbulent flux

\bar{ρ u_{x}^{″} Z ″}

⁠. As can be seen, the globally closed model is able to predict the counter-gradient transport of the turbulent flux, displaying a very good performance in the near injection region. However, this model has poor prediction capabilities within the recess region. Moreover, the proposed model displays a significantly worsened performance after the mixing layers collapse near

z \approx 5 ϕ_{H}

⁠. This model is presented to illustrate possible global closure strategies which address the fluctuating behavior between lean and rich modes formalized in manifested in .

The discussed modeling strategies have been focused on the first-order turbulent flux $\bar{ρ u_{i}^{″} Z ″}$ ⁠. Nonetheless, they can be applied as well for the term $\bar{ρ u_{i}^{″} Z^{″ 2}}$ ⁠, which is relevant for the transport of the mixture fraction variance. The different models display identical strong points and shortcomings when used to predict the field $\bar{ρ u_{i}^{″} Z^{″ 2}}$ and the results are not reproduced within this document for the sake of brevity.

V. CONCLUSIONS

Turbulent mixing in a high-pressure methane–oxygen recessed injector was investigated using Direct Numerical Simulations (DNS). The operating conditions were chosen to resemble those characteristics of modern space propulsion engines within the existing computational power limitations. The resolution requirements and the statistical error of the obtained database were thoroughly assessed to certify the fidelity of the results. It was observed that turbulent diffusion is the dominant mixing mechanism through most of the domain, apart from near-wall regions, where molecular diffusivity plays a relevant role. Moreover, mixing is strongly anisotropic, with the axial component being negligible compared to the radial one. Negative turbulent diffusivity was observed in specific regions, where the mixture rapidly oscillates from lean to fuel-rich configurations due to the vortex shedding process driven by Kelvin–Helmholtz instabilities. The dominance of these large-scale fluctuations yields the observed negative turbulent diffusivity. The subsequent counter-gradient diffusion produces an excessive advection of the oxidizer toward the originally fuel-rich regions, creating large lean pockets. This process could be behind the observed blanching in fuel-rich rocket engines. However, the conducted simulations did not consider the combustion wall region, which is where this process occurs. Hence, further investigations are required to confirm the implication of vortex-shedding mechanisms in the blanching onset. Conventional gradient models struggle to capture the transport of the average mixture fraction in the regions where vortex shedding-mixing is the driving turbulent mixing source. An extension of the gradient model to anisotropic configurations was proposed to capture the anisotropic turbulent mixing dynamics. The devised model provides significant improvements within the recess section of the injector and near the faceplate, where a fuel-rich mixture recirculates. Nonetheless, it fails in capturing counter-gradient diffusion similarly to the original isotropic model version. A novel RANS formulation using mixing-based filters was proposed as well. This new approach is based on the decomposition of the conservation equations in lean and fuel-rich conditions. The suggested model is able to accurately predict turbulent diffusion after injection in the combustion chamber, and it presents a decent performance within the recess segment. However, this approach entails significant implementation challenges since the available RANS solvers are conceived for global averages. One possible strategy to link the developed filtered RANS model to global averages is presented for the sake of completeness. Such an approach enables the closure of global quantities while using a set of equations suited to capture the dominance of large-scale fluctuations. Nonetheless, further investigations are required in this direction.

ACKNOWLEDGMENTS

The authors thank Super-MUC NGen for providing the computational resources for performing the numerical simulations and its post-processing. In particular, the authors are extremely grateful toward Martin Ohlerich for his continued support, which allowed this team to make the most of this HPC framework.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Daniel Martinez-Sanchis: Conceptualization (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead). Andrej Sternin: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Agnes Jocher: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Oskar Haidn: Methodology (supporting); Supervision (equal).

DATA AVAILABILITY

The authors will provide data upon reasonable request

APPENDIX: MODELING OF MEAN FILTERED VALUES

1. Demostration of relationship between Reynolds and Favre averages of mixing filters

The relationship between Reynolds and average filters embedded in Eq. is a fundamental part of the proposed approach. We first take advantage of the property

\bar{λ} + \bar{φ} = 1

⁠, to rewrite

\bar{λ}

as

\bar{λ} = \frac{1}{1 + \frac{\bar{φ}}{\bar{λ}}} .

(A1)

Hence, we shall determine the ratio

\bar{φ} / \bar{λ}

dependent on closed variables. Using the definitions for

{\bar{ρ}}_{θ}

\{\begin{matrix} {\bar{ρ}}_{λ} = \bar{ρ λ} / \bar{λ} \\ {\bar{ρ}}_{φ} = \bar{ρ φ} / \bar{φ} \end{matrix} \to \frac{\bar{φ}}{\bar{λ}} = \frac{\bar{ρ φ} {\bar{ρ}}_{λ}}{{\bar{ρ}}_{φ} \bar{ρ λ}} = \frac{\bar{ρ φ} {\bar{ρ}}_{λ} \bar{ρ}}{{\bar{ρ}}_{φ} \bar{ρ λ} \bar{ρ}} = \frac{{\bar{ρ}}_{λ} \tilde{φ}}{{\bar{ρ}}_{φ} \tilde{λ}} .

(A2)

This relationship can be inserted in (A1) to obtain the final relationship,

\bar{λ} = \frac{1}{1 + \frac{{\bar{ρ}}_{λ} \tilde{φ}}{{\bar{ρ}}_{φ} \tilde{λ}}} = \frac{\tilde{λ} {\bar{ρ}}_{φ}}{\tilde{λ} {\bar{ρ}}_{φ} + {\bar{ρ}}_{λ} \tilde{φ}} .

(A3)

An analogous procedure can be conducted to obtain

\bar{φ}

⁠.

2. Modeling of mean filtered values

To completely close the set of filtered RANS equations, it is necessary to obtain the Favre-average of the lean (⁠

λ

⁠) and fuel-rich

(φ)

filters. In the case of the flamelet approach, it is typically assumed that the mixture fraction PDF follows a beta distribution, which can be expressed as

{\tilde{f}}_{β} (Z) = \frac{Γ (a + b)}{Γ (a) Γ (b)} Z^{a - 1} {(1 - Z)}^{b - 1} .

(A4)

With

Γ (c) = \int_{0}^{+ \infty} t^{c - 1} e^{- t} d t

being the gamma function. This function is an extension of

Γ (n) = (n - 1)! \forall n \in N

defined in the complex domain apart from non-positive integers.⁸⁸ Moreover, the parameters

a

and

b

in (B1) are shape factors, which provide flexibility to the function. The corresponding cumulative distribution function (CDF) can be written as

{\tilde{F}}_{β} (Z) = \frac{B (x; a, b)}{B (a, b)} = I_{x} (a, b),

(A5)

where

B (x; a, b)

is the incomplete beta function,

B (x; a, b) = \int_{0}^{x} t^{a - 1} {(1 - t)}^{b - 1} d t .

(A6)

The shape parameters

a

and

b

can be easily determined using the average mixture fraction and its variance,

a = \tilde{Z} (\frac{\tilde{Z} (1 - \tilde{Z})}{\tilde{Z^{″ 2}}} - 1),

(A7)

and

b = \tilde{Z} (\frac{\tilde{Z} (1 - \tilde{Z})}{\tilde{Z^{″ 2}}} - 1) (\frac{1}{\tilde{Z}} - 1) .

(A8)

Hence, using the closed statistics

\tilde{Z}

⁠, and

\tilde{Z^{″ 2}}

⁠, it is possible to estimate the Favre-average of a filter

θ

as

\tilde{θ} \approx \int_{0}^{1} {\tilde{f}}_{β} (Z) θ (Z) d Z .

(A9)

This methodology was applied to the current dataset. The predicted and observed fields for the Favre-average lean filter

\tilde{λ}

are displayed in Fig. 16. As can be seen, the beta distribution assumption is able to capture the integral quantity

\tilde{λ}

with minimal error. The maximum found offset between predicted and observed average filter values is in the order of

| {\tilde{λ}}_{pred} - \tilde{λ} | \approx 0.01

⁠.

FIG. 16.

Predicted (a) and observed (b) Favre-averaged field of the lean filter λ ̃ through the domain of the performed simulations.

View large Download slide

Predicted (a) and observed (b) Favre-averaged field of the lean filter $\tilde{λ}$ through the domain of the performed simulations.

Other average quantities can be determined using the beta-distribution assumption. Generally, the Favre-average of any turbulent scalar conditional to the filter

θ

can be estimated as

{\tilde{q}}_{θ} = \frac{\int_{0}^{1} θ (Z) \bar{q (Z)} d Z}{\tilde{θ}},

(A10)

where

\bar{q (Z)}

denotes the average value of the fluctuating quantity

q

conditional to the mixture fraction. As long as this function is known, the expression in (B7) provides excellent prediction capabilities, such as the one displayed in Fig. 16. For example, the average mixture fraction through a filter can be easily obtained since it is obvious that

\bar{Z (Z)} = Z

⁠. Moreover, in the case of infinitely fast chemistry several quantities can be retrieved from tabulated functions such as

\bar{ρ (Z)}

or

\bar{T (Z)}

⁠.

REFERENCES

1.

D.

Haeseler

,

V.

Bombelli

,

P.

Vuillermoz

,

R.

Lo

,

T.

Marée

, and

F.

Caramelli

, “

Green propellant propulsion concepts for space transportation and technology development needs

,” in

2nd International Conference on Green Propellants for Space Propulsion

(

Chia Laguna

,

Cagliari

,

2004

).

Google Scholar

2.

H. N. M.

Ciezki

and

L.

Werling

, “

Trends in research and development on green chemical propulsion for orbital systems

,” in

7th International Conference on Recent Advances in Space Technologies, RAST 2015

,

Istanbul, Turkey

,

2015

.

Google Scholar

3.

C.

Scharlemann

, “

GRASP: Status and future of green propellants

,” in

Space Propulsion Conference

,

Bordeaux, France

,

2012

.

Google Scholar

4.

M. J. L.

Turner

, “

Liquid propellant rocket engines

,” in

Rocket and Spacecraft Propulsion

(

Springer

,

2009

), pp.

67

–

108

.

Google Scholar

5.

R. L.

Sackheim

and

R. K.

Masse

, “

Green propulsion advancement—Challenging the maturity of monopropellant hydrazine

,” in

49th AIAA Joint Propulsion Conference

(

AIAA

,

San Jose, CA, USA

,

2013

).

Google Scholar

Crossref

6.

O. J.

Haidn

,

M. P.

Celano

,

M.

Luo

,

C.

Roth

,

S.

Silvestri

, and

N. A.

Slavinskaya

, “

On methane/oxygen combustion for rocket applications

,” in

International Symposium on Innovation and Prospects of Liquid Propulsion

(

IEEE

,

Xi'an

,

2016

).

Google Scholar

7.

H.

Burkhardt

,

M.

Sippel

,

A.

Herbertz

, and

J.

Klevanski

, “

Kerosene vs methane: A propellant tradeoff for reusable liquid booster stages

,”

J. Spacecr. Rockets

41

(

5

),

762

–

769

(

2004

).

https://doi.org/10.2514/1.2672

Google Scholar

Crossref

8.

N.

Perakis

,

D.

Rahn

,

O. J.

Haidn

, and

D.

Eiringhaus

, “

Heat transfer and combustion simulation of seven-element O₂/CH₄ rocket combustor

,”

J. Propul. Power

35

(

1

),

1

–

18

(

2019

).

https://doi.org/10.2514/1.B37402

Google Scholar

Crossref

9.

A.

Sternin

,

N.

Perakis

,

M. P.

Celano

, and

O.

Haidn

, “

CFD-analysis of the effect of a cooling film on flow and heat transfer characteristics in a GCH4/GOX rocket combustion chamber

,” in

Space Propulsion 2018

,

Sevilla

,

2018

.

Google Scholar

10.

D.

Maestro

,

B.

Cuenot

,

A.

Chemnitz

,

T.

Sattelmayer

,

C.

Roth

,

O. J.

Haidn

,

Y.

Daimon

,

G.

Frank

,

H.

Müller

,

J.

Zips

,

M.

Pfitzner

,

R.

Keller

,

P.

Gerlinger

,

H.

Riedmann

, and

L.

Selle

, “

Numerical investigation of flow and combustion is a single-element GCH4/GOX rocket combustor: Chemistry modeling and turbulence-combustion interaction

,” AIAA Paper No. 2016-4995,

2016

.

Google Scholar

11.

D.

Martínez-Sanchis

,

S.

Banik

,

A.

Sternin

,

D.

Sternin

,

O.

Haidn

, and

M.

Tajmar

, “

Analysis of turbulence generation and dissipation in shear layers of methane-oxygen diffusion flames using direct numerical simulations

,”

Phys. Fluids

34

,

045121

(

2022

).

https://doi.org/10.1063/5.0087887

Google Scholar

Crossref

12.

D.

Martínez-Sanchis

,

A.

Sternin

,

O. J.

Haidn

, and

M.

Tajmar

, “

Turbulent combustion statistics in a diffusion flame for space propulsion applications

,”

Phys. Fluids

34

(

12

),

125115

(

2022

).

https://doi.org/10.1063/5.0130537

Google Scholar

Crossref

13.

R. W.

Bilger

, “

The structure of diffusion flames

,”

Combust. Sci. Technol.

13

,

155

–

170

(

1976

).

https://doi.org/10.1080/00102207608946733

Google Scholar

Crossref

14.

N.

Peters

, “

Laminar diffusion flamelet models in non-premixed turbulent combustion

,”

Prog. Energy Combust. Sci.

10

(

3

),

319

–

339

(

1984

).

https://doi.org/10.1016/0360-1285(84)90114-X

Google Scholar

Crossref

15.

B.

Cuenot

, “

The flamelet model for non-premixed combustion

,” in

Turbulent Combustion Modelling

(

Springer Nature

,

2011

), pp.

43

–

61

.

Google Scholar

Crossref

16.

N.

Slavinskaya

,

M.

Abbasi

,

J. H.

Starcke

, and

O.

Haidn

, “

Methane skeletal mechanism for space propulsion applications

,” AIAA Paper No. 2016-4781,

2016

.

Google Scholar

17.

D.

Martinez-Sanchis

, “

A flame control method for direct numerical simulations of reacting flows in rocket engines

,”

M.S. thesis

,

Institute of Space Propulsion

,

2021

.

Google Scholar

18.

A.

Bee

and

M.

Börner

, “

Laminar burning speeds and flammability limits of CH₄/O₂ mixtures with varying N₂ dilution at sub-atmospheric conditions

,”

Combust. Sci. Technol.

195

,

1910

–

1929

(

2023

).

https://doi.org/10.1080/00102202.2021.2006191

Google Scholar

Crossref

19.

D. B.

Morgan

and

A. C.

Kobayashi

, “

Main chamber combustion and cooling technology study

,” NASA Report No. NASA-CR-184345,

1989

.

Google Scholar

20.

M.

Murphy

,

R.

Anderson

,

D. C.

Rousar

, and

J. A.

Van Kleeck

, “

Effects of oxygen/hydrogen combustion chamber environment on copper alloys

,” in

Advanced Earth-to-Orbit Propulsion Technology

(

NASA

,

Huntsville

,

1986

).

Google Scholar

21.

L.

Ogbuji

, “

A table-top technique for assessing the blanching resistance of Cu alloys

,”

Oxid. Met.

63

(

5

),

383

–

399

(

2005

).

https://doi.org/10.1007/s11085-005-4393-2

Google Scholar

Crossref

22.

R.

Cook

,

E.

Fryk

, and

J.

Newell

, “

SSME main combustion chamber life prediction

,” Report No. NASA-CR-168215,

1983

.

Google Scholar

23.

D.

Martinez-Sanchis

, “

Investigation and mathematical formalization of roughening processes in combustion chambers of liquid rocket engines

,” Ph.D. thesis, Technical University of Munich, Garching,

2019

.

Google Scholar

24.

K.

Bubenheim

,

C.

Wilhelmi

,

S.

Beyer

, and

S.

Schmidt-Wimmer

, “

ERBURIG test facility: The next step of material testing for H₂/O₂ rocket combustion chambers

,” in

5th European Conference for Aeronautics and Space Sciences (Eucass)

(

EUCASS

,

Munich

,

2013

).

Google Scholar

25.

M. F.

Pouliquen

, “

HM60 cryogenic rocket engine for future European launchers

,”

J. Spacecr. Rockets

21

(

4

),

346

–

353

(

1984

).

https://doi.org/10.2514/3.25661

Google Scholar

Crossref

26.

O. J.

Haidn

, “

Advanced rocket engines

,” Report No. RTO-EN-AVT-150-06, Lampoldshausen, Germany,

2008

.

Google Scholar

27.

F.

Hötte

,

C.

Von Sethe

,

T.

Fiedler

,

M. C.

Haupt

,

O. J.

Haidn

, and

M.

Rohdenburg

, “

Experimental lifetime study of regeneratively cooled rocket chamber walls

,”

Int. J. Fatigue

138

,

105649

(

2020

).

https://doi.org/10.1016/j.ijfatigue.2020.105649

Google Scholar

Crossref

28.

D.

Denker

,

A.

Attili

,

M. N. K.

Gauding

, and

M.

Bode

, “

A new modeling approach for mixture fraction statistics based on dissipation elements

,”

Proc. Combust. Inst.

38

(

2

),

2681

–

2689

(

2021

).

https://doi.org/10.1016/j.proci.2020.06.208

Google Scholar

Crossref

29.

H.

Wang

,

P.

Zhang

, and

J.

Tao

, “

Examination of probability distribution of mixture fraction in LES/FDF modelling of a turbulent partially premixed jet flame

,”

Combust. Theory Modell.

26

,

320

–

337

(

2022

).

https://doi.org/10.1080/13647830.2021.2016980

Google Scholar

Crossref

30.

D.

Martinez-Sanchis

,

A.

Sternin

,

D.

Balbuena-Silvestre

, and

O. J.

Haidn

, “

Mixture fraction statistics in methane-oxygen turbulent combustion for space propulsion

,”

Aerosp. Sci. Technol.

139

(

3

),

108355

(

2023

).

https://doi.org/10.1016/j.ast.2023.108355

Google Scholar

Crossref

31.

T.

Poinsot

and

D.

Veynante

, “

Turbulent non premixed flames

,” in

Theoretical and Numerical Combustion

(

Aquaprint

,

Bordeaux, France

,

2012

), pp.

287

–

348

.

Google Scholar

32.

W. P.

Jones

and

J. H.

Whitelaw

, “

Calculation methods for reacting turbulent flows: A review

,”

Combust. Flame

48

,

1

–

26

(

1982

).

https://doi.org/10.1016/0010-2180(82)90112-2

Google Scholar

Crossref

33.

Y.

Hardalupas

,

M.

Tagawa

, and

M. K. P.

Taylor

, “

Characteristics of countergradient heat transfer in nonpremixed swirling flame

,” in

Developments in Laser Techniques and Applications to Fluid Mechanics

(

Springer

,

Berlin

,

1996

), pp.

159

–

184

.

Google Scholar

34.

K. H.

Luo

, “

On local countergradient diffusion in turbulent diffusion flames

,” in

28th Symposium (International) on Combustion

(

Combustion Institute

,

2000

), pp.

489

–

498

.

Google Scholar

Crossref

35.

S. B.

Pope

, “

The statistical description of turbulent flows

,” in

Turbulent Flows

(Cambridge University Press,

2000

), pp.

34

–

82

.

Google Scholar

Crossref

36.

M.

Eckert

, “

Chaos and turbulence

,” in

The Turbulence Problem

(Springer Cham,

2019

), pp.

75

–

85

.

Google Scholar

Crossref

37.

A.

Favre

,

Problems of Hydrodynamics and Continuum Mechanics

(

SIAM

,

Philadelphia

,

1969

).

Google Scholar

38.

S. B.

Pope

, “

The scales of turbulent motion

,” in

Turbulent Flows, Cambridge

(

Cambridge University Press

,

2000

), pp.

182

–

263

.

Google Scholar

Crossref

39.

P. K.

Yeung

and

S. B.

Pope

, “

Lagrangian statistics from direct numerical simulations of isotropic turbulence

,”

J. Fluid Mech.

207

,

531

–

586

(

1989

).

https://doi.org/10.1017/S0022112089002697

Google Scholar

Crossref

40.

R. W.

Bilger

, “

Turbulent jet diffusion flames

,”

Prog. Energy Combust. Sci.

1

(

2–3

),

87

–

109

(

1976

).

https://doi.org/10.1016/0360-1285(76)90022-8

Google Scholar

Crossref

41.

S. B.

Pope

, “

Turbulent-viscosity models

,” in

Turbulent Flows

(

Cambridge University Press

,

2000

), pp.

358

–

386

.

Google Scholar

Crossref

42.

W. P.

Jones

and

B. E.

Launder

, “

The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence

,”

Int. J. Heat Mass Transfer

16

(

6

),

1119

–

1130

(

1973

).

https://doi.org/10.1016/0017-9310(73)90125-7

Google Scholar

Crossref

43.

D. C.

Wilcox

,

Turbulence Modelling for CFD

(

DCW

,

La Cañada, CA

,

2002

).

Google Scholar

44.

J.

O'Brien

,

C. A. Z.

Towery

,

P. E.

Hamlington

,

M.

Ihme

,

A.

Poludnenko

, and

J.

Urzay

, “

The cross-scale physical-space transfer of kinetic energy in turbulent premixed flames

,”

Proc. Combust. Inst.

36

(

2

),

1967

(

2017

).

https://doi.org/10.1016/j.proci.2016.05.005

Google Scholar

Crossref

45.

C. A. Z.

Towery

,

A. Y.

Poludnenko

,

J.

Urzay

,

J.

O'Brien

,

M.

Ihme

, and

P. E.

Hamlington

, “

Spectral kinetic energy transfer in turbulent premixed reacting flows

,”

Phys. Rev. E

93

,

053115

(

2016

).

https://doi.org/10.1103/PhysRevE.93.053115

Google Scholar

Crossref

PubMed

46.

J.

O'Brien

,

J.

Urzay

,

M.

Ihme

,

P.

Moin

, and

A.

Saghafian

, “

Subgrid-scale backscatter in reacting and inert supersonic hydrogen-air turbulent mixing layers

,”

J. Fluid Mech.

743

,

554

–

584

(

2014

).

https://doi.org/10.1017/jfm.2014.62

Google Scholar

Crossref

47.

D.

Martinez-Sanchis

,

A.

Sternin

,

D.

Sternin

,

O.

Haidn

, and

M.

Tajmar

, “

Analysis of periodic synthetic turbulence generation and development for direct numerical simulations applications

,”

Phys. Fluids

33

,

125130

(

2021

).

https://doi.org/10.1063/5.0071002

Google Scholar

Crossref

48.

C.

Morsbach

and

M.

Franke

, “

Analysis of a synthetic turbulence generation method for periodic configurations

,” in

Direct and Large-Eddy Simulation XI

(

Springer International Publishing

,

2019

), pp.

169

–

174

.

Google Scholar

Crossref

49.

M. L.

Shur

,

P. R.

Spalart

,

M. K.

Strelets

, and

A. K.

Travin

, “

Synthetic turbulence generators for RANS-LES interfaces in zonal simulations of aerodynamic and aeroacoustic problems

,”

Flow, Turbul. Combust.

93

,

63

–

92

(

2014

).

https://doi.org/10.1007/s10494-014-9534-8

Google Scholar

Crossref

50.

R. H.

Kraichnan

, “

Diffusion by a random velocity field

,”

Phys. Fluids

13

,

22

–

31

(

1970

).

https://doi.org/10.1063/1.1692799

Google Scholar

Crossref

51.

S.

Silvestri

,

M. P.

Celano

,

G.

Schlieben

,

C.

Kirchberger

, and

O.

Haidn

, “

Characterization of a GOX-GCH4 single element combustion chamber

,” in

Space Propulsion

,

Cologne, Germany

,

2014

.

Google Scholar

52.

S.

Silvestri

,

F. G. M.

Winter

,

C. M.

Palma

,

G.

Schlieben

, and

O.

Haidn

, “

Investigation on recess variation of a shear coax injector in a GOX/GCH₄ rectangular combustion chamber

,” in

7th European Conference for Aeronautics and Space Science (EUCASS)

,

Milano

,

2017

.

Google Scholar

53.

N.

Chakraborty

,

M.

Katragadda

, and

R. S.

Cant

, “

Statistics and modelling of turbulent kinetic energy transport in different regimes of premixed combustion

,”

Flow, Turbul. Combust.

87

,

205

–

235

(

2011

).

https://doi.org/10.1007/s10494-010-9312-1

Google Scholar

Crossref

54.

N.

Chakraborty

, “

Influence of thermal expansion on fluid dynamics of turbulent premixed combustion and its modelling implications

,”

Flow, Turbul. Combust.

106

,

753

–

848

(

2021

).

https://doi.org/10.1007/s10494-020-00237-8

Google Scholar

Crossref

55.

F.

Zhang

,

H.

Bonart

,

T.

Zirwes

,

P.

Habisreuther

, and

H.

Bockhorn

, “

Direct numerical simulations of chemically reacting flows with the public domain code OpenFOAM

,” in

High Performance Computing in Science and Engineering '16

(

Springer

,

2015

), Vol.

14

, pp.

221

–

236

.

Google Scholar

Crossref

56.

T.

Zirwes

,

F.

Zhang

,

P.

Habisreuther

,

J. A.

Denev

,

H.

Bockhorn

, and

D.

Trimis

, “

Implementation and validation of a computationally efficient DNS solver for reacting flows in OpenFOAM

,” in

14th World Congress on Computational Mechanics

,

Virtual Congress

,

2021

.

Google Scholar

Crossref

57.

H.

Weller

,

G.

Tabor

,

H.

Jasak

, and

C.

Fureby

, “

A tensorial approach to computational continuum mechanics using object-oriented techniques

,”

Comput. Phys.

12

(

6

),

620

–

631

(

1998

).

https://doi.org/10.1063/1.168744

Google Scholar

Crossref

58.

H.

Weller

,

G.

Tabor

,

H.

Jasak

, and

C.

Fureby

,

OpenFOAM

(

openCFD ltd

,

2017

).

Google Scholar

59.

P. W.

Mcdonald

, “

The computation of transonic flow through two-dimensional gas turbine cascades

,” in

ASME 1971 International Gas Turbine Conference and Products Show

(

ASME

,

1971

).

Google Scholar

Crossref

60.

R. W.

MacCormack

and

A. J.

Paullay

, “

Computational efficiency achieved by time splitting of finite difference operators

,” AIAA Paper No. 72-154,

1972

.

Google Scholar

61.

F.

Zhang

,

T.

Zirwes

,

P.

Habisreuther

, and

H.

Bockhorn

, “

Effect of unsteady stretching on the flame local dynamics

,”

Combust. Flame

175

,

170

–

179

(

2017

).

https://doi.org/10.1016/j.combustflame.2016.05.028

Google Scholar

Crossref

62.

T.

Zirwes

,

T.

Häber

,

Z.

Feichi

,

H.

Kosaka

,

A.

Dreizler

,

M.

Steinhausen

,

C.

Hasse

,

A.

Stagni

,

D.

Trimis

,

R.

Suntz

, and

H.

Bockhorn

, “

Numerical study of quenching distances for side-wall quenching using detailed diffusion and chemistry

,”

Flow, Turbul. Combust.

106

,

649

–

679

(

2021

).

https://doi.org/10.1007/s10494-020-00215-0

Google Scholar

Crossref

63.

F.

Zhang

,

T.

Zirwes

,

T.

Häber

,

H.

Bockhorn

,

D.

Trimis

, and

R.

Suntz

, “

Near wall dynamics of premixed flames

,”

Proc. Combust. Inst.

38

(

2

),

1955

–

1964

(

2020

).

https://doi.org/10.1016/j.proci.2020.06.058

Google Scholar

Crossref

64.

T.

Zirwes

,

Z. F.

Sontheimer

,

A.

Abdelsamie

,

F. E.

Hernández-Pérez

,

O. T.

Stein

,

H. G.

Im

,

A.

Kronenburg

, and

H.

Bockhorn

, “

Assessment of numerical accuracy and parallel performance of OpenFOAM and its reacting flow extension EBIdnsFoam

,”

Flow, Turbul. Combust.

(published online).

https://doi.org/10.1007/s10494-023-00449-8

65.

D.

Goodwin

,

H.

Moffat

, and

R.

Speth

, “

Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics and transport processes version 2.3.0b

,”

2017

.

Google Scholar

66.

R.

Kee

,

M.

Coltrin

, and

P.

Glarborg

,

Chemically Reacting Flow: Theory and Practice

(

Wiley

, London,

2005

).

Google Scholar

67.

T.

Zirwes

,

F.

Zhang

,

P.

Habisreuther

,

M.

Hansinger

,

H.

Bockhorn

,

M.

Pfitzner

, and

D.

Trimis

, “

Quasi-DNS dataset of a piloted flame with inhomogeneous inlet conditions

,”

Flow, Turbul. Combust.

104

,

997

–

1027

(

2020

).

https://doi.org/10.1007/s10494-019-00081-5

Google Scholar

Crossref

68.

T.

Poinsot

and

D.

Veynante

, “

Introduction to turbulent combustion

,” in

Theoretical and Numerical Combustion

(

R.T. Edwards, Inc.

,

2005

), pp.

125

–

181

.

Google Scholar

69.

A. Y.

Poludnenko

and

E. S.

Oran

, “

The interaction of high-speed turbulence with flames: Global properties and internal flame structure

,”

Combust. Flame

157

(

5

),

995

–

1011

(

2010

).

https://doi.org/10.1016/j.combustflame.2009.11.018

Google Scholar

Crossref

70.

H.

Helmholtz

, “

On discontinuous movements of fluids

,”

Philos. Mag. Ser.

36

(

4

),

337

–

346

(

1868

).

https://doi.org/10.1080/14786446808640073

Google Scholar

Crossref

71.

W.

Thomson

, “

Hydrokinetic solutions and observations

,”

Philos. Mag.

42

(

281

),

362

–

377

(

1871

).

Google Scholar

Crossref

72.

J. W.

Miles

, “

On the stability of heterogeneous shear flows

,”

J. Fluid Mech.

10

(

4

),

496

–

508

(

1961

).

https://doi.org/10.1017/S0022112061000305

Google Scholar

Crossref

73.

S. A.

Thorpe

, “

Experiments on the instability of stratified shear flows: Miscible fluids

,”

J. Fluid Mech.

46

(

2

),

299

–

319

(

1971

).

https://doi.org/10.1017/S0022112071000557

Google Scholar

Crossref

74.

B. J.

Daly

and

F. H.

Harlow

, “

Transport equations in turbulence

,”

Phys. Fluids

13

,

2634

–

2649

(

1970

).

https://doi.org/10.1063/1.1692845

Google Scholar

Crossref

75.

A.

Attili

and

F.

Bisetti

, “

Statistics and scaling of turbulence in a spatially mixing layer at Reλ = 250

,”

Phys. Fluids

24

(

3

),

035109

(

2012

).

https://doi.org/10.1063/1.3696302

Google Scholar

Crossref

76.

A.

Attili

,

F.

Bisetti

,

M. E.

Mueller

, and

H.

Pitsch

, “

Formation, growth, and transport of soot in a three-dimensional turbulent non-premixed jet flame

,”

Combust. Flame

161

(

7

),

1849

–

1865

(

2014

).

https://doi.org/10.1016/j.combustflame.2014.01.008

Google Scholar

Crossref

77.

H.

Kawashima

,

K.

Kobayashi

, and

T.

Tomita

, “

A combustion instability phenomenon on a LOX/methane subscale combustor

,” AIAA Paper No. 2010-7082,

2010

.

Google Scholar

78.

C. H. K.

Williamson

, “

Defining a universal and continuous Strouhal-Reynolds number relationship for the laminar vortex shedding of a circular cylinder

,”

Phys. Fluids

31

(

10

),

2742

–

2744

(

1988

).

https://doi.org/10.1063/1.866978

Google Scholar

Crossref

79.

A.

Roshko

, “

Perspectives on bluff body aerodynamics

,”

J. Wind Eng. Ind. Aerodyn.

49

(

1–3

),

79

–

100

(

1993

).

https://doi.org/10.1016/0167-6105(93)90007-B

Google Scholar

Crossref

80.

C.

Norber

, “

Flow around a circular cylinder: Aspects of fluctuating lift

,”

J. Fluids Struct.

15

(

3–4

),

459

–

469

(

2001

).

https://doi.org/10.1006/jfls.2000.0367

Google Scholar

Crossref

81.

J. H.

Lienhard

,

Synopsis of Lift, Drag, and Vortex Frequency Data for Rigid Circular Cylinders

(

Technical Extension Service, Washington State University

,

Pullman, WA, USA

,

1966

).

Google Scholar

82.

E.

Achenbach

, “

Influence of surface roughness on the cross-flow around a circular cylinder

,”

J. Fluid Mech.

46

(

2

),

321

–

335

(

1971

).

https://doi.org/10.1017/S0022112071000569

Google Scholar

Crossref

83.

V. P.

Zhukov

and

M. F.

Feil

, “

Numerical simulations of the flame of a single coaxial injector

,”

Int. J. Aerosp. Eng.

2017

,

5147606

(

2011

).

https://doi.org/10.1155/2017/5147606(2017)

Google Scholar

84.

M.

Hansinger

,

P.

Breda

,

J.

Zips

,

C.

Traxinger

, and

M.

Pfitzner

, “

Hybrid RANS/LES simulation of a GOX/GCH₄ 7-element rocket combustor using a non-adiabatic flamelet method

,” Report No. SFB-SS2017-12, Sonderforschungsbereich/Transregio 40, Munich,

2017

.

Google Scholar

85.

N.

Perakis

,

O. J.

Haidn

, and

M.

Ihme

, “

Investigation of CO recombination in the boundary layer of CH₄/O₂ rocket engines

,”

Proc. Combust. Inst.

38

(

4

),

6403

–

6411

(

2021

).

https://doi.org/10.1016/j.proci.2020.07.080

Google Scholar

Crossref

86.

T.

Seitz

,

A.

Lechtenberg

, and

P.

Gerlinger

, “

Rocket combustion chamber simulations using high-order methods

,” in

Future Space-Transport-System Components under High Thermal and Mechanical Loads

(

Springer

,

2020

), pp.

381

–

394

.

Google Scholar

Crossref

87.

S.

Boulal

,

N.

Fdida

,

L.

Mateuszewski

,

L.

Vingert

, and

M.

Martin-Benito

, “

Flame dynamics of a subscale rocket combustor operating with gaseous methane and gaseous, subcritical or transcritical oxygen

,”

Combust. Flame

242

,

112179

(

2022

).

https://doi.org/10.1016/j.combustflame.2022.112179

Google Scholar

Crossref

88.

P. J.

Davis

, “

Leonhard Euler's integral: A historical profile of the gamma function

,”

Am. Math. Mon.

66

(

10

),

849

–

869

(

1959

).

https://doi.org/10.2307/2309786

Google Scholar

Crossref

Published open access through an agreement withTechnische Informationsbibliothek

Analysis of turbulent mixing in a methane–oxygen recessed injector for space propulsion

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Turbulence characterization

B. Turbulent mixing

C. Gradient assumption

D. Mixing-based filtering

III. SIMULATION SETUP

A. Simulation strategy

B. Numerical solver and resolution requirements

IV. RESULTS

A. Turbulent mixing and flame development

B. Turbulent diffusion dynamics

C. Numerical modeling

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: MODELING OF MEAN FILTERED VALUES

1. Demostration of relationship between Reynolds and Favre averages of mixing filters

2. Modeling of mean filtered values

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

$P$	Pressure	$20 bar$
$T_{u}$	Temperature of the unburnt reactants	$300 K$
$(O / F)$	Global oxidizer-to-fuel ratio	$1.99$
${(O / F)}_{st}$	Stoichiometric oxidizer-to-fuel ratio	$3.9898$
$ϕ_{CH 4}$	Hydraulic diameter of methane injector	$0.05 mm$
$ϕ_{O_{2}}$	Hydraulic diameter of oxygen injector	$0.1 mm$
$ϕ_{H}$	Hydraulic diameter of the entire injection element	$0.3 mm$
$δ_{L 0}$	Laminar flame thickness at stoichiometric conditions	$5.505 μ m$
$s_{L 0}$	Laminar flame speed at stoichiometric conditions	$2.5735 m / s$

Analysis of turbulent mixing in a methane–oxygen recessed injector for space propulsion

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Turbulence characterization

B. Turbulent mixing

C. Gradient assumption

D. Mixing-based filtering

III. SIMULATION SETUP

A. Simulation strategy

B. Numerical solver and resolution requirements

IV. RESULTS

A. Turbulent mixing and flame development

B. Turbulent diffusion dynamics

C. Numerical modeling

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: MODELING OF MEAN FILTERED VALUES

1. Demostration of relationship between Reynolds and Favre averages of mixing filters

2. Modeling of mean filtered values

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only