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

: An approach based on the OpenFOAM library has been developed to solve a high-speed, multicomponent mixture of a reacting, compressible flow. This work presents comprehensive validation of the newly developed solver, called compressibleCentralReactingFoam , with different supersonic flows, including shocks, expansion waves, and turbulence–combustion interaction. The comparisons of the simulation results with experimental and computational data confirm the fidelity of this solver for problems involving multicomponent high-speed reactive flows. The gas dynamics of turbulence– chemistry interaction are modeled using a partially stirred reactor formulation and provide promising results to better understand the complex physics involved in supersonic combustors. A time-scale analysis based on local Damköhler numbers reveals different regimes of turbulent combustion. In the core of the jet flow, the Damköhler number is relatively high, indicating that the reaction time scale is smaller than the turbulent mixing time scale. This means that the combustion is controlled by turbulent mixing. In the shear layer, where the heat release rate and the scalar dissipation rate have the highest value, the flame is stabilized due to finite rate chemistry with small Damköhler numbers and a limited fraction of fine structure. This solver allows three-dimensional gas dynamic simulation of high-speed multicomponent reactive flows relevant to practical combustion applications.


Introduction
The supersonic airbreathing or supersonic combustion ramjet (scramjet) engine is the most promising technology for significant reduction of flight times for long-distance civilian flights and the first stages of space systems for putting payloads 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 that 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 air 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, on the order of milliseconds [4], which is insufficient to provide a comprehensive analysis for designing efficient scramjet engines.A numerical approach remains the primary analysis tool to cover on code development.Among the open-source codes, the OpenFOAM library is a wellestablished 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 highspeed 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 under conditions relevant to scramjet operation.Below, short descriptions of these models are provided.Section 2.1 presents the governing equations of supersonic multicomponent gas flow with chemical reactions.Subgrid flow equations and closure of the combustion model are discussed in Sections 2.2 and 2.3, respectively.Section 2.4 briefly discusses the numerical method and algorithm of the compressibleCentralReactingFoam solver.In Section 3.1.1,validations of 1D shock tube problems are solved.Section 3.1.2presents simulations of 3D problems showing comparisons with experiments.The validation of turbulence-combustion interactions in the model scramjet combustor is presented in Section 3.1.3.The influence of turbulence models on scramjet combustion and analysis of turbulent-combustion interaction are presented in Sections 3.2 and 3.3, respectively.Section 4 presents conclusions.

Governing Equations
The governing equations of a reactive multicomponent mixture of gases consist 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] Dynamics 2024, 4 138 Here, the symbol "~" denotes Favre-filtered quantities, and "−" indicates filtered quantities; ρ is the gas mixture density, u is the velocity vector, p is the pressure, S is the viscous stress tensor, E is the total energy, h is the heat flux vector, .ω T is combustion heat release, Y i is the specific mass fraction, j i is the diffusive flux, and .ω i is the species reaction rate of i th species.M is a number of species, and due to the conservation ∑ i Y i = 1, the last species fraction is defined from this equation.The terms of subgrid flow physics in Equations ( 1)-( 4) are denoted by SGS upper indexes.The pressure of the gas mixture is equal to where R u is the universal gas constant, T is the temperature, and W i is the molar mass of i th species.The viscous stress tensor is defined as where D D is the deviatoric part of the rate strain tensor and µ is the dynamic viscosity, which is modeled using Sutherland's law where A s , T s are Sutherland's coefficient and temperature [30].The deviatoric part of the rate strain tensor is defined as The total energy is defined as where C p,i is the heat capacity at constant pressure, which is defined from JANAF polynomials [49].Following [2], one may postulate that the gas mixture behaves as a linear viscous fluid with Fourier heat conditions and Fickian species for diffusive flux, such as and where κ = µ/Pr is the thermal diffusivity, Pr is the Prandtl number, D i = µ/Sc i is the species diffusivities, and Sc i is the Schmidt number.The combustion heat release is defined as .
where ∆h 0 f ,i is the formation enthalpy of i th species.The reaction rate is defined as .
where ω 0 ij is the reaction rate of elementary reaction, and N > 1 when detailed chemical mechanism is used.

Subgrid Flow Equations
Subgrid pressure fluctuations and dissipation terms will be neglected.The subgrid flow physics is concealed in the subgrid stress tensor and flux vectors b which results from filtering the nonlinear convective terms.
The LES Equations ( 1)-( 4) are closed by providing submodels for The most common subgrid models used for high-speed reactive flows are a class of Smagorinsky models [43,44] of the form where µ k is the subgrid viscosity and k is the subgrid kinetic energy, where Sc t ≈ 0.7 and Pr t ≈ 0.8 are the turbulent Schmidt and Prandtl numbers, respectively.To close k and µ k in the Smagorinsky model (SMG) [44], they are formulated as and where ∆ is the filter size, and c D and c l are model coefficients.
The second subgrid model used in this work is the Local Dynamic k-equation (LD-kEqn) [45,46].The LDkEqn turbulence 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 by 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 exist; also, the LDkEqn model is designed to provide better predictions near walls.For the LDkEqn, a modeled transport equation is solved for k, and subgrid viscosity is modeled as The model coefficients c k and c ε are evaluated dynamically using scale similarity.
The third subgrid model used in this work is the Wall Adapted Local Eddy-Viscosity (WALE) 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 (ν 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.

Closure of Combustion Model
Finally, to completely close the LES Equations ( 1)-( 4), the filtered reaction rate .ω 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.A finite rate chemistry model, namely the LES Partially Stirred Reactor (LES-PaSR) model [48], is used in this work.
The solver used 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 ( * ) and surrounding structure (0), interacting through the balance equations [50,51].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 that is ideally regarded as a Perfectly Stirred Reactor (PSR), exchanging mass and energy with surroundings.The filtered reaction rates for PSR are defined as . .
the stochiometric coefficients are defined as The molar concentration is defined as where k f j and k rj being the Arrhenius forward and reverse rates of the jth reaction step, β j and T a,j are the Arrhenius parameters.The resulting reaction rate in the PaSR model is defined [52] as .
which can be simplified as .
. ω 0 i are most negligible outside of the fine structure regions, where the fraction of reactive zone (fine structure) γ * is calculated based on geometrical analysis [2] Dynamics 2024, 4 where τ c is a global representative of chemical time step and is the turbulent mixing time that is implemented in the OpenFOAM, µ e f f is the sum of laminar and turbulent dynamic viscosity, and ∼ ε is the turbulent dissipation rate.

Numerical Methods
In this section, a brief overview of the numerical methods implemented in the com-pressibleCentralReactingFoam solver is presented.The time derivatives in Equations ( 1)-( 4) are discretized by a simple Euler scheme.The finite volume method (FVM) is used in OpenFOAM to discretize Equations ( 1)-( 4).In the FVM, the solution domain is subdivided into a finite number of contiguous control volumes.The conservation equations are applied in integral form to each control volume, and then the volume integrals are transformed into surface integrals by using Gauss's theorem.Since the variable values are calculated only at the centroid of each control volume, interpolation is used to express variable values at the control volume surfaces.This is a critical part of the algorithm to approximate variables on the surfaces of the control volume, especially for the high-speed combustion with shocks, contact discontinuities, and rarefaction waves.The convective fluxes are reconstructed using a second-order (flux limiter-based) total variance diminishing (TVD) scheme [39,40].The discretization of viscous diffusion fluxes is approximated by using the second-order central differencing scheme.For the solution of the conservation Equations ( 1)-( 4), the following procedure is adopted [41]: (i) define convective fluxes for all Equations ( 1)-( 4) based on second-order (flux limiterbased) total variance diminishing (TVD) scheme [40], and viscous diffusion fluxes based on discretization by using the second-order central differencing scheme; (ii) solve explicit density Equation (1) and define new values ρ n+1 : (iii) using updated ρ n+1 , solve explicitly the inviscid momentum equation and define updated velocity u I , where the time derivative represents that due solely to inviscid fluxes: (iv) update the primitive variable , by solving implicit momentum Equation (2) with diffusive terms: (vi) solve explicit inviscid energy equation based on updated density and velocities to obtain an updated value for the total energy (viii) the new temperature T n+1 and pressure p n+1 are updated from equations of state.Note that equations with diffusive terms (i.e., (35), (36) and (38)) are formulated as implicit matrix equations.They are solved using the Preconditioned Bi-conjugate Gradient Stabilized (PBiCGStab) solver with preconditioners: Diagonal Incomplete Cholesky (DIC) for momentum equations, and Diagonal Incomplete LU (DILU) for species and energy equations [53].Explicit solvers are stable under the Courant-Friedrichs-Lewy (CFL) conditions [54], whereas Implicit solvers are unconditionally stable.The use of effective preconditioning techniques can enhance the convergence.
In summary, the new compressibleCentralReactingFoam solver is a combination of rho-CentralFoam and reactingFoam.The rhoCentralFoam was originally designed to solve compressible, supersonic gas flow for non-reactive systems.The reactingFoam can solve multicomponent reactive incompressible flows.The compressibleCentralReactingFoam solver developed in this work inherits the features of rhoCentralFoam and reactingFoam such that it can be used for supersonic combustion simulation with shock and rarefaction waves.

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.

Shock Tube Problem
A shock tube is a well-known benchmark problem (also known as the Sod shock tube problem [55]) for the validation of compressible codes.It is an attractive test case for model validation due to the simple one-dimensional geometry of the 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 problems 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 discontinuities, and rarefaction waves spreading in different directions (See Figure 1).
Dynamics 2024, 4, FOR PEER REVIEW conditions [54], whereas Implicit solvers are unconditionally stable.The use of effec preconditioning techniques can enhance the convergence.
In summary, the new compressibleCentralReactingFoam solver is a combination of CentralFoam and reactingFoam.The rhoCentralFoam was originally designed to solve c pressible, supersonic gas flow for non-reactive systems.The reactingFoam can solve mu component reactive incompressible flows.The compressibleCentralReactingFoam solver veloped in this work inherits the features of rhoCentralFoam and reactingFoam such th can be used for supersonic combustion simulation with shock and rarefaction waves.

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

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

•
Inert Multicomponent Shock Tube The first validation case is a multicomponent inert shock tube simulation without chemical source term in order to avoid shock-induced autoignition of the gas mixtur homogeneous premixed gas mixture of 20%H2/10%O2/70%Ar (in mole%) is used in

• 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 2 /10%O 2 /70%Ar (in mole%) is used in the simulation to verify the numerical accuracy of the model implementations for thermal and mass transport of a 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 ∆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 the 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. [56].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.
Dynamics 2024, 4, FOR PEER REVIEW 9 the 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 (lowpressure) 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. [56].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% H2/10% O2/70% Ar.The length of the tube is equal to  = 0.12 , with the initial conditions corresponding to the left and right sides as follows: (  ,   ,   ) = (0.072, 0, 7173) and (  ,   ,   ) = (0.18075, −487.34,35 594), where ,  , and  are measured in (/ 3 ), (/), and (), respectively.The length of the tube is discretized by 400 uniform elements with 0.01  timestep.The simulation used the reaction mechanism of Conaire et al. [57] for hydrogen, which consists of nine ( 9

• 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 2 /10% O 2 /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: (ρ l , u l , p l ) = (0.072, 0, 7173) and (ρ r , u r , p r ) = (0.18075, −487.34,35,594), where ρ, u , 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 µs timestep.The simulation used the reaction mechanism of Conaire et al. [57] for hydrogen, which consists of nine ( 9 The decay process of the initial multicomponent gas mixture and how the shock wave propagates in the tube were investigated in detail by comparing the profiles of inert and reactive mixtures of H 2 /O 2 /Ar. Figure 4 shows the simulation results for the temperature profiles of inert and reactive mixtures at various simulation times. shown in Figure 3 for 230  simulation time.Comparisons of the 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  are explained that Martínez-Ferrer et al. [34] sed a seventhorder 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 propagates in the tube were investigated in detail by comparing the profiles of inert and reactive mixtures of H2/O2/Ar.Figure 4 shows the simulation results for the temperature profiles of inert and reactive mixtures 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   = 8100 m/s and reaches the wall at around 75 .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   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  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 approaches 150 .The gas temperature approaches an equilibrium state as the simulation time increases to 220 .

Simulation of Ladenburg Jet Problem
Simulation of a well-known experiment of an underexpanded jet by Ladenburg et al. [58] is conducted to verify code accuracy for the shock formation of supersonic jets.A twodimensional axisymmetric simulation of the problem using rhoCentralFoam was reported previously [41,59].In the current work, a 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 in a frozen state, and hence, the gas mixture is considered be inert.The computational 3D mesh (see Figure 5) has a 5 mm inlet radius and 10 mm free surface radius, and is 30 mm in 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.72 • 10 5  , 316.6 / and 247.1 .The conditions at the free stream surface are 1.01 • 10 5 , 0 / and 297 .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.A 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 set by 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 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 at around 75 µ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 µ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 approaches 150 µs.The gas temperature approaches an equilibrium state as the simulation time increases to 220 µs.

Simulation of Ladenburg Jet Problem
Simulation of a well-known experiment of an underexpanded jet by Ladenburg et al. [58] is conducted to verify code accuracy for the shock formation of supersonic jets.A twodimensional axisymmetric simulation of the problem using rhoCentralFoam was reported previously [41,59].In the current work, a 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 in a frozen state, and hence, the gas mixture is considered be inert.The computational 3D mesh (see Figure 5) has a 5 mm inlet radius and 10 mm free surface radius, and is 30 mm in length.The following inlet conditions for pressure, axial velocity and temperature are used to ensure a sonic condition with Mach number of one (M = 1) at the inlet surface: 2.72•10 5 Pa, 316.6 m/s and 247.1 K.The conditions at the free stream surface are 1.01•10 5 Pa, 0 m/s and 297 K. Non-reflected boundary conditions are used at the outlet.The simulations are performed Dynamics 2024, 4 145 using a mesh with total hexahedral elements equal to 1.72 million.A 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 set by 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).This is shown in Figure 6b, in the free jet boundary above which the gas flow is at rest (i.e., M = 0).In Figure 6a     Figure 7 shows a comparison of the current simulation with experimental data of Ladenburg et al. [58] for density iso-lines.The simulation results have 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) show 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 the jet boundary and slip lines.Figure 7 shows a comparison of the current simulation with experimental data of Ladenburg et al. [58] for density iso-lines.The simulation results have 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 / 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,  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 the jet boundary and slip lines.

Simulation of DLR Scramjet Combustor
The combustion experiment of Waidmann et al. [60] was selected to validate the com-pressibleCentralReactingFoam 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), with the incoming air

Simulation of DLR Scramjet Combustor
The combustion experiment of Waidmann et al. [60] was selected to validate the com-pressibleCentralReactingFoam 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), with the incoming air flow at Mach number two (M = 2).The experimental system has been modeled using 2D and 3D computational studies by Oevermann [61] and Haung et al. [35], respectively.Oevermann [61] used RANS (i.e., k-ε turbulence model) with a laminar flamelet, whereas Huang et al. [35] used LES with a PaSR sub-grid scale combustion model to simulate the scramjet.Both computational studies produced reasonable agreement with the experimental data.In the current work, a 2D LES simulation is conducted with a PaSR sub-grid scale combustion model to validate the 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 the DLR test rig.The combustor has a height of 50 mm and width of 45 mm at air inflow.The divergence angle of the upper wall is 3 0 to compensate for the expansion of the boundary layer.m H2 ≈ 1.26 g/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 [62] is used.The reaction rates are expressed in Arrhenius form as described in Equations ( 24)-( 28) using the preexponent constant A f j , the temperature exponent β j , and the activation temperature T aj .
Dynamics 2024, 4, FOR PEER REVIEW 13 flow at Mach number two ( = 2).The experimental system has been modeled using 2D and 3D computational studies by Oevermann [61] and Haung et al. [35], respectively.Oevermann [61] used RANS (i.e., k-ε turbulence model) with a laminar flamelet, whereas Huang et al. [35] used LES with a PaSR sub-grid scale combustion model to simulate the scramjet.Both computational studies produced reasonable agreement with the experimental data.In the current work, a 2D LES simulation is conducted with a PaSR sub-grid scale combustion model to validate the 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.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 [62] is used.The reaction rates are expressed in Arrhenius form as described in Equations ( 24)-( 28) using the preexponent constant   , the temperature exponent   , and the activation temperature   .

• Simulation of Cold Flow DLR Combustor
Initially, a cold flow case 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 aerodynamic features such as shock waves, expansion flow at Mach number two ( = 2).The experimental system has been modeled using 2D and 3D computational studies by Oevermann [61] and Haung et al. [35], respectively.Oevermann [61] used RANS (i.e., k-ε turbulence model) with a laminar flamelet, whereas Huang et al. [35] used LES with a PaSR sub-grid scale combustion model to simulate the scramjet.Both computational studies produced reasonable agreement with the experimental data.In the current work, a 2D LES simulation is conducted with a PaSR sub-grid scale combustion model to validate the 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 the DLR test rig.The combustor has a height of 50 mm and width of 45 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   = 2:   = 730 /,   = 340 ,   = 10 5  with the mass rate of airflow ̇  = 1.0 /.The injector of fuel ( 2 ) is a slot in the center of the wedge-shaped strut with a height ∆ 2 = 0.26  that gives the surface area of injection ∆ 2 = 11.7  2 .The following boundary conditions are used to establish fuel injector Mach number  2 = 1:  2 = 1200 /,  2 = 250 ,  2 = 10 5  with the hydrogen mass flow ̇ 2 ≈ 1.26 /.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 [62] is used.The reaction rates are expressed in Arrhenius form as described in Equations ( 24)-( 28) using the preexponent constant   , the temperature exponent   , and the activation temperature   .

• Simulation of Cold Flow DLR Combustor
Initially, a cold flow case 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 aerodynamic features such as shock waves, expansion

• Simulation of Cold Flow DLR Combustor
Initially, a cold flow case 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 aerodynamic 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 an oblique shock reflected from the upper wall, whereas label (7) is the jet-shear layer formed from the Kelvin-Helmholtz instability.The oblique shocks that Dynamics 2024, 4 148 are reflected from the upper and lower walls interact with the jet-shear layer, to be reflected again without penetration, as shown Figure 10b.
Figure 11a shows the normalized contours of density gradient |/  |. 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 is 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  −  turbulence model did not generate shearlayer instabilities as reported by Oevermann [61].This finding confirms that  models are limited when predicting such phenomena and  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. [60].Figure 12a shows 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.Figure 11a shows the normalized contours of density gradient |∇ρ/ρ 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 is 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 − ε turbulence model did not generate shearlayer instabilities as reported by Oevermann [61].This finding confirms that RANS models are limited when predicting such phenomena and LES models are more appropriate in scramjet combustor simulations.
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 an oblique shock reflected from the upper wall, 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 be reflected again without penetration, as shown Figure 10b.
Figure 11a shows the normalized contours of density gradient |/  |. 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 is 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  −  turbulence model did not generate shearlayer instabilities as reported by Oevermann [61].This finding confirms that  models are limited when predicting such phenomena and  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. [60].Figure 12a shows 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.Figure 12 compares the cold flow simulation results for pressure distributions with the experimental measurements of Waidmann et al. [60].Figure 12a shows 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  • Simulation of Reacting Flow DLR Combustor Simulations were performed for reactive flow with a detailed hydrogen kinetic mechanism for experimental conditions of Waidmann et al. [60].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

• Simulation of Reacting Flow DLR Combustor
Simulations were performed for reactive flow with a detailed hydrogen kinetic mechanism for experimental conditions of Waidmann et al. [60].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 = 1500 K were implemented.Figure 13 compares a Schlieren image 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. [60].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 developed behind the hydrogen injection for the reactive flow case.

• Simulation of Reacting Flow DLR Combustor
Simulations were performed for reactive flow with a detailed hydrogen kinetic mechanism for experimental conditions of Waidmann et al. [60].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   = 1500 were implemented.Figure 13 compares a Schlieren image 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. [60].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 developed behind the hydrogen injection for the reactive flow case.Figure 14 shows the mean axial velocity along the center line of the DLR combustor.The numerical solution had good agreement with the experimental data.

Influence of Turbulence Models on Scramjet Combustion
A discussion on the influence of different turbulence models on supersonic combustion is presented in this section.Scramjet simulations for the DLR Scramjet experiment [60] were performed with the following turbulence models described in Section 2.2: SMG, LDkEqn, and WALE. Figure 15 compares the experimental data for temperature at three cross-sections along the combustor with the simulation results obtained with three turbulence models.The numerical results shown in Figure 15

Influence of Turbulence Models on Scramjet Combustion
A discussion on the influence of different turbulence models on supersonic combustion is presented in this section.Scramjet simulations for the DLR Scramjet experiment [60] were performed with the following turbulence models described in Section 2.2: SMG, LDkEqn, and WALE. Figure 15 compares the experimental data for temperature at three crosssections along the combustor with the simulation results obtained with three turbulence models.The numerical results shown in Figure 15 were obtained at 1 ms simulation time.The figure also shows the location of the temperature measurements in the DLR Scramjet combustor.The results close to the strut (i.e., x = 62 mm) show that the SMG and WALE models had similar profiles.Further downstream (i.e., x = 216 mm), all three models had similar profiles.Overall, the WALE model had slightly better agreement with the experimental data.

Influence of Turbulence Models on Scramjet Combustion
A discussion on the influence of different turbulence models on supersonic combustion is presented in this section.Scramjet simulations for the DLR Scramjet experiment [60] were performed with the following models described in Section 2.2: SMG, LDkEqn, and WALE. Figure 15 compares the experimental data for temperature at three cross-sections along the combustor with the simulation results obtained with three turbulence models.The numerical results shown in Figure 15 were obtained at 1 simulation time.The figure also shows the location of the temperature measurements in the DLR Scramjet combustor.The results close to the strut (i.e., x = 62 mm) show that the SMG and WALE models had similar profiles.Further downstream (i.e., x = 216 mm), all three models had similar profiles.Overall, the WALE model had slightly better agreement with the experimental data.[12,30,61,63].The modeling results of Berglund and Fureby [63] closely followed the experimental profile except at x = 62 mm.Overall, the present model had better agreement with the experimental data at all locations, but it did show relatively higher diffusivity at downstream locations (i.e., x = 108 mm and 216 mm).[12,30,61,63].The modeling results of Berglund and Fureby [63] closely followed the experimental profile except at x = 62 mm.Overall, the present model had better agreement with the experimental data at all locations, but it did show relatively higher diffusivity at downstream locations (i.e., x = 108 mm and 216 mm).Keys: symbols-experimental data [60]; lines-simulation.Color keys: black-current simulation with WALE model; red-Berglund and Fureby [63]; green-Genin and Menon [12]; blue-Zhang et al. [30]; purple-Oevermann [61].

Analysis of Turbulence-Combustion Interaction in Scramjet Combustor
Based on the analysis of different turbulence models, the WALE model was 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 the WALE model.Figure 17 shows the Mach number (Ma), rate of heat release per unit volume (Qdot), instantaneous temperature (T), and normalized contour lines of density gradient (|/ max |) obtained from the simulation.The shock structure, shown in Figure 17a,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  [60]; lines-simulation.Color keys: black-current simulation with WALE model; red-Berglund and Fureby [63]; green-Genin and Menon [12]; blue-Zhang et al. [30]; purple-Oevermann [61].

Analysis of Turbulence-Combustion Interaction in Scramjet Combustor
Based on the analysis of different turbulence models, the WALE model was 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 the WALE model.Figure 17 shows the Mach number (Ma), rate of heat release per unit volume (Qdot), instantaneous temperature (T), and normalized contour lines of density gradient (|∇ρ/ρ max |) obtained from the simulation.The shock structure, shown in Figure 17a,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 17b 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 15a.

Analysis of Turbulence-Combustion Interaction in Scramjet Combustor
Based on the analysis of different turbulence models, the WALE model was selected as the most suitable for 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 the WALE model.Figure 17 shows the Mach number (Ma), rate of heat release per unit volume (Qdot), instantaneous temperature (T), and normalized contour lines of density gradient (|/ max |) obtained from the simulation.The shock structure, shown in Figure 17a,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 17b 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 15a.Figure 19 shows the species mass fraction contours for , ,  2 and  2 .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 that correlates with the heat release rate.The 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 the species mass fraction contours for OH, H, HO 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 that correlates with the heat release rate.The field of H 2 O correlates with the field of instantaneous temperature, while the field of HO 2 correlates with OH and Qdot on the shear layer.
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 that correlates with the heat release rate.The 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 20 shows mixture fraction () and scalar dissipation rate ().The mixture fraction measures the local fuel/oxidizer ratio and is calculated as  = ( •   / 0 −   / + 1)/( + 1), where  is the equivalence ratio.The scalar dissipation rate is calculated as  = 2( • ), where  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 [65].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 20b.A time-scale analysis based on the estimation of the local Damköhler Number (Da) was also performed to examine the reactive zones in more detail.The Da number was quantified as the ratio between the turbulent mixing time scale  * (Equation (32)) and the chemical time scale   :  =  * /  [66,67].The Borghi diagram [68] shown in , where φ is the equivalence ratio.The scalar dissipation rate is calculated as χ = 2D(∇Z•∇Z), 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 [65].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 20b.
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 that correlates with the heat release rate.The 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 20 shows mixture fraction () and scalar dissipation rate ().The mixture fraction measures the local fuel/oxidizer ratio and is calculated as  = ( •   / 0 −   / 0 + 1)/( + 1), where  is the equivalence ratio.The scalar dissipation rate is calculated as  = 2( • ), where  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 [65].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 20b.A time-scale analysis based on the estimation of the local Damköhler Number (Da) was also performed to examine the reactive zones in more detail.The Da number was quantified as the ratio between the turbulent mixing time scale  * (Equation (32)) and the chemical time scale   :  =  * /  [66,67].The Borghi diagram [68] shown in A time-scale analysis based on the estimation of the local Damköhler Number (Da) was also performed to examine the reactive zones in more detail.The Da number was quantified as the ratio between the turbulent mixing time scale τ * (Equation (32)) and the chemical time scale τ c : Da = τ * /τ c [66,67].The Borghi diagram [68] shown in Figure 21a helps to recognize different regimes of turbulent combustion, where l is the characteristic size of flow geometry, l F is the laminar flame thickness, u ′ is characteristic velocity in the order of turbulence intensity, and S L is the laminar flame speed.Figure 21b,c show the Damköhler number (Da) on a logarithmic scale and the fraction of fine structure (γ * ), respectively.The bright white lines indicate the locations where Da = 1.The Damköhler number is relatively high (Da > 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 a combustion regime with fast chemistry.Thus, the fraction of reactive cells (γ * ) (see Equation (31) for details) is mostly small in this zone.In contrast, in the shear layer, where the rates of heat release (see Figure 17b) and scalar dissipation (see Figure 20b) have the highest values, the flame is controlled by finite rate chemistry with Da < 1 and 0 < γ * < 1.This combustion regime is characterized as a partially stirred reactor.In the core flow, small pockets where Da < 1 are also seen, indicating the combustion zone is in the "island formation" regime (see Figure 14.2 in Warnatz, et al. [66]).
tion (31) for details) is mostly small in this zone.In contrast, in the shear layer, where the rates of heat release (see Figure 17b) and scalar dissipation (see Figure 20b) have the highest values, the flame is controlled by finite rate chemistry with  < 1 and 0 <  * < 1.This combustion regime is characterized as a partially stirred reactor.In the core flow, small pockets where  < 1 are also seen, indicating the combustion is in the "island formation" regime (see Figure 14.2 in Warnatz, et al. [66]).

Conclusions
A computational validation of a new compressibleCentralReactingFoam code for modeling combustion in supersonic flows was presented.The validation cases included various complexities of high-speed flows with shocks, wave expansion, and turbulence-combustion interactions.Comparisons of the simulation results with data in the literature demonstrated the numerical fidelity of the solver for reactive multicomponent gas mixtures at supersonic conditions.In addition, a detailed analysis of gas dynamics and turbulence-combustion interaction under supersonic conditions was presented to gain a deeper understanding of the complex physical-chemical phenomena involved.A comparative analysis of the influence of three turbulent subgrid models on scramjet combustion was presented.The comparison of simulation results with the experimental data showed that the WALE subgrid scale is most suitable for scramjet modeling.The DLR scramjet simulation results showed that the flame was mainly stabilized in the combustor due to the formation of a subsonic bubble behind the strut.The supersonic flame structure was 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 revealed different regimes of turbulence combustion in the scramjet.In the core of the jet flow, the Da number was relatively high ( > 1), indicating that the chemical time scale was smaller than the turbulent mixing time scale.This means that the combustion regime was dominated by fast chemistry.The results also showed that the heat release rate and the scalar dissipation rate had the highest values in the shear layer, and the flame was stabilized due to finite rate chemistry with  < 1 and 0 <  * < 1.Thus, the new solver presented in this paper can

Conclusions
A computational validation of a new compressibleCentralReactingFoam code for modeling combustion in supersonic flows was presented.The validation cases included various complexities of high-speed flows with shocks, wave expansion, and turbulence-combustion interactions.Comparisons of the simulation results with data in the literature demonstrated the numerical fidelity of the solver for reactive multicomponent gas mixtures at supersonic conditions.In addition, a detailed analysis of gas dynamics and turbulence-combustion interaction under supersonic conditions was presented to gain a deeper understanding of the complex physical-chemical phenomena involved.A comparative analysis of the influence of three turbulent subgrid models on scramjet combustion was presented.The comparison of simulation results with the experimental data showed that the WALE subgrid scale is most suitable for scramjet modeling.The DLR scramjet simulation results showed that the flame was mainly stabilized in the combustor due to the formation of a subsonic bubble behind the strut.The supersonic flame structure was 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 revealed different regimes of turbulence combustion in the scramjet.In the core of the jet flow, the Da number was relatively high (Da > 1), indicating that the chemical time scale was smaller than the turbulent mixing time scale.This means that the combustion regime was dominated by fast chemistry.The results also showed that the heat release rate and the scalar dissipation rate had the highest values in the shear layer, and the flame was stabilized due to finite rate chemistry with Da < 1 and 0 < γ * < 1.Thus, the new solver presented in this paper can be used for high-speed simulations of multicomponent reactive gas mixtures in supersonic combustors.

Figure 1 .
Figure 1.Shock tube and spreading shock wave (red), the contact surface (black), and expansion fan (blue).
) chemical species (H2, O2, H, O, OH, HO2, H2O2, H2O, ), 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  , and temperature:   ,   , and  are shown in Figure3for 230  simulation time.Comparisons of the current solution (solid lines) in Figure3with 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  are explained that Martínez-Ferrer et al.[34] sed a seventhorder accurate Weighted Essentially Non-Oscillatory (WENO) scheme to discretize the non-linear advective terms.
) chemical species (H 2 , O 2 , H, O, OH, HO 2 , H 2 O 2 , H 2 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 µs simulation time.Comparisons of the 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 µs are explained that Martínez-Ferrer et al. [34] sed a seventh-order accurate Weighted Essentially Non-Oscillatory (WENO) scheme to discretize the non-linear advective terms.

Figure 4 .
Figure 4. Temperature profiles of shocks/detonation waves for (a) inert and (b) reactive mixture at various times.

Figure 4 .
Figure 4. Temperature profiles of shocks/detonation waves for (a) inert and (b) reactive mixture at various times.
, 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.thecircles 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 5 .
Figure 5. Computational 3D mesh for simulation of underexpanded jet of Ladenburg experiment.Inflow is shown by the yellow surface, the free stream inlet surface is violet, and the free stream surface is red.

Figure 6 .
Figure 6.(a) Schematics of an underexpanded air jet.Keys: 1 expansion waves, 2 oblique shock, 3 normal shock, 4 reflected shock, 5 reflected expansion waves.(b) Contours of Mach numbers at the plane through the center of computational domain (see Figure 5).

Figure 7
Figure7shows a comparison of the current simulation with experimental data of Ladenburg et al.[58] for density iso-lines.The simulation results have good agreement with the experimental data, especially the location and height of the Mach disk.The ex-

Figure 5 .
Figure 5. Computational 3D mesh for simulation of underexpanded jet of Ladenburg experiment.Inflow is shown by the yellow surface, the free stream inlet surface is violet, and the free stream surface is red.

Figure 5 .
Figure 5. Computational 3D mesh for simulation of underexpanded jet of Ladenburg experiment.Inflow is shown by the yellow surface, the free stream inlet surface is violet, and the free stream surface is red.

Figure 6 .
Figure 6.(a) Schematics of an underexpanded air jet.Keys: 1 expansion waves, 2 oblique shock, 3 normal shock, 4 reflected shock, 5 reflected expansion waves.(b) Contours of Mach numbers at the plane through the center of computational domain (see Figure 5).

Figure 6 .
Figure 6.(a) Schematics of an underexpanded air jet.Keys: 1 expansion waves, 2 oblique shock, 3 normal shock, 4 reflected shock, 5 reflected expansion waves.(b) Contours of Mach numbers at the plane through the center of computational domain (see Figure 5).

Dynamics 2024, 4 ,
FOR PEERREVIEW  12    barrel shock and reflected shock).The experimental data in Figure7(i.e., lower panel) show 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

Figure 7 .
Figure 7. Density distribution of the supersonic jet measured (reprinted with permission from Ladenburg [58]) (lower panel), compared with current computational results (upper panel).The values of the density contours are provided in kg/m 3 .

Figure 7 .
Figure 7. Density distribution of the supersonic jet measured (reprinted with permission from Ladenburg [58]) (lower panel), compared with current computational results (upper panel).The values of the density contours are provided in kg/m 3 .
The following boundary conditions for air are established to ensure the Mach number M air = 2: u air = 730 m/s, T air = 340 K, p air = 10 5 Pa with the mass rate of airflow .m air = 1.0 kg/s.The injector of fuel (H 2 ) is a slot in the center of the wedge-shaped strut with a height ∆y H2 = 0.26 mm that gives the surface area of injection ∆s H2 = 11.7 mm 2 .The following boundary conditions are used to establish fuel injector Mach number M H2 = 1: u H2 = 1200 m/s, T H2 = 250 K, p H2 = 10 5 Pa with the hydrogen mass flow.

Figure 8
shows the geometry of the DLR test rig.The combustor has a height of 50 mm and width of 45 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   = 2:   = 730 /,   = 340 ,   = 10 5  with the mass rate of airflow ̇  = 1.0 /.The injector of fuel ( 2 ) is a slot in the center of the wedge-shaped strut with a height ∆ 2 = 0.26  that gives the surface area of injection ∆ 2 = 11.7  2 .The following boundary conditions are used to establish fuel injector Mach number  2 = 1:  2 = 1200 /,  2 = 250 ,  2 = 10 5  with the hydrogen mass flow ̇ 2 ≈ 1.26 /.

Figure 9 .
Figure 9. Computational 2D mesh for simulation of turbulence combustion in DLR.

Figure 9 .
Figure 9. Computational 2D mesh for simulation of turbulence combustion in DLR.

Figure 9 .
Figure 9. Computational 2D mesh for simulation of turbulence combustion in DLR.

Figure 10 .
Figure 10.(a) Cold flow Schlieren image of the channel with hydrogen injection (reprinted with permission from [60,61]).(b) Computational density gradient contours corresponding to the experimental image in (a).Labels 1 to 7 in (b) refer to various flow structures discussed in the text.

Figure 10 .
Figure 10.(a) Cold flow Schlieren image of the channel with hydrogen injection (reprinted with permission from [60,61]).(b) Computational density gradient contours corresponding to the experimental image in (a).Labels 1 to 7 in (b) refer to various flow structures discussed in the text.

Figure 10 .
Figure 10.(a) Cold flow Schlieren image of the channel with hydrogen injection (reprinted with permission from [60,61]).(b) Computational density gradient contours corresponding to the experimental image in (a).Labels 1 to 7 in (b) refer to various flow structures discussed in the text.

15 Figure 12 .
Figure12compares the cold flow simulation results for pressure distributions with the experimental measurements of Waidmann et al.[60].Figure12ashows the distribution of pressure on the lower (solid line) and upper (dashed line) walls along with experimental data obtained on the lower wall (symbols).Figure12bcompares the pressure along the center line with the experimental data.The model predictions are able to capture the pressure profiles observed experimentally.Dynamics 2024, 4, FOR PEER REVIEW 15

Figure 13 .
Figure 13.(a) Schlieren photograph of the channel flow with hydrogen injection (reprinted with permission from [60,61]); (b) computational density gradient contours corresponding to the experimental conditions.

Figure 13 .
Figure 13.(a) Schlieren photograph of the channel flow with hydrogen injection (reprinted with permission from [60,61]); (b) computational density gradient contours corresponding to the experimental conditions.

Figure 14 Figure 14 .
Figure 14 shows the mean axial velocity along the center line of the DLR combustor.The numerical solution had good agreement with the experimental data.Dynamics 2024, 4, FOR PEER REVIEW 16 were obtained at 1 simulation time.The figure also shows the location of the temperature measurements in the DLR

Figure 14 .
Figure 14.Mean axial velocity profile along the middle of the channel at y = 25 mm.Keys: symbolsexperimental data [60]; line-current simulation.

Figure 14 .
Figure 14.Mean axial velocity profile along the middle of the channel at  = 25 .Keys: symbols-experimental data [60]; line-current simulation.

Figure 16
Figure16compares the current WALE model simulation results with other DLR Scramjet simulation data reported in the literature[12,30,61,63].The modeling results of Berglund and Fureby[63] closely followed the experimental profile except at x = 62 mm.Overall, the present model had better agreement with the experimental data at all locations, but it did show relatively higher diffusivity at downstream locations (i.e., x = 108 mm and 216 mm).

Figure 16
Figure16compares the current WALE model simulation results with other DLR Scramjet simulation data reported in the literature[12,30,61,63].The modeling results of Berglund and Fureby[63] closely followed the experimental profile except at x = 62 mm.Overall, the present model had better agreement with the experimental data at all locations, but it did show relatively higher diffusivity at downstream locations (i.e., x = 108 mm and 216 mm).

Figure 17 .
Figure 17.Instantaneous distribution of: (a) Mach number; (b) normalized heat release rate   = /10 10 ; (c) temperature; (d) normalized contour lines of density gradient.The direction of diffusive spreading of fuel and oxidizer helps to analyze turbulencechemistry interactions.It can be estimated by the interaction between the gradients of fuel (  ) and oxygen (  ) mass fractions.Takeno and co-workers [64] 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 [64].The TFI is defined as   = (  •   )/ (|  | |  |) .Positive and negative values of TFI indicate non-premixed or premixed modes, respectively [64]. Figure 18 shows the TFI calculated for the DLR Scramjet conditions.Pockets of non-premixed zones shown in red are surrounded by large premixed zone shown in dark blue.

Figure 17 . 18 Figure 18 .
Figure 17.Instantaneous distribution of: (a) Mach number; (b) normalized heat release rate Qdot Norm = Qdot/10 10 ; (c) temperature; (d) normalized contour lines of density gradient.The direction of diffusive spreading of fuel and oxidizer helps to analyze turbulencechemistry interactions.It can be estimated by the interaction between the gradients of fuel (∇Y F ) and oxygen (∇Y O ) mass fractions.Takeno and co-workers [64] 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 [64].The TFI is defined as G FO = (∇Y F •∇Y O )/(|∇Y F ||∇Y O |).Positive and negative values of TFI indicate non-premixed or premixed modes, respectively [64]. Figure 18 shows the TFI calculated for the DLR Scramjet conditions.Pockets of non-premixed zones shown in red are surrounded by large premixed zone shown in dark blue.Dynamics 2024, 4, FOR PEER REVIEW 18

Figure 20
Figure20shows mixture fraction (Z) and scalar dissipation rate ( χ).The mixture fraction measures the local fuel/oxidizer ratio and is calculated asZ = (φ•Y F /Y F0 − Y O /Y O0 + 1)/ (φ + 1), where φ is the equivalence ratio.The scalar dissipation rate is calculated as χ = 2D(∇Z•∇Z), 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[65].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 Figure20b.