Development and Validation of a Compressible Reacting Gas-Dynamic Flow Solver for Supersonic Combustion

1. Introduction

Supersonic airbreathing or supersonic combustion ramjet (scramjet) engine is the most promising technology for significant reduction flight times for long-distance civil flights and the first stages of space systems for putting payload into orbit. The future of access-to-space systems is faced with many technical and economic challenges. Designing and testing scramjet systems requires the ability to adequately understand and control the complex flow physics of these systems. The high-speed flow in the propulsion systems leads to many difficulties in achieving and maintaining combustion [1,2]. One of the significant features of the scramjet is the supersonic flow conditions in the combustor lead to very short residence times for fuel and air, typically of the order of milliseconds; therefore, the fuel must mix with the supersonic flow and burn entirely within a short time. Thus, future hypersonic propulsion systems' success will largely depend on efficient injection, mixing, and combustion processes inside the supersonic combustion chamber [1,3].

Computational fluid dynamics (CFD) and experimental testing have seen significant advances through the years in terms of reliability and accuracy. However, in the case of supersonic and hypersonic regimes, realistic high-enthalpy conditions can only be maintained and thus measured in experimental facilities for a very short period, in the order of milliseconds [4] which is insufficiently to provide a comprehensive analysis designing efficient scramjet engines. The numerical approach remains the primary analysis tool to cover the operating regimes of a scramjet consistently. CFD methods targeting the solution of the flow field inside supersonic engines are faced with the challenges of addressing internal flow dominated by complex shock wave patterns, shock wave boundary layer interaction with its separation, and most importantly, supersonic turbulence-combustion interaction, which, in many cases, is characterized by finite rate chemistry. Inaccurate prediction of shock position and separation locations can lead to unreliable results when reactions are activated, risking the success of the entire design process.

The computational approaches to simulate turbulent flows by solving compressible Navier-Stokes equations are Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), and Reynolds Averaged Navier-Stokes (RANS). RANS is very useful in many real-world industrial applications because RANS simulations are computationally less expensive than DNS and LES. In RANS, fewer equations must be solved, and they are often solved with limited computational resources. However, because the turbulent flows using RANS are not fully resolved, turbulence-combustion interactions are also unresolved. Moreover, the reacting flow is highly nonlinear, unstable, and non-stationary. RANS will not adequately model turbulence-combustion interactions in scramjet [5] due to these constraints. The Direct Numerical Simulation (DNS) is a detailed computational approach that demands a fine grid to resolve all turbulent length scales and small time-steps comparable to chemical time scales. Hence, DNS would require tremendous computational resources for the simulations of the high-speed flows of practical interest. On the other hand, Large-Eddy Simulation (LES) resolves large scales of the turbulent flow, i.e., turbulent eddies larger than the smallest mesh size, while the smallest scales are modeled [6,7,8]. Thus, this feature of LES leads to a significant reduction in computational resources compared to DNS. During the past decade, LES, which is a compromise between the DNS [9] and the RANS [10] methods, has been extensively used to study high-speed combustion applications [11,12,13]. In LES, turbulence scales larger than the grid spacing are resolved, while subgrid scales and their effects on the large scales are modeled [14]. LES has been shown [15,16,17,18,19] to be computationally practical for modeling unsteady phenomena, including transient combustion events. One of the critical aspects of modeling unsteady turbulent combustion in high-speed jets is appropriately accounting for the influence of turbulence on the rate of chemical reactions. Transient turbulent combustion will also depend on accurately modeling the interaction between the jet mixing layer and the surrounding air.

DNS and LES of compressible turbulent reactive flows require highly accurate numerical schemes, which must capture shock waves, describe multicomponent mixtures and detailed chemical kinetics, and resolve a broad range of length scales. Numerical simulation of combustion in supersonic flow fundamentally differs from combustion in the subsonic regime. In supersonic combustion, chemical reactions occur under conditions of high pressure, temperature, and density. These conditions can significantly alter chemical kinetics and reaction rates, leading to different reaction pathways and species formation compared to subsonic conditions. Ignition and flame stabilization processes in supersonic combustion devices may differ significantly from those observed in the subsonic regime [3]. Combustion processes involve highly nonlinear chemical reactions and many species. Furthermore, turbulence increases the mixing of fuel and oxidizer and enhances combustion [20]. In turn, the combustion can intensify the turbulence via flow expansion [20]. As a result, the turbulent fluctuations can be very large [21].

Various CFD solvers in the literature are used to simulate high-speed airbreathing engines. These codes can be broadly categorized into three main categories: commercial, open-source, and in-house solvers, each with its own advantages and disadvantages. The commercial codes such as Fluent [22], StarCCM++ [23], NASA codes REACTMB [24], WIND-US [25], TAU-code code [26], VULCAN [27], SU2 [28] have been used successfully in solving the above-mentioned problems. Commercial solvers can be expensive, with limited access to the source code. Open-source solvers are typically freely available, with access to the source code, allowing for extensive customization and the development of user-specific models or algorithms. The open-source community often collaborates on code development. Among the open-source codes, the OpenFOAM library is a well-established open platform in the scientific community [29]. It is noteworthy that there are several in-house codes reported in the literature for high-speed compressible reacting flows [13,14,30,31,32,33,34,35,36,37]. These codes have been used to solve many problems involving turbulence-combustion interactions. Some of them are built on the OpenFOAM foundation. However, none of these codes are readily and freely available, i.e., they are not yet fully open-source [38].

The OpenFOAM library has high-fidelity reactive flow solvers (e.g., reactingFoam) that were validated for a variety of turbulent combustion applications under subsonic conditions [10]. As this solver lacks appropriate model construction for high-speed flows, a density-based compressible flow solver, namely rhoCentralFoam, was developed for high-speed flow modeling based on the central-upwind schemes of Kurganov and Tadmor [39,40]. Greenshields et al. [41] detailed this solver and validation against standard test cases. The solver was validated for compressible gas flow, including shocks and wave expansion. However, the rhoCentralFoam is limited to non-reacting flows. Therefore, in the current work, an in-house code, compressibleCentralReactingFoam, was developed by combining rhoCentralFoam with reactingFoam in the latest version of the OpenFOAM (i.e., v2306 [42]), to simulate high-speed compressible multicomponent reactive flows.

The challenges of scramjet simulation are due to multiscale (in terms of space and time) phenomena of non-steady, non-equilibrium turbulent combustion in the supersonic flow. Thus, the model needs to account for complex mixing phenomena, non-equilibrium transfer of turbulence energy, and the interaction between turbulence, and chemical kinetics. This paper demonstrates a physics-based methodology/computational tool, namely the compressibleCentralReactingFoam solver, to model the complex physical phenomena of reactive multicomponent gas flow at high Mach numbers. Moreover, the present work aims to quantitatively assess the accuracy and applicability of the newly developed solver in modeling supersonic combustion. A set of validation test cases are presented to demonstrate the capabilities of the solver ranging from simple one-dimensional shock-tube problems to complex three-dimensional scramjet combustor. In addition, three popular subgrid models - Smagorinsky model (SMG) [43,44], the Localized Dynamic k-equation Model (LDkEqn) [45,46], and Wall Adapted Local Eddy-Viscosity (WALE) [47] are evaluated for simulating supersonic turbulent combustion at conditions relevant to scramjet operation. Below, the short descriptions of these models are provided.

Subsection 2.1 presents the governing equations of supersonic multicomponent gas flow with chemical reactions and the basic approach to closure turbulent combustion equations. In Subsubsection 3.1.1, validations of 1D shock tube problems are solved. Subsubsection 3.1.2 presents simulations of 3D problems showing comparisons with experiments. The problem of turbulence-combustion interactions in the model scramjet combustor is presented in Subsubsection 3.1.3. Section 4 presents conclusions.

2. Materials and Methods

2.1. Governing Equations and Computational Solver of OpenFOAM Library

The governing equation of a reactive multicomponent mixture of gases consists of the continuity, the momentum, the energy balance, and the balance of mass fraction equations. A LES formulation is used to solve the equations for turbulent fluid flow. After Favre filtering, these equations can be written as [14,48]

$${\partial }_t\left(\overline{\rho }\right) +\mathbf{\nabla }·\left(\overline{\rho }\tilde{\mathit{u}}\right)=0,$$

(1)

$${\partial }_t\left(\overline{\rho }\tilde{\mathit{u}}\right)+\mathbf{\nabla }·\left(\overline{\rho }\tilde{\mathit{u}}\otimes \tilde{\mathit{u}}\right)=-\mathbf{\nabla }\overline{p}+\mathbf{\nabla }·\left(\overline{\mathit{S}}-\mathit{B}\right),$$

(2)

$${\partial }_t\left(\overline{\rho }{\tilde{\mathit{h}}}_s\right)+\mathbf{\nabla }·\left(\overline{\rho }\tilde{\mathit{u}}{\tilde{h}}_s\right)=\overline{\mathit{S}}·\mathbf{\nabla }\tilde{\mathit{u}}+{\partial }_t\overline{p}+\mathbf{\nabla }\overline{p}·\tilde{\mathit{u}}+\mathbf{\nabla }·\left(\overline{\mathit{h}}-{\mathit{b}}_{\mathit{h}}\right)+ {\overline{\dot{\omega }}}_T,$$

(3)

$${\partial }_t\left(\overline{\rho }{\tilde{Y}}_i\right)+\mathbf{\nabla }·\left(\overline{\rho }\tilde{\mathit{u}}{\tilde{Y}}_i\right)=\mathbf{\nabla }·\left({\overline{\mathit{j}}}_i-{\mathit{b}}_i\right)+{\overline{\dot{\omega }}}_i .$$

(4)

Here, the symbol “~” denotes Favre-filtered quantities, and “$-$“indicates filtered quantities; $\rho$ is the gas mixture density, $\mathit{u}$ is the velocity vector, $p$ is the pressure, $\mathit{S}$ is the viscous stress tensor, $h_s={{\displaystyle \sum }}_i{Y}_i{{\displaystyle \int }}_{T_0}^T{C}_{p,i}dT$ is the sensible enthalpy, $T$ is the temperature, $\mathit{h}$ is the heat flux vector, $Y_i$ is the species mass fraction, ${\dot{\omega }}_T=-{\displaystyle \sum }_{i=1}^M\left({\dot{\omega }}_i\Delta h_{f,i}^0\right)$ is the combustion heat release, $\Delta h_{f,i}^0$ is the formation enthalpy of $i^{th}$ species, ${\dot{\omega }}_i$ is the species reaction rate, ${\mathit{j}}_i$ is the species mass flux. Subgrid pressure fluctuations and dissipation terms will be neglected.

The subgrid flow physics is concealed in the subgrid stress tensor $\mathit{B}=\overline{\rho }\left(\widetilde{\mathit{u}·\mathit{u}}-\tilde{\mathit{u}}·\tilde{\mathit{u}}\right)$and flux vectors ${\mathit{b}}_i=\overline{\rho }\left(\widetilde{\mathit{u}{Y}_i}-\tilde{\mathit{u}}\widetilde{Y_i}\right)$, ${\mathit{b}}_{\mathit{h}}=\overline{\rho }\left(\widetilde{\mathit{u}{h}_s}-\tilde{\mathit{u}}{\tilde{h}}_s\right)$, which results from filtering the nonlinear convective terms. Following [2], one may postulate that the gas mixture behaves as linear viscous fluid with Fourier heat conditions and Fickian species diffusion, such as ${\overline{\mathit{j}}}_i\approx D_i\mathbf{\nabla }\widetilde{Y_i}$, $D_i=\mu /S{c}_i$ is the species diffusivities, $\overline{p}\approx \overline{\rho }{\mathcal{R}}_u\tilde{T}{{\displaystyle \sum }}_{i=1}^M\left(Y_i/M_i\right)$, ${\mathcal{R}}_u$ is the universal gas constant, $M_i$ is the molar mass of $i^{th}$ species, $\overline{\mathit{S}}\approx 2\mu {\tilde{\mathit{D}}}_D$ is the viscous stress tensor, $\mu$ is the dynamic viscosity which is modeled using Sutherland’s law $\mu =A_s\sqrt{T}/\left(1+T_s/T\right)$, where $A_s, T_s$ are Sutherland’s coefficient and temperature[30], ${\tilde{\mathit{D}}}_D=\left(\mathbf{\nabla }\tilde{\mathit{u}}+\mathbf{\nabla }{\tilde{\mathit{u}}}^T\right)/2-\left(\mathbf{\nabla }·\tilde{\mathit{u}}\right)\mathit{I}$/3 is the deviatoric part of the rate strain tensor, $\overline{\mathit{h}}\approx \kappa \mathbf{\nabla }\tilde{T}$, $\kappa = \mu /Pr$ is the thermal diffusivity, $\mathit{S}{c}_i$ and $Pr$ are the Schmidt and Prandtl numbers, respectively.

The LES equation (1)-(4) is closed by providing submodels for $\mathit{B}$, ${\mathit{b}}_i$, ${\mathit{b}}_h$. The most common subgrid models used in this paper for high-speed reactive flows are class of Smagorinsky models [43,44] of the form $\mathit{B}=-2{\mu }_k{\tilde{\mathit{D}}}_D+\left(2/3\right)\overline{\rho }k\mathit{I}$, where ${\mu }_k$ is the subgrid viscosity and $k$ is the subgrid kinetic energy, ${\mathit{b}}_i=-\left({\mu }_k/S{c}_t\right)\nabla {\tilde{Y}}_i$, and ${\mathit{b}}_h=-\left({\mu }_k/P{r}_{\mathrm{t}}\right) \nabla {\tilde{h}}_s$, where $S{c}_t\approx 0.7$ and $P{r}_{\mathrm{t}}\approx 0.8$ are the turbulent Schmidt and Prandtl numbers, respectively. To closure $k$ and ${\mu }_k$ in the Smagorinsky model (SMG) [44], they are formulated as ${\mu }_k=c_D\overline{\rho }{·}^2\left|\tilde{\mathit{D}}\right|$ and $k=c_l{·}^2{\left|\tilde{\mathit{D}}\right|}^2$, where $c_D$ and $c_l$ are model coefficients.

The second subgrid model used in this work is the Local Dynamic k-equation (LDkEqn) [45,46]. The LDkEqn turbulent model is an advanced subgrid-scale (SGS) turbulence model used in Large Eddy Simulation to improve the prediction of turbulent flows. Instead of using fixed model coefficient constants (as in the SMG), the LDkEqn model calculates them dynamically, adjusting the model coefficients to match the resolved flow characteristics. This adaptability allows the model to capture variations in turbulence intensity and anisotropy better. The LDkEqn model exhibits a scale-adaptive behavior, which means it can adapt to different turbulence scales in the flow. This adaptability is particularly beneficial in LES, where a wide range of turbulent structures exists; also, the LDkEqn model is designed to provide better predictions near walls. For the LDkEqn, a modeled transport equation ${\partial }_t\left(\overline{\rho }k\right)+\nabla ·\left(\overline{\rho }k\tilde{\mathit{u}}\right)=-\mathit{B}·\tilde{D}+\mathbf{\nabla }·\left({\mu }_k\mathbf{\nabla }k\right)-c_{\epsilon }\overline{\rho }{k}^{3/2}/·$ is solved for k and subgrid viscosity is modeled as ${\mu }_k=c_k\overline{\rho }· k^{1/2}$. The model coefficients $c_k$ and $c_{\epsilon }$ are evaluated dynamically using scale similarity.

The third subgrid model used in this work is the Wall Adapted Local Eddy-Viscosity (WALE) subgrid model [47]. The WALE model is designed to improve the accuracy of LES by providing a more accurate representation of the subgrid-scale turbulence in the near-wall region, where traditional SGS models may struggle to capture complex flow phenomena. The WALE model is specifically designed to address the near-wall region of turbulent flows. In this region, small-scale turbulence structures are not well-resolved on the computational grid, and traditional SGS models may introduce excessive dissipation or fail to capture these structures accurately. The WALE model computes a subgrid-scale turbulent viscosity (${\nu }_t$) that varies spatially and temporally based on the local flow characteristics. The specific coefficients used in the WALE model are $C_k=0.094$; $C_e=1.048$; $C_w=0.325$.

Finally, to completely close the LES equations (1)-(4), the filtered reaction rates. ${\overline{\dot{\omega }}}_i$, incorporating the turbulence chemistry interactions, must be provided. The reaction rate provides the effect of combustion chemistry on the flow and how the turbulence affects the combustion chemistry. The LES Finite Rate Chemistry model (LES-FRC), specifically the LES Partially Stirred Reactor (LES-PaSR) model [48], is to be used in this work.

This study adopts the PaSR subgrid combustion model for the sequential and simultaneous processes of micro-mixing and chemical reactions. It is based on the conjecture that any turbulent flow can be divided into the fine structure $\left(\ast \right)$ and surrounding structure $\left(0\right)$, interacting through the balance equations [49,50]. Most dissipation and mixing, as well as chemical reactions, occur in the fine structure, where reactants are mixed at the molecular level. Therefore, each LES cell can be viewed as a partially stirred reactor containing a fine structure, ideally regarded as a Perfectly Stirred Reactor (PSR), exchanging mass and energy with surroundings. The filtered reaction rates for PSR are defined as $\overline{{\dot{\omega }}_i^{*}}=M_i{\displaystyle \sum }_{\mathit{j}=1}^M{P}_{i\mathit{j}}\overline{\widehat{\dot{{\omega }_{\mathit{j}}}}}$, where $\widehat{\dot{{\omega }_{\mathit{j}}}}=k_{f\mathit{j}}{{\displaystyle \prod }}_{k=1}^M{\left(X_k\right)}^{P_{k\mathit{j}}^{\prime }}-k_{r\mathit{j}}{{\displaystyle \prod }}_{k=1}^M{\left(X_k\right)}^{P_{k\mathit{j}}^{″}}$, $P_{i\mathit{j}}=P_{i\mathit{j}}^{″}-P_{i\mathit{j}}^{\prime }$ are the stochiometric coefficients, $\left(X_k\right)=\left(\rho Y_k/M_k\right)$ is molar concentration, $k_{fj}$ and $k_{rj}$ being the Arrhenius forward and reverse rates of the $jth$ reaction step, $k_{fj}=A_{fj}{T}^{{\beta }_j}\exp \left(-T_{a,j}/T\right)$, ${\beta }_j$ and $T_{a,j}$ are the Arrhenius parameters. The resulting rates in the PaSR model are [51], ${\overline{\dot{\omega }}}_i={\gamma }^{*}\overline{{\dot{\omega }}_i^{*}}+\left(1-{\gamma }^{*}\right)\overline{{\dot{\omega }}_i^0}$, which can be simplified as ${\overline{\dot{\omega }}}_i\approx {\gamma }^{*}\overline{{\dot{\omega }}_i^{*}}$, since $\overline{{\dot{\omega }}_i^0}$ are most negligible outside of the fine structure regions, where the fraction of reactive zone (fine structure) ${\gamma }^{*}$ is calculated based on geometrical analysis [2] ${\gamma }^{*}\approx {\tau }_c/\left({\tau }^{*}+{\tau }_c\right)$, ${\tau }_c$ is a global representative of chemical time step and ${\tau }^{*}=\sqrt{{\mu }_{eff}/\left(\overline{\rho }\tilde{\epsilon }\right)}$ is the turbulent mixing time that is implemented in the OpenFOAM, ${\mu }_{eff}$ is the sum of laminar and turbulent dynamic viscosity, $\tilde{\epsilon }$ is the turbulent dissipation rate.

3. Results

3.1. Solver Validation

The capability of compressibleCentralReactingFoam to solve compressible flows with multicomponent reacting gas mixtures will be demonstrated for various problems. These include 1D shock tube problems, 2D and 3D geometries involving underexpanded jet formation, and 2D turbulent reacting gas flow in a scramjet combustor.

3.1.1. Shock Tube Problem.

Shock tube is a well-known benchmark problem (also known as Sod shock tube problem [52]) for the validation of compressible codes. It is an attractive test case for model validation due to the simple one-dimensional geometry of computational area and the possibility to simulate a wide range of flow physics, including shock and rarefaction waves, contact discontinuities, and complex wave interactions. The available analytical and computational solutions in the literature allow for the assessment of the code's accuracy and reliability. The shock tube problem for two cases have been solved: (1) a shock tube with an inert multicomponent gas mixture, and (2) a shock tube with a chemically reactive multicomponent gas mixture. In these validations, the initial discontinuity of the gas in the shock tube is decayed into different gas configurations such as shock waves, contact discontinuous, and rarefaction waves spreading into different directions (See Figure 1).

• Inert Multicomponent Shock Tube

The first validation case is a multicomponent inert shock tube simulation without the chemical source term in order to avoid shock-induced autoignition of the gas mixture. A homogeneous premixed gas mixture of 20%H₂/10%O₂/70%Ar (in mole%) is used in the simulation to verify the numerical accuracy of the model implementations for thermal and mass transport of multicomponent mixture. The shock tube is 0.10 m in length, and it is discretized with 400 uniform elements. The time step is equal to $\Delta t={10}^{-7} s$. In the middle of the tube is a fixed diaphragm that separates the gas between the left side ($T_l = 400 K$ and $P_l = 8000 Pa$) and the right side ($T_r = 1200 K$ and $P_r = 80,000 Pa$). At time $t=0,$ the diaphragm is eliminated, and the shock wave starts to propagate from right to left and expansion wave moves in the opposite direction. The initial configuration of the gas mixture generates the propagation of shock and discontinuous waves to the left (low-pressure) and a rarefaction (expansion) wave to the right side (high-pressure) of the tube. Figure 2 shows the current simulation results (solid line) compared with the simulation results of Huang et al. [53]. The modeling results have good agreement. It is noteworthy that the compressibleCentralReactingFoam solver's accuracy is 2nd order for both time and spatial convective terms.

• Reactive Multicomponent Shock Tube

A reactive multicomponent gas mixture in the shock tube is also modeled using the conditions of Martínez-Ferrer et al. [34] for comparison. The molar composition of the gas mixture is 20% H₂/10% O₂/70% Ar. The length of the tube is equal to $L=0.12 m$ with the initial conditions corresponding to the left and right sides as follows: $\left({\rho }_l,u_l,p_l\right)=\left(0.072, 0, 7173\right)$ and $\left({\rho }_r,u_r,p_r\right)=\left(0.18075, -487.34, 35 594\right)$, where $\rho , u , \mathrm{and} p$ are measured in ($kg/m^3$), ($m/s$), and ($Pa$), respectively. The length of the tube is discretized by 400 uniform elements with 0.01 $\mu s$ timestep. The simulation used the reaction mechanism of Conaire et al. [54] for hydrogen which consists of nine (9) chemical species (H₂, O₂, H, O, OH, HO₂, H₂O₂, H₂O, $Ar$), and 18 elementary reactions. A wall boundary condition is implemented on the left side of the tube, and on the right side, a non-reflected boundary condition is implemented. The initial configuration of non-stable discontinuity is decayed into a shock propagated to the tube's left side. The solutions of spreading discontinuities in the shock tube for velocity, mass fraction of $H$, and temperature: $U_x$, $Y_H$, and $T$ are shown in Figure 3 for 230 $\mu s$ simulation time. Comparisons of current solution (solid lines) in Figure 3 with computational data of Martínez-Ferrer et al. [34] (symbols) show good agreement. The visible discrepancies are seen in the vicinities of local maximums around x = 0.06 m at 230 $\mu s$ that is explained of using by Martínez-Ferrer et al. [34] a seventh-order accurate Weighted Essentially Non-Oscillatory (WENO) scheme to discretize the non-linear advective terms.

The decay process of the initial multicomponent gas mixture, and how the shock wave propagate in the tube was investigated in detail by comparing the profiles of inert and reactive mixture of H₂/O₂/Ar. Figure 4 shows the simulation results for the temperature profiles of inert and reactive mixture at various simulation times.

The shock starts to propagate from its initial position in the middle of the tube from right to left toward the wall (i.e., x = 0) with a velocity $v_{rl}$ = 8100 m/s and reaches the wall around 75 $\mu s$. After reflection from the wall, the shock moves in the opposite direction from left to right. The reflected shock moves with a velocity of 4500 m/s that is almost half of $v_{rl}$ as the gas flows right to left with a velocity of 490 m/s. Due to the interaction of the shock with the wall, the temperature of the reflected shock is increased from 750 K to 1200 K for both inert and reactive cases. As seen in Figure 4b, the shock-induced autoignition starts to occur (i.e., a slight increase in temperature at x = 0) at 120 $\mu s$ for the reactive case as the simulation time is longer than the chemical induction time to generate combustion radicals. Due to the combustion heat release, the gas temperature increases to 2100 K as the simulation time is increased to 150 $\mu s$. The gas temperature approaches equilibrium state as the simulation time is further increased to 220 $\mu s$.

3.1.2. Simulation of Ladenburg Jet Problem.

Simulation of a well-known experiment of underexpanded jet problem by Ladenburg et al. [55] is conducted to verify code accuracy for the shock formation of supersonic jets. Two-dimensional axisymmetric simulation of the problem using rhoCentralFoam was reported previously [41,56]. In the current work, 3D simulation of the Ladenburg experiment is performed with the compressibleCentralReactingFoam solver. As the maximum temperature of the experiment is no more than 300 K, the chemical reactions will be at frozen state and hence, the gas mixture is considered be inert. The computational 3D mesh (see Figure 5) has 5 mm inlet radius, 10 mm free surface radius, and 30 mm length. The following inlet conditions for pressure, axial velocity and temperature are used to ensure a sonic condition with Mach number of one (1) at the inlet surface: $2.72x{10}^5 Pa$, $316.6 m/s$ and $247.1 K$. The conditions at the free stream surface are: $1.01x{10}^5 Pa$, $0 m/s$ and $297 K$. Non-reflected boundary conditions are used at the outlet. The simulations are performed using a mesh with total hexahedral elements equal to 1.72 million. The sketch showing wave structures and the simulation flow of Mach number contours is shown in Figure 6a. The pressure at the inlet of the computational region is the pressure of the reservoir, which is greater than the ambient pressure; thus, the flow is defined as underexpanded. When air accelerates from high inlet pressure to the atmosphere, the gas flow declines due to its expansion and acceleration. The gas jet expands to the atmospheric pressure through an expansion fan (Figure 6a). It is seen in Figure 6b in the free jet boundary above which the gas flow is at rest (i.e., M=0). In Figure 6a, numbers 1-5 in the circles indicate different gas structures, such as expansion and shock waves. The expansion waves (denoted by 1) are reflected from the free boundary as compression waves. The accumulated compression waves at a point then form an oblique shock wave, shown by the solid line (denoted by 2). The formed oblique shock waves are again reflected off the center line of the flow. This reflection causes the formation of a Mach reflection before the intersection point with the centerline, producing a triple shock point where the compression waves, oblique shock, and centerline coincide. The incoming oblique shock wave (denoted by 2) is the first Mach reflection. The second shock (denoted by 3) is a Mach reflection, perpendicular to the centerline termed a normal (Mach disk) shock. The third Mach reflection (denoted by 4) is another oblique shock wave that is reflected off the constant pressure-free air-jet boundary as an expansion fan (denoted by 5). The escalation of the pressure and temperature occurs as the flow passes through the normal shock or Mach disk. So, the Mach disk is the high-pressure region in an underexpanded air jet that is formed through shock waves and expansion waves created by the vast difference between the inlet pressure and the ambient pressure.

Figure 7 shows the comparison of the current simulation with experimental data of Ladenburg et al. [55] for density iso-lines. The simulation results have a good agreement with the experimental data, especially the location and height of the Mach disk. The experimental results show that a weak shock is produced due to the expansion of air from the nozzle orifice (close to density $3.8 kg/m^3$) and then extends toward the nozzle axis. Thus, it creates a Mach disk feature and a triple point (i.e., the intersection of Mach disk, barrel shock and reflected shock). The experimental data in Figure 7 (i.e., lower panel) shows that the position of the triple point is at 1.7 mm radial distance and 13.3 mm axial distance. The computed radial and axial distance of the triple point are 1.65 mm and 13.5 mm, respectively. Thus, the computational results yield less than 3% error compared to the experimental data. In summary, the computational results are in good agreement with the experimental data as well as the physics of the shock formation, expansion wave behavior, and their interactions with jet boundary and slip lines.

3.1.3. Simulation of DLR Scramjet Combustor

The combustion experiment of the Waidmann et al. [57] is selected to validate the compressibleCentralReactingFoam solver to model the turbulence combustion interactions in a scramjet combustor. The experimental data were obtained for hydrogen combustion using the Scramjet test facility at the German Aerospace Center (DLR) at Mach number two. The experimental system was modeled using 2D and 3D computational studies by Oevermann [58] and Haung et al. [35], respectively. Oevermann [58] used RANS (i.e., k-ε turbulence model) with laminar flamelet whereas Huang et al. [35] used LES with PaSR sub-grid scale combustion model to simulate the Scramjet. Both computational studies produced reasonable agreement with the experimental data. In the current work, 2D LES simulation is conducted with PaSR sub-grid scale combustion model to validate compressibleCentralReactingFoam solver. Although the two-dimensional assumption may introduce some uncertainty when directly comparing it to the experiment, the 2D computational study will help to verify the accuracy of the model implementation in compressibleCentralReactingFoam qualitatively.

Figure 8 shows the geometry of DLR test rig. The combustor has a height of $50 \mathrm{mm}$ and width of $45 \mathrm{mm}$ at air inflow. The divergence angle of the upper wall is $3^0$ to compensate for the expansion of the boundary layer. The following boundary conditions for air are established to ensure the Mach number ${\mathrm{M}}_{\mathrm{air}}\approx 2$: ${\mathrm{u}}_{\mathrm{air}}=730 \mathrm{m}/\mathrm{s}$, ${\mathrm{T}}_{\mathrm{air}}=340 \mathrm{K}$, ${\mathrm{p}}_{\mathrm{air}}={10}^5\mathrm{Pa}$ with the mass rate of airflow ${\dot{\mathrm{m}}}_{\mathrm{air}}\approx 1.0 \mathrm{kg}/\mathrm{s}$. The injector of fuel (${\mathrm{H}}_2$) is a slot in the center of the wedge-shaped strut with a height $\Delta {\mathrm{y}}_{\mathrm{H}2}=0.26 \mathrm{mm}$ that gives the surface area of injection $\Delta {\mathrm{s}}_{\mathrm{H}2}=11.7{ \mathrm{mm}}^2$. The following boundary conditions are used to establish fuel injector Mach number ${\mathrm{M}}_{\mathrm{H}2}\approx 1$: ${\mathrm{u}}_{\mathrm{H}2}=1200 \mathrm{m}/\mathrm{s}$, ${\mathrm{T}}_{\mathrm{H}2}=250 \mathrm{K}$, ${\mathrm{p}}_{\mathrm{H}2}={10}^5\mathrm{Pa}$ with the hydrogen mass flow ${\dot{\mathrm{m}}}_{\mathrm{H}2}\approx 1.26 \mathrm{g}/\mathrm{s}$. The no-slip boundary conditions are implemented on the solid walls. The computational region is discretized on 100,000 hexahedrons (See Figure 9). A detailed kinetic mechanism for hydrogen combustion with 9 species and 18 elementary reactions [59] is used.

• Simulation of Cold Flow DLR Combustor

Initially, the cold flow is simulated. Figure 10 shows a Schlieren photograph of the channel flow and computational data of density gradient contours corresponding to the experiment. Comparisons of these two pictures show that the current compressible solver is able to accurately capture the aerodynamics features such as shock waves, expansion waves, interaction, and reflection of these structures from the walls and jet-shear layer. The computational results shown in Figure 10b also identify various flow characteristics with labels 1 to 7. Label (1) shows the leading shocks formed on the tip of the edge (see Figure 10). Label (2) refers to the expansion waves due to the supersonic flow past the corners of the edge. Label (3) is the reflected shock of (1) from the upper wall whereas label (4) refers to the shocks formed due to collision of gas flow passing through the corner of the edge. Label (5) is an oblique shock formed from merging of the reflected shock and shocks in (4). Lebel (6) is a reflected from the upper wall oblique shock whereas label (7) is the jet shear-layer formed from the Kelvin-Helmholtz instability. The oblique shocks that are reflected from the upper and lower walls interact with the jet shear layer to get reflected again without penetration as shown Figure 10b.

Figure 11a shows the normalized contours of density gradient $|\mathbf{\nabla }\rho /{\rho }_{\max }$|. Figure 11b shows the Mach numbers of the cold flow. Shocks are generated from the tip of the wedge and are further reflected from the lower and upper walls of the combustor. Also, the expansion wave on the upper wall at the starting divergence is seen. Gas flow passing the corners of the wedge are subjected to expansion by forming two expansion waves. The jet vortex structure arises from the Kelvin-Helmholtz instability of the jet-shear layer [14]. It is noteworthy that the simulation with the $k-\epsilon$ turbulent model did not generate shear-layer instabilities as reported by Oevermann [58]. This finding confirms that $RANS$ models are limited when predicting such phenomena and $LES$ models are more appropriate in scramjet combustor simulations.

Figure 12 compares the cold flow simulation results for pressure distributions with the experimental measurements of Waidmann et al. [57]. Figure 12a show the distribution of pressure on the lower (solid line) and upper (dashed line) walls along with experimental data obtained on the lower wall (symbols). Figure 12b compares the pressure along the center line with the experimental data. The model predictions are able to capture the pressure profiles observed experimentally.

• Simulation of Reacting Flow DLR Combustor

Simulations were performed for reactive flow with detailed hydrogen kinetic mechanism for experimental conditions of Waidmann et al. [57]. Combustion was initiated with the help of an igniter during the test as no autoignition was observed experimentally. Our initial simulation also showed that combustion did not occur without an external ignition source. Hence, ignition zones at the two corners of the strut with temperature $T_{ign}=1500K$ are implemented. Figure 13 compares Schlieren's photograph of the reactive flow with the computational data for the density contours obtained in the current work for the experimental condition of Waidmann et al. [57]. The comparison shows similar gas dynamic features such as shock waves, expansion waves, interaction, and reflection of the flow structures from the walls and jet-shear layer. Significant gas dynamic differences can be seen between cold flow (Figure 10) and hot flow (Figure 13). For example, expansion waves (i.e., labeled 2 in Figure 10) at the corners of the wedge in cold flow are replaced with the shocks generated due to the combustion in Figure 13. It should be noted that a considerable subsonic region is developed behind the hydrogen injection for reactive flow.

• Influence of Turbulent Models on Scramjet Combustion

A discussion on the influence of different turbulent models on supersonic combustion is presented in this section. Scramjet simulations for the DLR Scramjet experiment [57] are performed with the following turbulent models described in Subsection 2.1: SMG, LDkEqn, and WALE. Figure 14 compares the experimental data for temperature at three cross sections along the combustor with the simulation results obtained with three turbulent models. The figure also shows the location of the temperature measurements in the DLR Scramjet combustor. The simulation results close to the strut (i.e., x = 62 mm) show that LDkEqn and WALE models have similar profiles. Further downstream (i.e., x = 216 mm), all three models have similar profiles. Overall, WALE model has slightly better agreement with the experimental data.

Figure 15 compares the current WALE model simulation results with other DLR Scramjet simulation data reported in the literature [12,30,58,60]. Modeling result of Berglund and Fureby, [60] closely follow the experimental profile except at x = 62 mm. Overall, the present model has better agreements with the experimental data at all locations, but it does show relatively higher diffusivity at downstream locations (i.e., x = 108 mm and 216 mm).

• Analysis of Turbulence-Combustion Interaction in Scramjet Combustor

Based on the analysis of different turbulent models, the WALE model is selected as the most suitable for Scramjet simulation. The effect of turbulence-combustion coupling on gas dynamics and flame structure is discussed in this section using the DLR Scramjet simulation results obtained with WALE model. Figure 16 shows the Mach number (Ma), rate of heat release per unit volume (Qdot), instantaneous temperature (T), and normalized contour lines of densities’ gradient ($|\nabla \rho /{\rho }_{\max }$|) obtained from the simulation. The shock structure, shown in Figure 16a,d, is characterized by strong interaction with the shear layers characterized by large density gradients. Behind the strut, a central recirculation bubble is generated with a subsonic region. The formation of the recirculation zone is very important for flame holding and, hence, stabilizing the flame in the scramjet combustor. The rate of heat release contours in Figure 16b show that combustion is intensive within the thin shear layers. The combustion is initiated within the shear layers separated by the recirculation zone close to the edges of the strut. This observation is supported by the temperature data presented in Figure 14a.

The direction of diffusive spreading of fuel and oxidizer helps to analyze turbulence-chemistry interactions. It can be estimated by the interaction between the gradients of fuel ($\mathbf{\nabla }{Y}_F$) and oxygen ($\mathbf{\nabla }{Y}_O$) mass fractions. Takeno and co-workers [61] suggested the use of a flame index (known as the Takeno Flame Index - TFI) based on scalar productions of the gradients of the fuel and oxidizer mass fraction, which help to distinguish between the premixed and non-premixed zones [61]. The TFI is defined as $G_{FO}=\left(\mathbf{\nabla }{Y}_F·\mathbf{\nabla }{Y}_O\right)/\left(\left|\mathbf{\nabla }{Y}_F\right| \left|\mathbf{\nabla }{Y}_O\right|\right)$. Positive and negative values of TFI indicate non-premixed or premixed modes, respectively [61]. Figure 17 shows the TFI calculated for the DRL Scramjet conditions. Pockets of non-premixed zones shown in red are surrounded by large premixed zone shown in dark blue.

Figure 18 shows the species mass fraction contours for $OH$, $H$, $H{O}_2$ and $H_{2}O$. It is observed [9] that the OH field is far from uniform in the mixture due to turbulence, which would have been expected from fast chemistry arguments. High OH regions are found where the mixture fraction iso-surfaces are highly convoluted, and low values are in the areas where the mixture fraction iso-surfaces are stretched and not wrinkled. The mass fraction of OH is relatively small at the beginning, while it is most prominent in the downstream regions, especially in the transition zone. The OH radical represents the existence of a flame front which correlates with heat release rate. Field of H2O correlates with the field of instantaneous temperature while the field of HO2 correlates with OH and Qdot on the shear layer.

Figure 19 shows mixture fraction ($Z$) and scalar dissipation rate ($\chi )$. The mixture fraction measures the local fuel/oxidizer ratio, and calculated as $Z=\left(\phi ·Y_F/Y_{F0}-Y_O/Y_{O0}+1\right)/\left(\phi +1\right)$ where $\phi$ is equivalence ratio. The scalar dissipation rate is calculated as $\chi =2D\left(\mathbf{\nabla }Z·\mathbf{\nabla }Z\right)$, where $D$ is the mixture diffusion coefficient. The scalar dissipation rate is essentially the rate of mixing between fuel and oxidizer [9] and represents the local mixing rate at the molecular level [62]. The scalar dissipation rate is critical for modeling the effect of turbulence on reaction rates as it has significant influence on non-premixed combustion. It often provides the connection between the molecular mixing and the combustion. In turbulent flows, the scalar dissipation is seen as a scalar energy dissipation. Its role is to destroy (dissipate) scalar variance (scalar energy) analogous to the dissipation of the turbulence. Unlike the kinetic energy dissipation, most of the scalar dissipation occurs at the finest scales. The scalar dissipation has higher values near the strut and decreases further downstream as shown in Figure 19b.

A time-scale analysis based on the estimation of the local Damköhler Number (Da) is also calculated to examine the reactive zones in more detail. The Da number, quantified as the ratio between the turbulent mixing time scale ${\tau }^{*}$ and the chemical time scale ${\tau }_c$: $Da={\tau }^{*}/{\tau }_c$ [63,64]. The Borghi diagram [63] shown in Figure 20a helps to recognize different regimes of turbulent combustion. Figure 20b,c show the Damköhler number (Da) on a logarithmic scale and the fraction of fine structure (${\gamma }^{*}$), respectively. The bright white lines indicate the locations where $Da=1$. The Damköhler number is relatively high ($Da\gg 1$) in the core of the jet flow. This indicates that the time scale of chemical reaction is smaller than the turbulent mixing time scale which means combustion regime with fast chemistry. Thus, the fraction of reactive cells (${\gamma }^{*}$) (see Section 0 for details) is mostly small in this zone. In contrast, in the shear layer where the rate of the heat release (see Figure 16b) and the scalar dissipation rate (see Figure 19b) have the highest values, the flame is controlled by finite rate chemistry with $Da\ll 1$ and ${\gamma }^{*}\approx 1.$In the core flow, it is also seen small pockets where $Da\ll 1$ indicating the combustion zone is in the “island formation” regime as shown in Figure 20a. These zones are characterized by the combustion regime with a finite rate chemistry with a rather high fraction of reactive cells (${\gamma }^{*}~1$) and are called partially-stirred reactor.

4. Conclusions

A computational validation of a new compressibleCentralReactingFoam code for modeling combustion in supersonic flows is presented. The validation cases include various complexities of high-speed flows with shocks, wave expansion, and turbulence-combustion interactions. Comparisons of the simulation results with the literature data demonstrated the numerical fidelity of the solver for reactive multicomponent gas mixtures at supersonic conditions. Moreover, detailed analysis on gas dynamics and turbulence-combustion interaction under supersonic conditions is presented to gain a deeper understanding of the complex physical-chemical phenomena involved. A comparative analysis on the influence of three turbulent subgrid models on scramjet combustion is presented. The comparison of simulation results with the experimental data shows that WALE subgrid scale is most suitable for Scramjet modeling. Moreover, the DLR Scramjet simulation results show that the flame is mainly stabilized in the combustor due to the forming of a subsonic bubble behind the strut. The supersonic flame structure is investigated through a comprehensive analysis of chemical species formation, scalar dissipation rate, flame index, heat release rate, fine structure, and Damköhler number. A time-scale analysis based on the local Damköhler number reveals different regimes of turbulence combustion in the Scramjet. In the core of the jet flow, the Da number is relatively high ($Da\gg 1$) indicating that the chemical time scale is smaller than the turbulent mixing time scale. This means that the combustion regime is dominated by fast chemistry. The results also show that the heat release rate and the scalar dissipation rate have the highest values in the shear layer, and the flame is stabilized due to finite rate chemistry with $Da\ll 1$ and ${\gamma }^{*}\approx 1$. Thus, the new solver presented in this paper can be used for high-speed simulation of multicomponent reactive gas mixtures in supersonic combustors.