Numerical Simulation of a GH₂/LOx Single Injector Combustor and the Effect of the Turbulent Schmidt Number

Abstract

A large-eddy simulation (LES) of a gaseous hydrogen/liquid oxygen (GH₂/LO_X) single-injector rocket combustor is performed in this study. The Redlich–Kwong–Peng–Robinson (RK–PR) equation of state is used to simulate the real-gas effect under high-pressure conditions, and the steady laminar flamelet model (SLFM) is implemented to simulate fast chemistry, such as a H₂/O₂ reaction. From the numerical simulation, the characteristics of time-averaged flow and flame fields are obtained, and their relationship with the real-gas effect is investigated. It is possible to investigate unsteady flame features and the mixing mechanism of propellants in detail by examining multiple snapshots of the field contour. Another purpose of the study is to investigate the differences in flow and flame structures according to the variation in the turbulent Schmidt number. By comparing the simulation result with the natural OH* emission image and temperature profiles from experimental data, the appropriate range of the turbulent Schmidt number for the simulation is obtained. Furthermore, this paper suggests the usefulness and validity of the current research by quantitatively comparing (i.e., temperature profiles) numerical results with those of existing literature.

1. Introduction

The conventional approach to increasing the propulsion performance of a liquid rocket engine for a launch vehicle is to increase combustion pressure. Achieving a high combustion pressure requires a supply of high-pressure propellant, which inevitably accompanies the liquid or supercritical states of the propellant. Many experiments have been conducted to visualize the behavior of real-gas fluid and the combustion process and to quantify various parameters in these high-pressure conditions. However, when fuel such as hydrogen is used as a propellant, quantifying the parameters is challenging due to the high combustion temperature (e.g., the insertion type sensor is susceptible to high temperatures), and the visualization technique becomes limited (e.g., shadowgraph, radical emission image). In this case, a numerical simulation can be an alternative and efficient analysis tool.

One of the most widely studied rocket combustors, which uses hydrogen and oxygen as propellants and operates in real-gas conditions, is the ONERA RCM-3 combustor. Juniper et al. [1] measured the natural OH* emission using the laser-induced fluorescence (LIF) technique for the A60 and C60 cases of the experiment, conducted at 60 bar. The longitudinal cross-section of the three-dimensional (3D) flame’s structure was restored to a two-dimensional image using Abel inversion. Furthermore, the oxygen jet’s structure was captured using a backlighting technique. Habiballah et al. [2] attempted to measure the flame temperature with the coherent anti-Stokes Raman scattering (CARS) technique for several cases of the experiment and obtained radial-temperature profiles at several axial locations.

Numerical simulations for the combustor were also performed. Poschner and Pfitzner [3] used the commercial code ANSYS-CFX for the simulation, and they used the eddy dissipation model (EDM) and the representative interactive flamelet (RIF) approaches as combustion models in an axisymmetric grid system with 340,000 cells. Their chemical reaction mechanisms included one-step chemistry and a detailed reaction mechanism containing five additional intermediate chemical species. The purpose of their study was to compare the temperature profile and OH field with experimental data and briefly investigate the turbulent Schmidt number effect on the flame field. Kim et al. [4] conducted a numerical simulation based on the steady laminar flamelet model (SLFM) with an axisymmetric grid system. A detailed reaction mechanism and Reynolds-Averaged Navier-Stokes (RANS) turbulence model were used for their study. By comparing their numerical result with the experimental data, it is found that the RANS result underestimated the temperature far downstream of injector since the RANS analysis did not capture the turbulent diffusion. The differences between the ideal-gas equation and the real-gas equation on the combustion field were compared, and the temperature profile and OH field were compared with experimental data.

Coclite and Cutrone [5] conducted a numerical simulation based on the flamelet progress-variable (FPV) approach in an axisymmetric grid system and compared their result with experimental data for the temperature profile and OH field using four different detailed H₂-O₂ reaction mechanisms—the difference was significant. Benmansour et al. [6] conducted their simulation for one-quarter of the combustor volume in a three-dimensional (3D) grid system with more than 900,000 cells. They used the EDM as a combustion model with both one-step chemistry and two-step chemistry as reaction mechanisms, and they used two RANS models. Seidl et al. [7] performed their simulation with the finite-rate chemistry (FRC) approach and a detailed reaction mechanism. Schmitt et al. [8] performed their simulation based on the equilibrium combustion model and large-eddy simulation (LES) in a 3D grid system using 11,000,000 cells. However, they performed only qualitative analysis on the combustion field. Most recently, Riedmann et al. [9] performed their simulation with three different reacting flow solvers and three different combustion models in a 3D grid system.

Of these studies, the axisymmetric simulations [3,4,5,7] and RANS studies [3,4,5,6,7,9] are not sufficient to capture the actual structure and dynamics of a 3D turbulent flame in the combustor. Furthermore, the use of the equilibrium combustion model [8,9] is not sufficient in that the model does not fully reflect the combustion process affected by the strain effect of the flow field. Moreover, although the variation in flame length and structure caused by the change in the turbulent Schmidt number is apparent, previous studies have excluded this number’s effect or a sufficient reason for selecting its value [4,5,6,7,8,9].

The turbulent Schmidt number is defined as the ratio between turbulent momentum diffusion and turbulent molecular diffusion. The change in this number determines whether a velocity shear layer, or, in contrast, a chemical reaction layer, develops predominantly in a flow, thus affecting flame shape. The appropriate value of this number for a jet flow without combustion is usually in the range 0.6–0.8 [10,11], and, for a combustion or gas turbine engine, 0.2–0.5 [12,13]. Poschner and Pfitzner [3] found that, in an aircraft turbine engine simulation, the temperature field can sometimes be more accurately predicted when the value is in the range of 0.3–0.5. Jiang and Campbell [14] found that, in a generic combustor, the range 0.4–0.5 produces a more accurate temperature profile inside the combustor. Ivancic et al. [15] measured the wall heat flux of a gaseous hydrogen/gaseous oxygen (GH₂/GO₂) rocket combustor in their simulation by changing the turbulent Prandtl number and turbulent Schmidt number from 0.6 to 1.1, respectively. However, none of these previous studies [10,11,12,13,14,15] reflected the real-gas effect.

In this study, numerical simulations of a gaseous hydrogen/liquid oxygen (GH₂/LOx) shear-coaxial injector combustor are performed. This study analyzes the dynamics of a turbulent flame by adopting a 3D grid system, the SLFM, and LES. After obtaining a reliable time-averaged field, the study analyzes the combustion flow field reflecting a real-gas effect. The Redlich–Kwong–Peng–Robinson (RK–PR) equation is used to reflect the real-gas effect in the simulation. Furthermore, the characteristics of the combustion flow field according to the change in the turbulent Schmidt number are presented.

A summary of the previous numerical studies is presented in

Table 1. Comparison of numerical studies for Rocket Combustor Modelling (RCM)-3 test case.

	CFD Solver	Grid System	Number of Cells	Combustion Model	Equation of State	$\mathit{S}{\mathit{c}}_{\mathit{t}}(\mathit{P}{\mathit{r}}_{\mathit{t}}$$)\mathbf{or}\mathit{L}\mathit{e}$	Turbulence Model	RCM-3 Case	Ref.
Poschner and Pfitzner	ANSYS CFX	Axi.	340,000	EDM, RIF	RK, PR	$S{c}_t(P{r}_t$) = 0.5, 0.7, 0.8, 0.9	Standard k-ε	A60	[3]
Kim et al.	in-house	Axi.	35,560	SLFM	SRK	Le = 1	RANS	C60	[4]
Coclite and Cutrone	in-house	Axi.	18,000	FPV	PR	Le = 1	Low Reynolds number k-ω	A60	[5]
Benmansour et al.	ANSYS Fluent	3d	794,270–915,349 (AMR)	EDM	polynomial fitting of NIST data	Le = 1	Standard k-ε, Realizable k-ε	A60	[6]
Seidl et al.	TASCOM3D	Axi.	63,500	FRC	SRK-C	N/A	q-ω, k-ω SST	A60	[7]
Schmitt et al.	AVBP	3d	11,000,000	Equilibrium	SRK	$P{r}_t(S{c}_t$) = 0.6	LES(WALE)	A60, C60	[8]
Riedmann et al.	Airbus DS Rocflam3, DLR TAU, ANSYS CFX	Axi., 3d	1,600,000 or less	Equilibrium, FRC, CFX Flamelet	BWR, real-gas data from UDF	N/A	Standard k-ε Spalart–Allmaras, k-ω SST	A10, A60	[9]

.

2. Materials and Methods

2.1. RCM-3 Combustor

The RCM-3 combustor operating on the ONERA’s Mascotte test bench is a small-scale liquid rocket combustor employing a single shear-coaxial injector [2]. The combustor simply consists of an injector and a combustion chamber, as depicted in Figure 1. In the experiment [2], the test bench is equipped with quartz windows, optical diagnostic instruments for visualization, an exhaust nozzle, a propellant supply system, and the combustor. The combustor employs a shear-coaxial type injector: oxygen is injected from the circular orifice of the injector in a compressible liquid state, and hydrogen is injected from the annular orifice in a supercritical state. The propellant supply conditions for the combustor, such as mass flow rate, pressure, and temperature, are summarized in

Table 2. Propellant supply conditions for RCM-3 C60 case.

${\dot{\mathit{m}}}_{{\mathbf{H}}_2}$ [g/s]	${\dot{\mathit{m}}}_{{\mathbf{O}}_2}$ [g/s]	${\mathit{T}}_{{\mathbf{H}}_2}$ [K]	${\mathit{T}}_{{\mathbf{O}}_2}$ [K]	${\mathit{p}}_{{\mathbf{H}}_2}$ [bar]	${\mathit{p}}_{{\mathbf{O}}_2}$ [bar]	${\mathit{\rho }}_{{\mathbf{H}}_2}$ [kg/m3]	${\mathit{\rho }}_{{\mathbf{O}}_2}$ [kg/m3]	${\mathit{U}}_{{\mathbf{H}}_2}[\mathbf{m}/\mathbf{s}]$	${\mathit{U}}_{{\mathbf{O}}_2}[\mathbf{m}/\mathbf{s}]$
70	100	287	85	60	60	5.51	1177.8	236	4.35

.

2.2. Numerical Setup

2.2.1. Governing Equations

The objective of this study is to perform a 3D LES of turbulent combustion based on the SLFM. Therefore, the governing equations for the numerical simulation are the continuity, momentum, mixture fraction, and mixture fraction variance equations as written in Equations (1)–(4), respectively. All equations are filtered by Favre filtering and expressed in Einstein notation with $i$,$j$ = 1–3. The turbulent Schmidt number in the diffusion terms of Equations (3) and (4) is defined as $S{c}_t$ = ${\nu }_t$/$D_t$.

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial x_i}=0$$

(1)

$$\frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=\frac{\partial \overline{p}}{\partial x_j}+\frac{\partial {\overline{\tau }}_{ij}}{\partial x_j}-\frac{\partial {\tau }_{ij}^{sgs}}{\partial x_j}$$

(2)

$$\frac{\partial \overline{\rho }\tilde{z}}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_j\tilde{z}}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\frac{{\mu }_t}{S{c}_t}\frac{\partial \tilde{z}}{\partial x_j}\right)$$

(3)

$$\frac{\partial \overline{\rho }{{\tilde{z}}^{″}}^2}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_j{{\tilde{z}}^{″}}^2}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\frac{{\mu }_t}{S{c}_t}\frac{\partial {{\tilde{z}}^{″}}^2}{\partial x_j}\right)+\frac{2{\mu }_t}{S{c}_t}{\left(\frac{\partial \tilde{z}}{\partial x_j}\right)}^2-\overline{\rho }\tilde{\chi }$$

(4)

2.2.2. Sub-Grid Scale Model

The sub-grid scale (SGS) stress tensor ${\tau }_{ij}^{sgs}$ in Equation (2) is modeled by the Smagorinsky turbulence model [16]. The Smagorinsky model assumes ${\tau }_{ij}^{sgs}$ based on the Boussinesq’s eddy viscosity hypothesis, as written in Equation (5). The term ${\mu }_{sgs}$ is expressed as a function of $\Delta$ and $k_{sgs}$ as in Equation (6) to more accurately represent the SGS eddy length scale and the velocity scale, where $C_k$ = 0.094. $k_{sgs}$ is expressed as Equation (7) based on the assumption that the generation and destruction of SGS eddy energy are in equilibrium—the local equilibrium hypothesis. The coefficients in Equation (7) are $A$ = $C_{ϵ}$/$\Delta$, $B$ = 2tr(${\tilde{S}}_{ij}$)/3 and $C$ = 2$C_k\Delta$[${\tilde{S}}_{ij}$–tr(${\tilde{S}}_{ij}){\delta }_{ij}$/3], respectively. Constant $C_{ϵ}$ = 1.048, and the operator tr is the trace of a second-rank tensor.

$${\tau }_{ij}^{sgs}=-2{\mu }_{sgs}{\tilde{S}}_{ij}+\frac{2}{3}{k}_{sgs}{\delta }_{ij}$$

(5)

$${\mu }_{sgs}=\rho C_k\Delta \sqrt{k_{sgs}}$$

(6)

$$k_{sgs}={\left(-\frac{B+\sqrt{B^2+4AC}}{2A}\right)}^2$$

(7)

The Smagorinsky model generally requires a near-wall treatment model to damp out the eddy viscosity near the wall of the domain, because it tends to overestimate eddy viscosity near the wall. Therefore, the van Driest model [17] is applied to the filter cutoff length $\Delta$ near the wall, as written in Equation (8).

$${\Delta }_{vD}=\min \left[\Delta ,\frac{\kappa }{C_{\Delta }}y\left(1-e^{-y^{+}/A^{+}}\right)\right]$$

(8)

where the von Kármán constant $\kappa$= 0.4187, $C_{\Delta }$ = 0.158, and $A^{+}$ = 26. The limit of $y^{+}$ is 0–500.

2.2.3. Equation of State

$$p=\frac{RT}{v-b}-\frac{a_c{\left(\frac{3}{2+T_r}\right)}^n}{\left(v+{\delta }_{1}b\right)\left(v+\frac{1-{\delta }_1}{1+{\delta }_1}b\right)}$$

(9)

The real-gas equation used in this study is the RK–PR equation proposed by Cismondi and Mollerup [18]. The RK–PR possesses the merits of the Soave–Redlich–Kwong (SRK) equation [19] and the Peng–Robinson (PR) equation [20], enabling it to accurately predict the fluid properties of both light molecules and low hydrocarbon species. Although the PR equation is proposed to more accurately predict the vapor pressure of the low hydrocarbon species, the equation performs poorly for particular chemical species in predicting a critical compressibility factor as the pressure increases to close to 1000 bar. Moreover, for hydrocarbons less than C4, the PR equation performs poorly in assessing the density compared to the SRK equation. In these scenarios, the RK–PR equation is a superior model.

Cismondi and Mollerup [18] proposed a fairly complex form of the equation, as defined in Equation (9), while curve fitting their cubic equation of state from existing literature. The numerator of the second term of the right-hand side of the equation is proposed for improving the vapor pressure assessment. The denominator of the second term is proposed to combine the SRK equation and PR equation and to reduce the density deviation. The coefficients $a_c$, $n$, ${\delta }_1$, and $b$ are described in detail in the literature.

In this study, before the numerical simulation is performed, the distributions of two parameters, i.e., density and isobaric specific heat of pure oxygen, according to the temperature change in high-pressure condition (60 bar) are calculated, and compared between the SRK, PR, and RK–PR equations. As a result, the PR equation had the least accurate prediction, while the SRK and RK–PR equations had similar distributions for the two parameters close to the National Institute of Standards and Technology (NIST) data [21]. This result is consistent with the result of Kang et al. [22], who studied cryogenic nitrogen under supercritical condition. They also compared the two thermodynamic properties of the three real-gas equations to the NIST data. Consequently, the SRK and RK–PR demonstrated similar performance.

In this study, the RK–PR equation is eventually selected, considering its wide range of applicability from light molecule fuel to hydrocarbon fuel. As shown in Figure 2, the distributions of density and isobaric specific heat of the RK–PR equation do not significantly differ from the NIST data except for the range 100–160 K, with a slight deviation.

2.2.4. Combustion Model

The SLFM assumes and models a turbulent flame as consisting of an ensemble of stretched laminar flamelets [23]. This approach is suitable for fast chemistry, such as a H₂-O₂ reaction, and it accurately reflects the strain effect of fluid flow. Moreover, for liquid rocket propulsion, combustion typically belongs to the flamelet regime [24,25]. Therefore, this approach is used as the combustion model for this study. The distribution of scalar properties, such as density, temperature, and chemical species mass fraction according to $\tilde{z}$, can be tabulated by numerical analysis of a one-dimensional laminar diffusion flame before the simulation for the RCM-3. In the tabulation process, a probability density function (pdf) is required to calculate the mean quantities of the fluctuating properties in a turbulent flow. For an axisymmetric jet flow, particularly, the pdf of the beta distribution ($\beta$-pdf) can be used. The $\beta$-pdf used in this study is written in Equations (10) and (11) [26]. Therefore, the turbulent scalar properties can be obtained from the integration of the pdf, as expressed in Equation (12).

$$\tilde{P}\left(z\right)=\frac{z^{\alpha -1}{\left(1-z\right)}^{\beta -1}}{\mathsf{\Gamma }\left(\alpha \right)\mathsf{\Gamma }\left(\beta \right)}\mathsf{\Gamma }\left(\alpha +\beta \right)$$

(10)

$$\alpha =\tilde{z}\gamma ,\beta =\left(1-\tilde{z}\right)\gamma ,\gamma =\frac{\tilde{z}\left(1-\tilde{z}\right)}{\widetilde{{z^{″}}^2}}$$

(11)

$$\tilde{\varphi }={\displaystyle \iint }\varphi \left(z,{\chi }_{st}\right)\tilde{P}\left(z\right)\tilde{P}\left({\chi }_{st}\right)dzd{\chi }_{st}$$

(12)

The scalar transport equations for $\tilde{z}$ and $\widetilde{{z^{″}}^2}$ are given in Equations (3) and (4), respectively. The scalar dissipation rate (SDR) $\tilde{\chi }$ in Equation (4) is defined as Equation (13) [27].

$$\tilde{\chi }=C_{\chi }\frac{\mu +{\mu }_t}{\overline{\rho }S{c}_t}{\left(\frac{\partial \tilde{z}}{\partial x_j}\right)}^2$$

(13)

where $C_{\chi }$ = 2.

The H₂-O₂ reaction mechanism used in the calculation of a one-dimensional laminar diffusion flame is obtained from Conaire et al. [28]. This mechanism contains 11 chemical species and 19 reaction steps. In this study, eight chemical species (H₂, O₂, H, O, H₂O, OH, HO₂, and H₂O₂) are used, three chemical species (N₂, Ar, and He) are excluded, and 12 reaction steps are used. This mechanism is rigorously tested by comparing with experimental data from existing literature and with conventional H₂-O₂ reaction mechanisms at an initial temperature range of 298–2700 K, a pressure range of 0.05–87 atm, and an equivalence ratio range of 0.2–6. Conaire et al.’s mechanism proposes superior performance under high-pressure condition exceeding tens of atmospheric pressure. For example, this mechanism functioned accurately in a numerical simulation of a GH₂-LOx rocket combustor under high pressure with cryogenic condition [7,29].

2.2.5. Computational Grid

Figure 3 illustrates the computational grid of the combustor. The propellant supply pipes are attached to the chamber to simulate the hydrodynamically developed velocity profile at the injector exit. The entire length of the computational domain is approximately 44 cm, and the diameter of the combustion chamber is 5 cm. The grid system of the combustor consists predominantly of hexahedral elements with a hint of prism elements. The average non-orthogonality of the grid system is 4°, such that the artificial diffusion of scalar properties is minimized. The grids near the wall are refined to capture the boundary layer. The number of cells reaches approximately 3,000,000, and the average size of the cells except for those in the boundary layer is approximately 0.7 mm. Ansys Meshing software is used as a grid generation tool for the present study. Grid convergence test is performed using three levels of grids, coarse, medium, and fine. Each grid has 890,000, 3,000,000, and 10,125,000 cells, respectively. Grid size is reduced by 2/3 in a one-dimensional sense, as the number of grids increase. Figure 4 is the comparison of time-averaged temperature and density variations along the center line of the combustor. The coarse grid has large difference at near field of injector and slightly over predicts the temperature at downstream. However, the difference between medium and fine grids is not so big in an entire region and agrees quite well overall. Therefore, the medium grid is selected for further studies by considering the computational efficiency.

The fast Fourier transform (FFT) result of the turbulent kinetic energy (TKE) measured at a probe point ($x$/$d_2$ = 1) in Figure 3 throughout the entire simulation time is shown in Figure 5. The FFT result illustrates that the grid system satisfies the Kolmogorov’s 5/3 law, such that the grid system has sufficient resolution to resolve various turbulence length scales.

2.2.6. Solver Setup

In this study, SLFMFoam_LES implemented on OpenFOAM, an open source platform, is used as the main reacting flow solver. This solver is developed based on both the SLFM and the PIMPLE algorithm, focusing on an LES of turbulent combustion. PIMPLE combines the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) steady-state flow algorithm [30] with the PISO (Pressure Implicit with Splitting of Operator) unsteady flow algorithm [31]. The detailed structure of SLFMFoam_LES is depicted in Figure 6. The solver has a structure that solves momentum, mixture fraction, and mixture fraction variance equations in order, and then calculates SDR. Then, the solver performs linear interpolations on the mixture fraction, mixture fraction variance, and SDR to obtain the temperature and density from the flamelet table. Because this solver is built upon the pressure-based segregated algorithm, as depicted in Figure 6, it solves the pressure equation based on the number of inner corrections and repeats the outer correction according to the nature of the SIMPLE algorithm. If the convergence limit is reached, the calculation moves on to the next time step.

For temporal discretization, the implicit Euler scheme with first-order accuracy is applied to the transient term of each equation. For spatial discretization, the minmod limiter [32] with second-order accuracy is applied to the convective term of each equation. The diffusion term of each equation uses the central-differencing scheme with second-order accuracy. For the convergence check, the absolute residual for the flow solver is set to 1 × 10⁻⁵ for all equations and this criterion is achieved during the calculations.

The boundary conditions are as follows. A fixed mass flow rate condition is imposed on both the hydrogen and oxygen inlets. A non-slip, adiabatic condition is applied to all walls, and a fixed pressure condition is set at the outlet. The van Driest model is used for the boundary layer at the walls. The turbulence intensities at the inlets are set to zero, because the experimental data on turbulence property does not exist.

The range of SDR constructed on the flamelet table is 105–10⁵. Figure 7 illustrates two scalar variables, i.e., temperature and density, according to the mixture fraction obtained from the table construction. Figure 7a is the temperature graph and it illustrates that the stoichiometric mixture fraction value is positioned near $z_{\mathrm{st}}$ = 0.12. At this location, the temperature reaches close to the constant-pressure adiabatic flame temperature of 3627 K. Figure 7b is the density graph and it illustrates a significant density variation near $z$ = 0–0.004. This mixture fraction range is where isobaric specific heat displays a singularity—as heat energy is generated, the fluid experiences a volumetric expansion rather than a rise in temperature in this range.

The simulation time is set to 1 s, which is sufficient for the flow to pass through the combustion chamber’s outlet more than 70 times and to calculate a reliable time-averaged flow field. The time step size is fixed at 5 × 10⁻⁷ s, which corresponds to the Courant number of 3. For the computation, 612 cores of an Intel Xeon Gold 6154 3GHz processor are used for the simulation, which corresponds to 58.75 teraFLOPS. The wall clock time is 67 h, which corresponds to 41,004 CPU-hr.

3. Results

3.1. Flow and Flame Field Characteristics

Based on the results, the time-averaged flow and flame structures are first investigated by examining the contour fields obtained from the simulation. Then, to observe the unsteady characteristics of the flow and flame, multiple contour snapshots are investigated in sequential order in time. Then, to observe the effect of the turbulent Schmidt number, a series of numerical simulations is conducted by increasing the turbulent Schmidt number’s interval by 0.1 in the range 0.2 to 0.8. The simulation result is validated by comparing with the experimental OH* emission image. In the latter section of the results, the temperature profiles are quantitatively compared between this study and the existing literature.

Figure 8 shows the longitudinal slices of the time-averaged temperature, velocity, and mixture fraction for the combustor. In the temperature contour, the iso-line of $\tilde{z}$ = 0.004 (white) is plotted along the edge of the oxygen jet core. From this line, the inside is the liquid state, and the outside is the supercritical state. In the temperature field in Figure 8a, the maximum time-averaged temperature of the flame reaches approximately 2400 K, occurring predominantly in the downstream region after the location at $x$/$d_2$ = 10. This location, $x$/$d_2$ = 10, is where the flame starts to expand rapidly. The flame, having passed through the transcritical region ($\tilde{z}$ = 0–0.004), and entered the supercritical state ($\tilde{z}$ > 0.004, similar state to gas), remains in the form of an elongated flame starting from the injector lip to the location of flame expansion ($x$/$d_2$ = 10). This elongation is caused by the low degree of mixing between the hydrogen and oxygen due to the low turbulent intensity between the two jets—even though the velocity of the hydrogen jet is very high. As the flow reaches the reattachment point behind the annular recirculation zone (ARZ), depicted in Figure 8b, the velocity of hydrogen becomes very low, and the hydrogen–oxygen mixing becomes more active, causing the flame to expand rapidly in the radial direction.

The time-averaged velocity field in Figure 8b illustrates that a very-high-speed hydrogen gas with a maximum time-averaged velocity of 340 m/s is injected from the annular orifice of the injector. Furthermore, the path of the jet is gradually curved in the radial direction toward the wall as the path is affected by the flame expansion. Near the reattachment point, the velocity of the jet is significantly slower and close to 0 m/s.

The mixture fraction field in Figure 8c illustrates that because the mixing between the hydrogen and oxygen jets is not active near the injector exit, the gradient of the mixture fraction at that location is steep and distinct. However, after the location of flame expansion, the intensity of mixing increases, and, finally, the mixture fraction becomes homogenous in the latter part of the combustor ($x/d_2$ > 30). Simultaneously, as illustrated in Figure 8a, the high-temperature flame (2400 K) gradually transforms into a medium-temperature flame (1500 K), before and after the location $x$/$d_2$ = 30.

3.2. Unsteady Flow and Flame Features

Figure 9 illustrates the process of mixing between the oxygen and hydrogen jets as the jets flow downstream starting from the injectors. This mixing is facilitated by two effects. The group of contours in Figure 9 illustrates the density field in log scale with the heat release rate overlapped on the density field in grayscale. After the specific points of simulation time, $t_0$ and $t_1$, are selected, the flow-field change according to uniform $\Delta t$ change from $t_0$ and $t_1$ is investigated.

Figure 9a illustrates the Kelvin–Helmholtz (K–H) instability is developing along the surface of the oxygen jet core due to the velocity difference between the high-speed hydrogen jet and the low-speed oxygen jet discharged from the injectors at $t_0$ ms. At $t_0$ + 1 ms, a portion of the oxygen jet is pressed down by the strong curling effect of the first penetrating-vortex, and at $t_0$ + 2 ms, the tip of the jet eventually breaks up by the sustained curling pressure from the same vortex. When breaking up, the mixing of the propellants intensifies, and the maximum heat release rate occurs at the jet separation region. Simultaneously, the second penetrating-vortex develops and moves down downstream. At $t_0$ + 3 ms, the second penetrating-vortex presses down the oxygen jet once more. This compression again produces a strong heat-release rate at the separation location. The tip of the jet splits once more at $t_0$ + 4 ms.

In Figure 9b, the intact oxygen jet flows at t₁ ms, and the recirculation zone described in Figure 8b is placed outside the hydrogen jet. At $t_0$ + 2 ms, the lower part of the oxygen jet begins to be sucked toward the lower reattachment point by the lower recirculation flow. The sucked portion of the jet reacts with the recirculation flow to produce an intermediate heat release rate. At $t_0$ + 4 ms, the upper part of the oxygen jet is weakly sucked toward the upper reattachment point by the upper recirculation flow. At $t_0$ + 6 ms, the oxygen jet is split by the combination of the recirculation and the K–H instability. At $t_0$ + 8 ms, the oxygen jet is returned to its intact state.

Therefore, strong heat generation is caused predominantly by propellant mixing due to the K–H instability and the recirculation flow effects, occurring primarily upstream ($x$/$d_2$ < 12) of the combustor. As depicted in Figure 9, in the region behind this strong heat generation, i.e., $x$/$d_2$ = 12–18, an effective mixing mechanism does not work, so the mixing rate there is the slowest, and the heat generation is the lowest in the combustor.

3.3. Effect of Turbulent Schmidt Number on Flow Field

Figure 10 shows the longitudinal slices of the time-averaged temperature field overlapped with the streamline for the turbulent Schmidt number in the range 0.2–0.8 with an interval 0.2, and the white iso-line corresponding to the time-averaged mixture fraction of 0.2 is plotted for each slice. This iso-line is the edge of the expansion zone of oxygen plotted in Figure 10, and this expansion zone (i.e., the inside of the white iso-line) is where the oxygen jet core suddenly begins to expand and stops expanding.

The exact location of the annular recirculation zone described in Section 3.1 is depicted. The annular recirculation zone consists of two smaller annular recirculation zones due to the protrusion of the injector, and its function is to pre-heat the cold oxygen and hydrogen jets throughout the combustion, thus contributing to flame stabilization in the combustor [33].

The change in the flow and flame structures with the increasing turbulent Schmidt number is indistinguishable between the four slices except for $S{c}_t$ = 0.2. For $S{c}_t$ = 0.2, the streamline is curved in the radial direction earlier than for other $S{c}_t$ cases. Simultaneously, the flame is located further upstream, and the flame is thicker than the others. As soon as the propellant is injected from the injector, the mixing starts faster due to the strong turbulent molecular diffusion, leading to the development of flame closer to the injector. However, the development of turbulent momentum diffusion is slower for the $S{c}_t$ = 0.2 case; consequently, it is difficult for oxygen to permeate into the annular recirculation zone, producing a thicker flame.

The remaining three slices are different from $S{c}_t$ = 0.2 in the flow and flame structures, and they have apparently similar structures. However, the streamline illustrates that, for $S{c}_t$ = 0.4, the third small recirculation zone exists, whereas in the cases of $S{c}_t$ = 0.6 to 0.8, the third small recirculation disappears. The generation of this third small recirculation is related to the size of the expansion zone of the oxygen jet: the size of the third recirculation zone increases with that of the expansion zone. The expansion zone’s size decreases as the turbulent Schmidt number increases. Another notable turbulent Schmidt number effect is the streamline pattern in the downstream region after the expansion location ($x$/$d_2$ > 10). As the turbulent Schmidt number increases, the streamline is less parallel to the $x$ -axis and is curved inwards. As the turbulent Schmidt number increases, the specific volume behind the expansion location decreases more drastically due to the strong mixing between the two propellants, i.e., unburnt hydrogen gas flowing toward the outlet via the reattachment point and the oxygen in the expansion zone. Figure 11 supports this fact by showing the kinematic turbulent viscosity is the highest at $S{c}_t$ = 0.8 over the flame expansion zone painted in bright blue ($x$/$d_2$ = 10–26).

3.4. Effect of Turbulent Schmidt Number on Flame Field

Figure 12 compares the time-averaged OH field of the simulation to the experimental OH image of Juniper et al. [1] near the injector exit according to the change in the turbulent Schmidt number. The experimental image is the Abel transform of the natural OH* emission image. This image also includes a backlighting image for the oxygen jet in the range $x$/$d_2$ = 0–13 to confirm the appearance of the oxygen jet. Simultaneously, the contour of the half-thickness of the oxygen jet for this study is overlapped on its OH field in grayscale, ranging from the minimum value of $\tilde{z}$ = 0 (black) to the maximum value of ${\tilde{z}}_{st}$ = 0.12 (white). Both the simulation and experiment OH results are non-dimensionalized based on their respective maximum OH values, and the non-dimensionalized legends are also placed on each corner.

The flame of the simulation for $S{c}_t$ = 0.2 differs significantly from that of the experimental image. The OH distribution of $S{c}_t$ = 0.2 is noticeably wider and thicker than the other three cases ($S{c}_t$ = 0.4, 0.6, 0.8). The half-thickness of the oxygen jet is also significantly thicker than that of the experimental image.

In contrast, the flames for $S{c}_t$ = 0.4–0.8 are similar. The intense OH region (bright red in the experiment) is distributed along the elongated region of the flame starting from the injector lip to roughly $x$/$d_2$ = 10, whereas the intense region of the simulation is distributed throughout the area after the flame expansion zone ($x$/$d_2$ = 10). The simulation results of Schmitt et al. [8] and Kim et al. [4] indicate that, like in this study, the intense region is also distributed in the regions after the expansion zone. This difference is thought to originate from the assumption of a constant turbulent Schmidt number. According to studies by Pohl et al. [34] and Xiao et al. [35], assuming a variable turbulent Schmidt number, the turbulent Schmidt number is variable in the region right after an injector lip or at the region where the fuel and oxidizer first start to mix. It is believed that the OH value for the present simulation in this region should be as high as the experimental value because the turbulent Schmidt number in this region is predominantly low, so the molecular diffusion by eddy is predominant, thus increasing propellant mixing.

For the flame spreading angle (the red dotted line in Figure 12) after the expansion location ($x$/$d_2$ = 10), a rapid change in the spreading angle is distinctly visible, and the angle is high in this simulation result, but the rapid change is ambiguous, and the angle is small in the experimental image (11.1°). Three simulation cases have nearly the same flame spreading angles after expansion (31.3°, 32.6°, and 33.5° for $S{c}_t$ = 0.4, 0.6, and 0.8, respectively).

The flame closure in the downstream region is common to both the experiment and the simulation, as depicted in Figure 12. The iso-line of the inner flame surface at $S{c}_t$ = 0.4 corresponds most closely with the experimental image, but the iso-line gradually deviates from the iso-line of the experiment as the turbulent Schmidt number increases, illustrating a significant difference compared to the experiment at $S{c}_t$ = 0.8. The iso-line of the flame surface is created from the iso-value 0.425 in the normalized time-averaged OH range (0 to 1). As the location of the flame closure moves upstream, the length of the time-averaged oxygen jet shortens accordingly.

Figure 13 is a close-up view of the OH field near the injector. The contour lines of the time-averaged jet position are depicted in the experimental OH image. In the experiment image, the central gray region is where the jet is continuously present along the line of sight, and the dark gray contour lines are, respectively, 0.875, 0.75, 0.625, 0.5, 0.375, and 0.25 times the maximum jet existence time as it is farther away from the central-axis. The jet’s contour lines are similarly represented in the simulation image. The white contour lines in the simulation image represent the oxygen mass fraction values of 0.875, 0.75, 0.625, 0.5, 0.375, and 0.25 as it moves farther away from the central-axis. At $S{c}_t$ = 0.2, the flame-spreading angle of the simulation image is larger and the distribution of contour lines is wider than those of the experiment. Therefore, the difference is larger compared to the other turbulent Schmidt number cases. At $S{c}_t$ = 0.8, the flame spreading angle is smaller, and the distribution of contour lines is narrower, than those of the experiment. In contrast, $S{c}_t$ = 0.4, 0.6 agree closely with the experiment. More precisely, $S{c}_t$ = 0.4 has the closest resemblance.

Consequently, the appropriate turbulent Schmidt number range for this study is 0.4–0.5, producing the most accurate resemblance to the experimental image from two perspectives: the inner flame surface’s size and the oxygen jet’s contour distribution.

Figure 14 illustrates the time-averaged temperature profiles along the central-axis for the different turbulent Schmidt numbers. As the turbulent Schmidt number increases, the position of the peak temperature (2400 K) moves upstream. The difference of the peak temperature position between the numbers 0.4 and 0.8 reaches $x$/$d_2$ = 3.1, which is the distance that is approximately three times the oxygen injector’s pipe diameter. This phenomenon quantitatively confirms that the size of the oxygen jet expansion zone described in Section 3.2 decreases with an increase in the turbulent Schmidt number. The expansion zone is where most of the oxygen jet breaks up, and the cause of the break-up is that the tip of jet potential core is the most affected by the excitation of the acoustic wave in the combustor and this leads to the largest fluctuation at the tip of jet [36]. Moreover, the presence of the first peak showing an intense turbulent mixing in Figure 11 supports this excitation phenomenon. As the turbulent Schmidt number increases, the mixing rate between the propellants in this expansion zone is further increased by turbulent mixing, which reduces its size.

Figure 15 illustrates the time-averaged temperature profiles—one along the central-axis and the other three along the radial direction at several locations ($x/d_2$ = 3, 10, and 20) for $S{c}_t$ = 0.4—and compares these profiles with those of previous studies. The experimental data are measured using CARS, and each probe point is averaged over 225 samples for 15 s.

In Figure 15, all researchers’ results are from their RANS simulations, except for this study. In Figure 15a, The EDM result of Poschner and Pfitzner [3] reaches a maximum temperature of 3850 K at $x/d_2$ = 30. They use a one-step reaction mechanism that does not contain intermediate chemical species in the combustion model. Benmansour et al. [6] use a one-step reaction mechanism, and Figure 15a illustrates that their peak temperature approaches 3760 K at $x/d_2$ = 19. Due to the generation of a higher TKE of the realizable k-ε turbulence model, the peak temperature is moved further upstream than that of Poschner and Pfitzner [3] and is located adjacent to the third measurement point. The third measurement point is close to the adiabatic flame temperature of 3627 K; thus, the measurement value is believed to be the instantaneous flame temperature rather than the time-averaged flame temperature.

Riedmann et al. [9] also point out that, from the perspective of turbulent flow, the temperature value appears to be somewhat overestimated. Therefore, it is questionable whether the temperature measurement method in the experiment is appropriate. The RIF result of Poschner and Pfitzner [3] includes five intermediate chemical species; consequently, the maximum temperature is lowered significantly to 2820 K, but the overall temperature distribution is still higher than the experiment. The result of the equilibrium model of Riedmann et al. more closely agrees with the experiment due to the use of a 3D grid system. The result of this study exhibit similar performance to those of Riedmann et al., except that the temperature distribution for this study also accurately reflects the second measurement points.

Figure 15b–d are the radial temperature profiles, with only 4–5 measurement points; although this complicates the comparative analysis, the results of this study are superior to those of the other RANS studies in terms of trend and deviation. In Figure 15b, the result of Seidl et al. [7] still illustrates a higher peak temperature. The results of the equilibrium model of Riedmann et al. [9] demonstrate that molecular diffusion is strong near the injector exit, so the flame thickness is wider than in the other studies. In Figure 15c, the RANS results also illustrate higher peak temperatures and wider flame thicknesses. In Figure 15d, the RANS results do not agree with the distribution of experimental points regardless of the combustion model or turbulence model.

Figure 16 illustrates the time-averaged density profile along the central-axis for the near injector field. For examining the length characteristics of the oxygen jet, the liquid length, and jet length of oxygen are defined as follows [37].

$$l_{\rho }^{ic}=x_{\rho }^{0.99}-x_0$$

(14)

$$l_{\rho }^{0.05}=x_{\rho }^{0.05}-x_0$$

(15)

The superscript ic in Equation (14) denotes an intact oxygen jet or pure liquid oxygen. In Equation (14), the liquid length of oxygen is defined as the $x$-coordinate corresponding to 99% of the density at the injector exit minus the $x$-coordinate of the injector exit. In Equation (15), the jet length of oxygen is the $x$-coordinate corresponding to 5% of the density at the injector exit minus the $x$-coordinate of the injector exit.

The fluid densities at the injector exit are $\rho$ = 828 kg/m³ for the oxygen and $\rho$ = 5.4 kg/m³ for the hydrogen, and, therefore, the liquid length of oxygen is calculated to be roughly 0.5, as depicted in Figure 17. Thus, the liquid length of oxygen is as short as half the oxygen injector diameter. This indicates that the oxygen jet rapidly mixes and reacts with the hydrogen jet at the injector lip as soon as it is injected, producing the resultant characteristics of the attached flame. In the range $x$/$d_2$ = 0.5–10 in Figure 16, the density gradient is moderately steep. The smaller the turbulent Schmidt number, the faster the density decreases because the hydrogen–oxygen mixing due to molecular diffusion is stronger than the turbulent diffusion in this range. However, in the range $x$/$d_2$ = 10–20, the effect of turbulent diffusion is stronger, so the larger the turbulent Schmidt number, the faster the density decreases. The jet length of oxygen is calculated as depicted in Figure 17. As described in Section 3.4, as the turbulent Schmidt number increases, the break-up becomes more active at the tip of the oxygen jet, shortening the jet length of oxygen.

4. Conclusions

In this study, a numerical simulation of turbulent combustion for the RCM-3 combustor—a representative case of small liquid rocket combustors with a single GH₂/LOx shear-coaxial injector—is performed. The average flow and flame structures of the combustor can be investigated by temporal-averaging the obtained flow field data, and the dynamics of the flow and flame can be investigated by analyzing the unsteady flow field data. The features of the flow and flame structures according to the change in the turbulent Schmidt number are also investigated.

First, several conclusions are drawn based on the time-averaged flow and flame field analysis. In contrast to the typical elongated flame appearing under ideal-gas condition, the flame under real-gas condition expands rapidly in the radial direction after a specific axial location. This location is where the oxygen jet breaks up intensely, and, simultaneously, the location of the reattachment point is coincident with this location such that the oxygen jet strongly spreads in the radial direction. Moreover, concerning the fluid state, this location is where liquid oxygen changes from a liquid to a supercritical state (close to a gaseous state), in which the specific volume of gas increases rapidly. These phenomena are likely the major causes of the rapid flame expansion.

Second, unsteady flow and flame field analyses are performed. The focus is placed on analyzing the mixing mechanism of the fuel and oxidizer. In the upstream region of the combustor, i.e., before the expansion of the oxygen jet occurs, the main mechanism of mixing is caused primarily by the K–H instability and recirculation flow. The hydrogen jet diffuses into or penetrates the oxygen jet by the interface vortex caused by the K-H instability, resulting in increased mixing and heat release rate. Simultaneously, a part of the oxygen jet is sucked toward the reattachment point and is mixed with recirculation flow, increasing the heat release rate around the suction region. The recirculation flow effect also contributes to the break-up of the oxygen jet.

Third, the flow and flame field feature according to the turbulent Schmidt number change is analyzed. Turbulent Schmidt numbers in the range 0.2–0.8 are analyzed, and the results for the range 0.2–0.3 demonstrate a significant difference in the appearance of the flame when comparing the OH field to the experimental image. The range 0.4–0.5 displays a close resemblance to the experimental image in terms of the inner flame surface’s size. The oxygen jet’s contour distribution for the range 0.4–0.5 also most closely resembles the experimental image. Comparing the streamline and temperature field in the range 0.4–0.8, the expansion zone of the oxygen jet decreases as the turbulent Schmidt number increases. The oxygen jet’s break-up is further increased by strong turbulent mixing.

Finally, the axial temperature profile along the central-axis and the radial temperature profiles at several $x$-coordinates of the simulation results—four temperature profiles in total—are compared with experimental data. In the central-axis temperature profile, when comparing with the RANS results from existing literature, the LES result of this study illustrates the smallest deviation from the experimental data. The remaining three radial profiles also illustrate the most accurate trends and smallest deviations of this study when compared to the RANS results. Therefore, at least a 3D LES is necessary for a numerical simulation of the RCM-3 combustor. The axial temperature difference along the central-axis according to the change in the turbulent Schmidt number is not large, but the profiles for range $S{c}_t$ = 0.4–0.5 have the smallest deviation from the experimental profile.