# Numerical Investigation on Flame Stabilization of Cavity-Based Scramjet Combustor Using Compressible Modified FGM Model

Authors: Qi Zhang; Weibing Zhu; Dongchao Yang; Hong Chen
Journal: International Journal of Chemical Engineering
 (2022)

---

## Abstract

The technology of flame stabilization on cavity-based scramjet combustor has great significance in the field of future spacecraft. In this paper, a compressible modified FGM model was established based on the idea of the flamelet model, which was adopted to simulate the unsteady combustion process of the hydrogen transverse jet in the upper cavity of the scramjet. The results show that the compressible modified FGM model can accurately reflect the flow field and the propagation process of the flame in the supersonic cavity, and can capture the fine shock structure in the flow field. The coupling effect of shock waves and shear layer cause the shear layer to quickly destabilize, resulting in the turbulence effect, which promotes the mixing of air and fuel. The boundary layer separation at the upper wall of the combustion chamber will reduce the stability of the shear layer.


---

## Body

## 1. Introduction

Hypersonic vehicle usually refers to a flight power device with a flight speed of more than Mach 5 [1, 2]. As a power device capable of providing high-speed flight, the hypersonic vehicle has been widely recognized and valued by countries all over the world [3, 4]. In order to achieve hypersonic flight, currently, the United States, Russia, France, the United Kingdom, Japan, and other countries have, respectively, proposed hypersonic propulsion technical solutions, and studied combined power engines based on scramjet (supersonic combustion ramjet). Among them, it is of great significance to break through and master the key technology of the scramjet engine, which is the necessary prerequisite for realizing hypersonic flight [5–13].The cavity-based scramjet engine has been widely recognized since it was put forward [14, 15], but the flame stabilization mechanism of the cavity has not been fully clarified so far. At present, the research on the flame stabilization of the cavity is mostly based on experimental research. Byrne et al. [16] and Micka et al. [17] carried out experiments on fuel injection in the cavity and outside, respectively, and observed the supersonic combustion process with OH/CH-PLIF for the first time. Song et al. [18] investigated lean blowoff characteristics of the cavity in an ethylene-fueled model scramjet combustor under inflow conditions of Ma = 2.52. Sun et al. [19, 20] and Fan et al. [21] used OH-based spontaneous radiation imaging and PLIF technology to study the flow field of transverse hydrogen injection upstream of the cavity. Meanwhile, in the experiment, they also photographed the distribution of flame structure, which has a high reference value for studying the formation mechanism of the flame upper cavity. Trudgian et al. [22] experimentally investigated the upstream wall of a scramjet combustor cavity, which was inclined, aiming to reduce cavity base drag and remove the superfluous cavity vortex structures. Rong et al. [23] carried out optical research on the ignition process of ethylene-fueled cavity-based scramjet with inlet velocity of Ma = 2.92, and recorded flame chemiluminescence and schlieren photography. However, limited by experimental conditions and insufficient understanding of supersonic combustion, it is difficult to obtain details of the unsteady combustion process from ignition to flame stabilization. Numerical simulation is suitable for adoption to study the flow characteristics of the cavity [24–26], which can make up for the limitations and shortcomings of the experiment.Due to the limitations of the turbulence model, the combustion model, and other factors, the numerical research mainly focuses on the analysis of pressure, temperature, composition, and Mach number in combustion chambers with different structures, seeking for methods to improve the efficiency of combustion and total pressure recovery coefficient. However, little attention has been paid to the flame form in the combustion chamber. The accuracy and precision have not yet reached the degree of replacing the experiment, and much work still needs to be further improved [27–30]. Therefore, many scholars strive to find and improve the numerical model to enhance the simulation accuracy and realize the accurate simulation of the flame form. Among them, the FGM (Flamelet-Generated Manifolds) model shows great potential in simulating flame form and can effectively reduce the computational cost [31, 32], which is favored by some scholars. Vreman et al. [33] employed FGM model to achieve direct numerical simulations of a lean premixed turbulent Bunsen flame with hydrogen addition. The results show that the FGM model can effectively improve the computational efficiency and realize the accurate simulation of the flame form. Honzawa et al. [34] adopted the FGM model and combined with the LES (large eddy simulation) model to simulate the combustion process of the NH3/CH4/air swirl burner, the simulation results of which were more consistent with the experiment. Liu et al. [35] utilized the FGM model to simulate flame diffusion in high-hydrogen micro-mixing model burner, which can accurately capture the position of internal and external shear layers of flame. At present, the FGM model is mostly used for the simulation of subsonic combustion. When simulating the supersonic combustion studied in this paper, the compressible effect of the flow field should be fully considered [36–38]. Therefore, it is necessary to modify the FGM model for compressibility.In this paper, the compressible modification of the FGM model, which can simulate the flame form more realistically, is achieved. A compressible modified FGM model suitable for supersonic combustion is established by improving the method of solving compressible flow field parameters. The integrity and stability of the model are verified by simulating the combustion flow field of a typical DLR scramjet combustion chamber. Combining with hybrid RANS/LES turbulence model, the unsteady combustion process of hydrogen transverse jet in a cavity-based scramjet is numerically simulated. Compared with the experimental results of Fan et al., it is found that the simulation can accurately capture the formation of flame, distribution of shock waves, and the phenomenon of “flame lifting” in a single cavity combustion chamber. In addition, the combustion stabilization mechanism of the cavity structure is analyzed by combining the phenomenon of shock waves and boundary layer, which provides a reliable theoretical basis for the design of the scramjet engine.


## 2. Numerical Approach

### 2.1. Governing Equation

According to the conservation of mass, momentum, and energy, the governing equations of gas-phase flow are established. The detailed derivation of the equations is based on the work of Versteeg and Malalasekera [39], where only the final results in matrix form are given. In the actual simulation of combustion, the Navier–Stokes equation system containing multi-species and conservative chemical reactions is usually used, and the gas transport process includes three phenomena, namely: diffusion phenomenon of species transport caused by concentration gradient, viscosity phenomenon of momentum transport caused by velocity gradient, and heat conduction phenomenon of energy transport caused by temperature gradient and species diffusion. In order to calculate the momentum, species, and energy transport intensity in the combustion flow field, the transport coefficients of each species need to be given. The governing equations of the gas-phase are written as follows:(1)∂Q∂t+∂E−Ev∂x+∂F−Fv∂y+∂G−Gv∂z=H,(2)Q=ρρuρvρwρeρYiE=ρuρuu+pρuvρuwuρe+pρuYiEv=0τxxτxyτxzuτxx+vτxy+wτxz−qxρiDim∂Yi∂x,(3)F=ρvρvuρvv+pρvwvρe+pρvYiFv=0τyxτyyτyzuτxy+vτyy+wτyz−qyρiDim∂Yi∂y,(4)G=ρwρwuρwvρww+pwρe+pρwYiGv=0τzxτzyτzzuτzx+vτzy+wτzz−qzρiDim∂Yi/∂zH=00000ωi.In Equations (1)–(4): i = 1, 2, …, Ns−1, Ns is the number of gas species. Q is an independent variable that is conserved; E, F, G, Ev, Fv, and Gv are, respectively, the inviscid and viscous fluxes in the x, y, and z directions of the Cartesian coordinate system; u, v , and w are the partial velocities along the x, y, and z directions of the coordinate axes; p is the pressure of the gas-phase; ρi is the density of the species i; ρ is the density of the gas-phase mixture; Yi is the mass fraction of the species i; ωi is the mass generation rate of species i. τi,j is the component of the viscous stress, as follows:(5)τxx=−23μ∇⋅U+2μ∂u∂x,τyy=−23μ∇⋅U+2μ∂u∂y,τzz=−23μ∇⋅U+2μ∂u∂z,τxy=τyx=μ∂u∂y+∂v∂x,τyz=τzy=μ∂v∂z+∂w∂y,τxz=τzx=μ∂w∂x+∂u∂z.μ is the laminar viscosity coefficient, given by the Surferland formula:(6)μ=μ0TT01.5T0+TsT+Ts,μ0 is the reference viscosity coefficient, T0 is the reference temperature, and Ts is the characteristic temperature. For the air, they are 17.161 × 10–6 (Pa∙s), 273.16 K, and 124 K, respectively. qi (i = x, y, z) represents the energy flux caused by thermal conductivity and species diffusion, as follows:(7)qx=−λ∂T∂x−ρ∑i=1NSDimhi∂Yi∂x,qy=−λ∂T∂y−ρ∑i=1NsDimhi∂Yi∂y,qz=−λ∂T∂z−ρ∑i=1NsDimhi∂Yi∂z.Dim is the diffusion coefficient of speciesi, as follows:(8)Dim=1−Xi∑i,i≠jXj/Dij.In this paper, the thermal conductivityλ and the species diffusion coefficient Dim are given by the laminar Prandtl number Pr and Schmidt number Sc: λ = cp∙μ/Pr, Dim = μ/(ρ∙Sc).The mole fraction of the species is calculated as follows:(9)Xi=Yi/Wi∑j=1NsYj/Wj.For calorically perfect gas, when translational kinetic energy and rotational energy are fully excited, the specific heat capacity and specific heat ratio are constant, and the internal energy is only a linear function of temperature. For thermally perfect gas involving combustion, the enthalpy is no longer a linear function of temperature. The thermodynamic function of each gas species is usually given in the form of a polynomial temperature, and coefficients of the polynomial are obtained by least square method fitting. In this paper, the polynomial form given by JANAF database [40] is used to calculate the cp,i and enthalpy hi of gas at a given temperature, as follows:(10)cp,i=A+B×T+C×T2+D×T3+ET4,hi−hi,298.15=A×T+B×T22+C×T33+D×T44−ET+F−H.The enthalpy and internal energy of the gas-phase mixture are defined as follows:(11)hT=∑i=1NsYihiT,eT=∑i=1NsYihiT−pρ.In addition, the multi-species gas-phase mixture satisfies the ideal gas equation of state, and the pressure of the mixture can be obtained according to Dalton partial pressure theorem:(12)p=R0T∑i=1NsρiWi.R0 is the general gas constant, 8.3145 kJ/(kg K). Ri is the gas constant of species i. The gas constant of the gas-phase mixture is defined as follows:(13)R=∑iYiRi=∑iYiR0Wi.The specific heat capacity at constant volume of gas-phase mixture is calculated as follows:(14)cp=∑i=1NsYicpiThe molar mass of the gas-phase mixture is calculated as follows:(15)W=∑XiWi=∑YiWi−1,According to equation (11), energy (internal energy e or enthalpy h) is a nonlinear equation of temperature T. After obtaining the energy field according to the energy equation, it is necessary to solve the temperature field iteratively. Generally, Newton iteration method is adopted [41].


### 2.2. FGM Model

The FGM (flamelet-generated Manifolds) model is developed from incompressible premixed combustion [42], which is essentially a model developed based on the idea of the flamelet model [43]. The idea of the FGM model comes from ILDM (Intrinsic low-dimensional Manifolds) [44]. ILDM simplifies the chemical reaction into fast reaction and slow reaction, and only considers the manifold with the smallest time scale of chemical reaction (Manifold). Therefore, the simulation effect of fast-reaction and high-temperature section is better, while the effect of slow-reaction and low-temperature section (The effects of diffusion and convection cannot be ignored.) is poor. In order to deal with the existence of diffusion and convection, the FGM model adopts the premixed flame calculation of 1D laminar flow based on nontension in both low-temperature and high-temperature sections. In other words, low-dimensional manifolds are created by using the solution of 1D flamelet equation, and the effect of diffusion and convection is taken into account, which greatly improves the simulation of low-temperature section and high-temperature section.The FGM model can be parameterized by any form of laminar flame, which can be calculated by 1D steady premixed flamelet or diffused flamelet model. According to the different laminar flame forms, the FGM model is divided into two sub-models, premixed model and diffusion model. Generally, the premixed FGM model is used in the turbulent partially premixed flame dominated by premixed combustion. The combustion of the scramjet studied in this paper is a nonpremixed turbulent partially premixed flame, which needs to adopt the diffusion FGM model, and the required query table is obtained through the calculation of one-dimensional steady-state diffusion flame.


#### 2.2.1. Establishment of Database

FGM is developed based on the flamelet model. The establishment method of its database is the same as that of the flamelet model. It needs to solve the state parameters of the laminar diffusion flame with one-dimensional structure first. This process can be realized by solving the state equations of the laminar flow collision flame (as shown in Figure1).Figure 1
Schematic of the planar laminar counter flow flame structure.Mixture fractionZ is an important concept in diffusion combustion. The mixture fraction represents the scalar of the fuel to oxidant mixing ratio, which is a function of space and time, and can be defined in many ways [45–49]. Masri gives a definition form based on the conservation of local elements [50], as follows:(16)Z=∑iyi1−yi0∑jyj1−yj0.In equation (16): yi0 and yi1 are the mass fractions of element i at the entrance of oxidant and fuel, respectively.According to the definition ofZ given in equation (16), the oxidant inlet Z = 0; the fuel inlet Z = 1. It can be guaranteed that Z is a conserved scalar without a source term. The transport equation is written as follows:(17)∂ρZ∂t+∇⋅ρUZ=∇⋅ρDZ∇Z.In equation (17): DZ is the diffusion coefficient of Z.The state of collision flame is also described by N-S equations given in Section2.1. Considering the geometric conditions and characteristics of the flame structure, equation (1) can be simplified appropriately. A new coordinate system with the normal direction of the flamelet as x is established, which is perpendicular to the flamelet and the other two directions are x1 and x2. The simplified system of equations is written as follows [48]:(18)∂ρu∂x=−ρK,(19)∂ρYi∂x=1Lei∂∂xλcp∂Yi∂x+ωi−ρKYi,i=1,…,Ns−1,(20)∂ρuh∂x=∂∂xλcp∂h∂x+∑Nshiλcp1Lei−1∂Yi∂x−ρKh.In Equations (18)–(20): x, u, and K are the space coordinate, velocity component, and flame stretch rate perpendicular to the flamelet, respectively. For the one-dimensional situation discussed in this section, K is the rate of change of u along the x direction, that is, K = ∂u/∂x.In addition, the transport equation for the mixture fractionZ is written as follows:(21)∂ρuZ∂x=∂∂xλcp∂Z∂x−ρKZ.Combined with the analysis of equations (16) and (21), Z (x) is a pure monotone function along the axis of the colliding flame (normal direction of the flame surface). Therefore, Z can be used to parameterize the flow field variables.In actual nonpremixed combustion, in order to express the state parameters in the form ofφZ in Z, the governing equation in the physical space (Figure 2) needs to be converted into the Z space (Figure 3). This process can be achieved by coordinate transformation in the type of Crocco [42].Figure 2
Variables in physical space.Figure 3
Variables in Z space.The transport equations of species and temperature are transformed as follows [51]:(22)∂Yi∂t−χ2∂2Yi∂Z2−ωiρ=0,(23)∂T∂t−χ2∂2T∂Z2+1ρcp∑i=1Nshiωi−1ρcp∂p∂t=0.In Equations (22) and (23): The dimension of Scalar Dissipation Rate χ is s−1, which characterizes the effect of chemical nonequilibrium effects. The equation of definition is written as follows:(24)χ=2D∇Z2.Equations (22) and (23) are the one-dimensional flamelet equations in Z space, and there is no convection term. After the initial and boundary conditions are determined, the distribution of chemical thermodynamic parameters in flamelet can be obtained by solving these two equations. The solution of the laminar diffusion flame represented by t, Z, and χ is defined as follows:(25)φ=φt,Z,χ.In the FGM model,χ is defined as Progress-variable C that characterizes the degree of combustion. In the case of steady-state combustion, equation (25) can be transformed into the following equation:(26)φ=φZ,C.In general, equation (26) is not directly solved. Instead, according to certain rules, equation (26) is made into a unified format database, which is known as the laminar FGM database.When the thickness of the flamelet is less than the minimum vortex scale (Kolmogorov scale) of turbulent flow (as shown in Figure4), the flamelet can still maintain a laminar state. In this case, the turbulent flame structure can be represented by an ensemble embedded in the turbulent flow field with a one-dimensional wrinkled laminar flamelet. The thermodynamic state parameters of the average flame structure in the turbulent field are given by Statistical average of the laminar flamelet, as follows:(27)φ¯=∫∫φZ,CPZ,CdZdC.Figure 4
Kolmogorov vortex and flamelet in the wrinkled laminar flamelet.In equation (27): φ¯ is the average thermodynamic state parameter required in turbulent flame; P (Z, C) is the Joint Probability Density Function (JPDF) of Z and C. Since φ has been obtained by the solution of equation (26), P (Z, C) only needs to know to calculate φ¯.It is generally assumed thatZ and C are independent of each other, so JPDF can be represented as the product of the Marginal Density-weighted Filter PDF (FPDF):(28)PZ,C=PZPC.By introducing the mean valueZ¯ and variance Z″2¯ of Z, and substituting equation (28) into equation (27), the chemical thermodynamic state of turbulent combustion can be calculated as follows:(29)φ¯=φ¯Z¯,Z″2¯,C.By integrating the calculated multi-group laminar database, the turbulent database under the condition of differentC (progress variables) can be obtained, and then the turbulent database can be made into a lookup table according to certain rules, which is known as the turbulent FGM database.The turbulent FGM database can be combined with CFD solvers to solve the combustion flow field. After the solver calculates the distribution ofZ and C, it uses the lookup table to obtain the components and other required parameters to complete the solution. Therefore, for the FGM model, in addition to solving the continuity equation, and the momentum equation and transport equation of Z, it is also necessary to solve the scalar transport equation of C:(30)∂ρC∂t+∇⋅ρUC=∇ρDC∇C+ωC.In equation (30): DC is the effective transport coefficient of the schedule variable; ωC is the chemical reaction source term of C.


#### 2.2.2. Compressible Modification of FGM Model

As far as the FGM model is concerned, the thermochemical state parameters in the flamelet database include but are not limited toZ, T, ρ, Yi, and ωi. The incompressible combustion gas state parameters T, ρ, and p are generally decoupled (regardless of the buoyancy flow). T and ρ can be calculated in advance and obtained from the lookup table. However, for compressible turbulent combustion in supersonic conditions, U, T, ρ, and p are strongly coupled and limited by the ideal gas state equation. A change in kinetic energy can cause a change in temperature, and the temperature in the database cannot reflect this coupling effect. If the T and ρ in the flow field are the same as in the case of incompressibility, the parameters of the flow field will introduce a large error. In addition, the flamelet database is based on the incompressible flamelet equation. The state points represented in the database are different from the supersonic combustion state. Therefore, in order to apply the traditional FGM model to supersonic combustion, it is necessary to modify the FGM model given above.At present, the modification methods of the FGM model can be summarized into two categories: one is to improve the generation method of the database, add more compressible factors to the flame control equation or database, and strive to build a more accurate database; the other is to adjust the query mode of the database and only use part of the parameters in the database, or to make some transformation of the parameters to match the state in supersonic combustion [52, 53]. Among the numerous correction algorithms, the method proposed by Oevermann [54] in 2000 has attracted great attention due to its simplicity, feasibility, and reasonable results. This method is a shock capture algorithm, which only uses the mass fraction of species in the database without temperature or other parameters. The temperature is obtained by the energy equation. In terms of computational efficiency, this processing method only solves one more energy equation than the original flamelet model, which still retains the advantages of the flamelet model. In fact, the huge amount of calculation and the rigidity of the chemical reaction source term brought by the multi-species equation in the combustion solution are difficult to deal with. Oevermann solved these two problems skillfully with the help of the flamelet model and the database query algorithm.Based on the above analysis, the flamelet model modification algorithm proposed by Oevermann was adopted to generate the compressible modified FGM model and conduct a numerical study on supersonic combustion. The specific implementation process is as follows:(1)
To solve the energy distributionh(T) in the flow field from the energy transport equation;(2)
To query the composition distributionYi (T) in the flow field through the database;(3)
To substituteh(T), Yi (T) and the thermodynamic parameters cp,i (T) of each species (equation (10)) into Equation (11). The temperature field T is obtained by the Newton iteration method.Compressibility of gas should be taken into full consideration when simulating the flow field supersonic combustion. With the compressible modification of the FGM model, it is not necessary to query the temperatureT from the FGM database when solving the temperature field, and the database only provides the parameters of the species. At the same time, a step of solving the energy equation is added. By substituting the solution of the energy equation and the species parameters obtained from database query into Equation (11) for iteration, the error caused by directly querying the database to obtain temperature T can be eliminated. Compressible modified FGM model can achieve accurate simulation of supersonic combustion while retaining the advantages of simulating the flame form.


#### 2.2.3. Verification of Compressible Modified FGM Model

In order to verify the stability and integrity of the model, the traditional flamelet model and the modified FGM model are used to numerically simulate the combustion flow field of the typical DLR scramjet combustion chamber. The structure of the combustion chamber is shown in Figure5 [55]. The preheated air expands into the combustion chamber through the nozzle. The conditions of air inflow and hydrogen injection remain the same as the experimental ones, as shown in Table 1. The combustion chamber has a three-dimensional structure and relatively dense nozzle holes are arranged in the span direction. After comparing the two-dimensional and three-dimensional simulations in detail, Luo et al. [56] found that in addition to the three-dimensional effect near the nozzle, the main flow field presents obvious two-dimensional characteristics. Therefore, the RANS method is adopted in this section to simulate the two-dimensional simplified structure, and the nozzle is equivalent to a slit in the same area. For the reaction mechanism of hydrogen and air used to generate laminar flamelet database, a 19-step reaction mechanism of 9 species (H2, O2, N2, H2O, OH, H, O, HO2, H2O2) proposed by Jachimowski [57] is adopted, as shown in Table 2.Figure 5
Schematic diagrams of the DLR scramjet combustion chamber.Table 1
Experiment conditions of the DLR scramjet combustion chamber.
AirH2Ma2.01.0U/(m/s)7301200T/K340250p/MPa0.10.1ρ/(kg/m3)1.0020.097YO20.2320YN20.7360YH2O0.0320YH201k/(m2/s2)102400ε/(1/s)722462963Table 2
H2-O2-N2 9-specie 19-step reaction mechanism [57].
ReactionA(mole/(s∙cm3))E(cal/mol)B (-)H2 + O2 ↔ HO2+H1.0 × 1014560000H + O2 ↔ OH + O2.6 × 1014168000O + H2 ↔ OH + H1.8 × 101089001OH + H2 ↔ H + H2O2.2 × 101351500OH + OH ↔ O + H2O6.3 × 101210900H + OH + M ↔ H2O + M2.2 × 10220−2H + H + M ↔ H2 + M6.4 × 10170−1H + O + M ↔ OH + M6.0 × 10160−0.6H + O2 + M ↔ HO2 + M2.1 × 1015−10000O + O + M ↔ O2 + M6.0 × 1017−18000HO2 + H ↔ OH + OH1.4 × 101410800HO2 + H ↔ H2O + O1.0 × 101310800HO2 + O ↔ O2+OH1.5 × 10139500HO2 + OH ↔ H2O + O28.0 × 101200HO2 + HO2 ↔ H2O2+O22.0 × 101200H + H2O2 ↔ H2+HO21.4 × 101236000O + H2O2 ↔ OH + HO21.4 × 101364000OH + H2O2 ↔ H2O + HO26.1 × 101214300H2O2 + M ↔ OH + OH + M1.2 × 1017455000A structured mesh is adopted, and the height of the first layer of the mesh on the wall is 1×10−6m, which can meet y+≈1 [58, 59]. In order to ensure the precision near the wall surface of the struts, mesh division, as shown in Figure 6, is adopted, and appropriate refining is carried out behind the struts and near the nozzle.Figure 6
A local magnification of the mesh.Figure7 shows the comparison of the combustion state between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. The simulation captures the reflected shock waves on the upper and lower walls due to the separation of the boundary layer as well as the reflected shock waves on the shear layer in the combustion zone. In addition, the position and shape of the shock waves, the position of the narrowest shear layer, and the width of the turbulent vortex region are basically consistent with the experiment.Figure 7
Comparison of the combustion between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. (a) Experiment [55]. (b) Simulation.
(a)(b)In terms of the combustion of hydrogen, the OH group is an important intermediate product. The area with higher OH group concentration is basically the same as the combustion exothermic zone. Therefore, the distribution of the OH group can be used to represent the distribution of the flame. Figure8 shows the mass fraction distribution of the OH group. It can be seen that the distributions of the OH group obtained by the two models are very different. The flame obtained by the flamelet model (Figure 8(b)) is a typical attached-flame. Its starting point is at the corner point of the trailing edge of the strut, which is closely attached to the wall near the two sides of the nozzle and extends downstream along the shear layer. The concentration of the stripe area is higher, and there is almost no OH group in the range of x = 110∼160 mm behind the strut. The flame of the FGM model is lift-off flame (Figure 8(a)), and its starting point is located near x = 120 mm, which proves that the compressible modified FGM model can better capture the phenomenon of “flame lifting” in the supersonic flow field.Figure 8
Distribution of OH group. (a) Compressible modified FGM model. (b) Flamelet model.
(a)(b)Figure9 shows the comparison between the velocity profile of the two models and experiment at three locations along the x-direction. In general, the velocity distribution at three different locations obtained by the two models is not significantly different. In the center area of the flame near x = 120 mm, the airflow in the shear layer decreases to subsonic velocity due to the exothermic combustion. Outside the shear layer, due to the narrowing of the effective channel of the air inlet, the air is compressed. The simulation results are slightly lower than the experimental values in the area near the nozzle. It is in good agreement with the experiment at x = 167 mm. At x = 249 mm, the velocity profile obtained by the FGM model is consistent with the experiment, while the trend obtained by the flamelet model is opposite, which highlights the advantages of the FGM model.Figure 9
Velocity profile of the two models and experiment. (a)x = 120 mm. (b) x = 167 mm. (c) x = 249 mm.
(a)(b)(c)


### 2.3. Numerical Setup

In this paper, a three-dimensional structured mesh is established to simulate the combustion of a cavity-based scramjet combustor. In order to accurately simulate the flow and flame at the nozzle and vicinity of the cavity, the mesh is refined at the entrance of the nozzle, the vicinity of the cavity, and the shear layer of the cavity. In order to ensure the accuracy of the hybrid RANS/LES turbulence model, it is recommended that the dimensionless distance between the first layer of mesh and the wall be kept aty+≈1. Based on the theory of dimensionless distance, a mesh with the total number of approximately 6 million is generated, as shown in Figure 10. X, Y, and Z represent flow direction, longitudinal direction, and span direction, respectively. As the flame close to the lower wall is stronger and the temperature gradient is larger, the mesh needs to be secondary refined, so that the mesh on the lower wall is finer than that on the upper wall.Figure 10
Mesh of the cavity-based scramjet combustor. (a) Front view of cavity. (b) Side view of cavity. (c) Top view of nozzle.
(a)(b)(c)The LES model is superior to the RANS model for some complex flows, such as nonequilibrium, three-dimensional separated flows. However, it is still difficult to apply LES to high Reynolds number near-wall flows, because a higher resolution of mesh is required to capture the small-scale vortex structures. In contrast, RANS is more suitable for the simulation of near-wall areas, because only the mesh in the normal direction of the wall needs to be refined, and the coarse meshes can be used in the parallel direction of the wall. Combining the advantages of the RANS model and LES model, the total number of mesh can be greatly reduced by adopting the hybrid RANS/LES turbulence model.The injection conditions and inflow conditions of the fuel are consistent with the experimental conditions [21], which are given in Tables 3 and 4. The outlet adopts supersonic extrapolation boundary. The solid wall is adiabatic and nonslip boundary. The treatment near the wall is the same as the k-ω SST model. A hybrid RANS/LES turbulence model coupled with k-ω SST model proposed by Gritskevich is used to solve the turbulent flow field [60]. The combustion model adopts nonpremixed compressible modified FGM model. The reaction mechanism between hydrogen and air is based on the 19-step reaction mechanism of 9 species (H2, O2, N2, H2O, OH, H, O, HO2, H2O2) given by Jachimowski [57], as shown in Table 3. The viscosity of each species of the gas is given according to equation (6), and the specific heat is determined according to the JANAF physical properties table.Table 3
H2 inlet conditions of experiment.
Nozzle diameter (mm)Ma (-)Total pressure (MPa)Total temperature (K)Flow (g/s)Equivalence ratio (-)2.01.01.52902.8260.097Table 4
Air inlet conditions of experiment.
MaStatic temperature (K)Static pressure (kPa)YO2YH2OYN22.52753.886.20.2320.0810.687


## 2.1. Governing Equation

According to the conservation of mass, momentum, and energy, the governing equations of gas-phase flow are established. The detailed derivation of the equations is based on the work of Versteeg and Malalasekera [39], where only the final results in matrix form are given. In the actual simulation of combustion, the Navier–Stokes equation system containing multi-species and conservative chemical reactions is usually used, and the gas transport process includes three phenomena, namely: diffusion phenomenon of species transport caused by concentration gradient, viscosity phenomenon of momentum transport caused by velocity gradient, and heat conduction phenomenon of energy transport caused by temperature gradient and species diffusion. In order to calculate the momentum, species, and energy transport intensity in the combustion flow field, the transport coefficients of each species need to be given. The governing equations of the gas-phase are written as follows:(1)∂Q∂t+∂E−Ev∂x+∂F−Fv∂y+∂G−Gv∂z=H,(2)Q=ρρuρvρwρeρYiE=ρuρuu+pρuvρuwuρe+pρuYiEv=0τxxτxyτxzuτxx+vτxy+wτxz−qxρiDim∂Yi∂x,(3)F=ρvρvuρvv+pρvwvρe+pρvYiFv=0τyxτyyτyzuτxy+vτyy+wτyz−qyρiDim∂Yi∂y,(4)G=ρwρwuρwvρww+pwρe+pρwYiGv=0τzxτzyτzzuτzx+vτzy+wτzz−qzρiDim∂Yi/∂zH=00000ωi.In Equations (1)–(4): i = 1, 2, …, Ns−1, Ns is the number of gas species. Q is an independent variable that is conserved; E, F, G, Ev, Fv, and Gv are, respectively, the inviscid and viscous fluxes in the x, y, and z directions of the Cartesian coordinate system; u, v , and w are the partial velocities along the x, y, and z directions of the coordinate axes; p is the pressure of the gas-phase; ρi is the density of the species i; ρ is the density of the gas-phase mixture; Yi is the mass fraction of the species i; ωi is the mass generation rate of species i. τi,j is the component of the viscous stress, as follows:(5)τxx=−23μ∇⋅U+2μ∂u∂x,τyy=−23μ∇⋅U+2μ∂u∂y,τzz=−23μ∇⋅U+2μ∂u∂z,τxy=τyx=μ∂u∂y+∂v∂x,τyz=τzy=μ∂v∂z+∂w∂y,τxz=τzx=μ∂w∂x+∂u∂z.μ is the laminar viscosity coefficient, given by the Surferland formula:(6)μ=μ0TT01.5T0+TsT+Ts,μ0 is the reference viscosity coefficient, T0 is the reference temperature, and Ts is the characteristic temperature. For the air, they are 17.161 × 10–6 (Pa∙s), 273.16 K, and 124 K, respectively. qi (i = x, y, z) represents the energy flux caused by thermal conductivity and species diffusion, as follows:(7)qx=−λ∂T∂x−ρ∑i=1NSDimhi∂Yi∂x,qy=−λ∂T∂y−ρ∑i=1NsDimhi∂Yi∂y,qz=−λ∂T∂z−ρ∑i=1NsDimhi∂Yi∂z.Dim is the diffusion coefficient of speciesi, as follows:(8)Dim=1−Xi∑i,i≠jXj/Dij.In this paper, the thermal conductivityλ and the species diffusion coefficient Dim are given by the laminar Prandtl number Pr and Schmidt number Sc: λ = cp∙μ/Pr, Dim = μ/(ρ∙Sc).The mole fraction of the species is calculated as follows:(9)Xi=Yi/Wi∑j=1NsYj/Wj.For calorically perfect gas, when translational kinetic energy and rotational energy are fully excited, the specific heat capacity and specific heat ratio are constant, and the internal energy is only a linear function of temperature. For thermally perfect gas involving combustion, the enthalpy is no longer a linear function of temperature. The thermodynamic function of each gas species is usually given in the form of a polynomial temperature, and coefficients of the polynomial are obtained by least square method fitting. In this paper, the polynomial form given by JANAF database [40] is used to calculate the cp,i and enthalpy hi of gas at a given temperature, as follows:(10)cp,i=A+B×T+C×T2+D×T3+ET4,hi−hi,298.15=A×T+B×T22+C×T33+D×T44−ET+F−H.The enthalpy and internal energy of the gas-phase mixture are defined as follows:(11)hT=∑i=1NsYihiT,eT=∑i=1NsYihiT−pρ.In addition, the multi-species gas-phase mixture satisfies the ideal gas equation of state, and the pressure of the mixture can be obtained according to Dalton partial pressure theorem:(12)p=R0T∑i=1NsρiWi.R0 is the general gas constant, 8.3145 kJ/(kg K). Ri is the gas constant of species i. The gas constant of the gas-phase mixture is defined as follows:(13)R=∑iYiRi=∑iYiR0Wi.The specific heat capacity at constant volume of gas-phase mixture is calculated as follows:(14)cp=∑i=1NsYicpiThe molar mass of the gas-phase mixture is calculated as follows:(15)W=∑XiWi=∑YiWi−1,According to equation (11), energy (internal energy e or enthalpy h) is a nonlinear equation of temperature T. After obtaining the energy field according to the energy equation, it is necessary to solve the temperature field iteratively. Generally, Newton iteration method is adopted [41].


## 2.2. FGM Model

The FGM (flamelet-generated Manifolds) model is developed from incompressible premixed combustion [42], which is essentially a model developed based on the idea of the flamelet model [43]. The idea of the FGM model comes from ILDM (Intrinsic low-dimensional Manifolds) [44]. ILDM simplifies the chemical reaction into fast reaction and slow reaction, and only considers the manifold with the smallest time scale of chemical reaction (Manifold). Therefore, the simulation effect of fast-reaction and high-temperature section is better, while the effect of slow-reaction and low-temperature section (The effects of diffusion and convection cannot be ignored.) is poor. In order to deal with the existence of diffusion and convection, the FGM model adopts the premixed flame calculation of 1D laminar flow based on nontension in both low-temperature and high-temperature sections. In other words, low-dimensional manifolds are created by using the solution of 1D flamelet equation, and the effect of diffusion and convection is taken into account, which greatly improves the simulation of low-temperature section and high-temperature section.The FGM model can be parameterized by any form of laminar flame, which can be calculated by 1D steady premixed flamelet or diffused flamelet model. According to the different laminar flame forms, the FGM model is divided into two sub-models, premixed model and diffusion model. Generally, the premixed FGM model is used in the turbulent partially premixed flame dominated by premixed combustion. The combustion of the scramjet studied in this paper is a nonpremixed turbulent partially premixed flame, which needs to adopt the diffusion FGM model, and the required query table is obtained through the calculation of one-dimensional steady-state diffusion flame.


### 2.2.1. Establishment of Database

FGM is developed based on the flamelet model. The establishment method of its database is the same as that of the flamelet model. It needs to solve the state parameters of the laminar diffusion flame with one-dimensional structure first. This process can be realized by solving the state equations of the laminar flow collision flame (as shown in Figure1).Figure 1
Schematic of the planar laminar counter flow flame structure.Mixture fractionZ is an important concept in diffusion combustion. The mixture fraction represents the scalar of the fuel to oxidant mixing ratio, which is a function of space and time, and can be defined in many ways [45–49]. Masri gives a definition form based on the conservation of local elements [50], as follows:(16)Z=∑iyi1−yi0∑jyj1−yj0.In equation (16): yi0 and yi1 are the mass fractions of element i at the entrance of oxidant and fuel, respectively.According to the definition ofZ given in equation (16), the oxidant inlet Z = 0; the fuel inlet Z = 1. It can be guaranteed that Z is a conserved scalar without a source term. The transport equation is written as follows:(17)∂ρZ∂t+∇⋅ρUZ=∇⋅ρDZ∇Z.In equation (17): DZ is the diffusion coefficient of Z.The state of collision flame is also described by N-S equations given in Section2.1. Considering the geometric conditions and characteristics of the flame structure, equation (1) can be simplified appropriately. A new coordinate system with the normal direction of the flamelet as x is established, which is perpendicular to the flamelet and the other two directions are x1 and x2. The simplified system of equations is written as follows [48]:(18)∂ρu∂x=−ρK,(19)∂ρYi∂x=1Lei∂∂xλcp∂Yi∂x+ωi−ρKYi,i=1,…,Ns−1,(20)∂ρuh∂x=∂∂xλcp∂h∂x+∑Nshiλcp1Lei−1∂Yi∂x−ρKh.In Equations (18)–(20): x, u, and K are the space coordinate, velocity component, and flame stretch rate perpendicular to the flamelet, respectively. For the one-dimensional situation discussed in this section, K is the rate of change of u along the x direction, that is, K = ∂u/∂x.In addition, the transport equation for the mixture fractionZ is written as follows:(21)∂ρuZ∂x=∂∂xλcp∂Z∂x−ρKZ.Combined with the analysis of equations (16) and (21), Z (x) is a pure monotone function along the axis of the colliding flame (normal direction of the flame surface). Therefore, Z can be used to parameterize the flow field variables.In actual nonpremixed combustion, in order to express the state parameters in the form ofφZ in Z, the governing equation in the physical space (Figure 2) needs to be converted into the Z space (Figure 3). This process can be achieved by coordinate transformation in the type of Crocco [42].Figure 2
Variables in physical space.Figure 3
Variables in Z space.The transport equations of species and temperature are transformed as follows [51]:(22)∂Yi∂t−χ2∂2Yi∂Z2−ωiρ=0,(23)∂T∂t−χ2∂2T∂Z2+1ρcp∑i=1Nshiωi−1ρcp∂p∂t=0.In Equations (22) and (23): The dimension of Scalar Dissipation Rate χ is s−1, which characterizes the effect of chemical nonequilibrium effects. The equation of definition is written as follows:(24)χ=2D∇Z2.Equations (22) and (23) are the one-dimensional flamelet equations in Z space, and there is no convection term. After the initial and boundary conditions are determined, the distribution of chemical thermodynamic parameters in flamelet can be obtained by solving these two equations. The solution of the laminar diffusion flame represented by t, Z, and χ is defined as follows:(25)φ=φt,Z,χ.In the FGM model,χ is defined as Progress-variable C that characterizes the degree of combustion. In the case of steady-state combustion, equation (25) can be transformed into the following equation:(26)φ=φZ,C.In general, equation (26) is not directly solved. Instead, according to certain rules, equation (26) is made into a unified format database, which is known as the laminar FGM database.When the thickness of the flamelet is less than the minimum vortex scale (Kolmogorov scale) of turbulent flow (as shown in Figure4), the flamelet can still maintain a laminar state. In this case, the turbulent flame structure can be represented by an ensemble embedded in the turbulent flow field with a one-dimensional wrinkled laminar flamelet. The thermodynamic state parameters of the average flame structure in the turbulent field are given by Statistical average of the laminar flamelet, as follows:(27)φ¯=∫∫φZ,CPZ,CdZdC.Figure 4
Kolmogorov vortex and flamelet in the wrinkled laminar flamelet.In equation (27): φ¯ is the average thermodynamic state parameter required in turbulent flame; P (Z, C) is the Joint Probability Density Function (JPDF) of Z and C. Since φ has been obtained by the solution of equation (26), P (Z, C) only needs to know to calculate φ¯.It is generally assumed thatZ and C are independent of each other, so JPDF can be represented as the product of the Marginal Density-weighted Filter PDF (FPDF):(28)PZ,C=PZPC.By introducing the mean valueZ¯ and variance Z″2¯ of Z, and substituting equation (28) into equation (27), the chemical thermodynamic state of turbulent combustion can be calculated as follows:(29)φ¯=φ¯Z¯,Z″2¯,C.By integrating the calculated multi-group laminar database, the turbulent database under the condition of differentC (progress variables) can be obtained, and then the turbulent database can be made into a lookup table according to certain rules, which is known as the turbulent FGM database.The turbulent FGM database can be combined with CFD solvers to solve the combustion flow field. After the solver calculates the distribution ofZ and C, it uses the lookup table to obtain the components and other required parameters to complete the solution. Therefore, for the FGM model, in addition to solving the continuity equation, and the momentum equation and transport equation of Z, it is also necessary to solve the scalar transport equation of C:(30)∂ρC∂t+∇⋅ρUC=∇ρDC∇C+ωC.In equation (30): DC is the effective transport coefficient of the schedule variable; ωC is the chemical reaction source term of C.


### 2.2.2. Compressible Modification of FGM Model

As far as the FGM model is concerned, the thermochemical state parameters in the flamelet database include but are not limited toZ, T, ρ, Yi, and ωi. The incompressible combustion gas state parameters T, ρ, and p are generally decoupled (regardless of the buoyancy flow). T and ρ can be calculated in advance and obtained from the lookup table. However, for compressible turbulent combustion in supersonic conditions, U, T, ρ, and p are strongly coupled and limited by the ideal gas state equation. A change in kinetic energy can cause a change in temperature, and the temperature in the database cannot reflect this coupling effect. If the T and ρ in the flow field are the same as in the case of incompressibility, the parameters of the flow field will introduce a large error. In addition, the flamelet database is based on the incompressible flamelet equation. The state points represented in the database are different from the supersonic combustion state. Therefore, in order to apply the traditional FGM model to supersonic combustion, it is necessary to modify the FGM model given above.At present, the modification methods of the FGM model can be summarized into two categories: one is to improve the generation method of the database, add more compressible factors to the flame control equation or database, and strive to build a more accurate database; the other is to adjust the query mode of the database and only use part of the parameters in the database, or to make some transformation of the parameters to match the state in supersonic combustion [52, 53]. Among the numerous correction algorithms, the method proposed by Oevermann [54] in 2000 has attracted great attention due to its simplicity, feasibility, and reasonable results. This method is a shock capture algorithm, which only uses the mass fraction of species in the database without temperature or other parameters. The temperature is obtained by the energy equation. In terms of computational efficiency, this processing method only solves one more energy equation than the original flamelet model, which still retains the advantages of the flamelet model. In fact, the huge amount of calculation and the rigidity of the chemical reaction source term brought by the multi-species equation in the combustion solution are difficult to deal with. Oevermann solved these two problems skillfully with the help of the flamelet model and the database query algorithm.Based on the above analysis, the flamelet model modification algorithm proposed by Oevermann was adopted to generate the compressible modified FGM model and conduct a numerical study on supersonic combustion. The specific implementation process is as follows:(1)
To solve the energy distributionh(T) in the flow field from the energy transport equation;(2)
To query the composition distributionYi (T) in the flow field through the database;(3)
To substituteh(T), Yi (T) and the thermodynamic parameters cp,i (T) of each species (equation (10)) into Equation (11). The temperature field T is obtained by the Newton iteration method.Compressibility of gas should be taken into full consideration when simulating the flow field supersonic combustion. With the compressible modification of the FGM model, it is not necessary to query the temperatureT from the FGM database when solving the temperature field, and the database only provides the parameters of the species. At the same time, a step of solving the energy equation is added. By substituting the solution of the energy equation and the species parameters obtained from database query into Equation (11) for iteration, the error caused by directly querying the database to obtain temperature T can be eliminated. Compressible modified FGM model can achieve accurate simulation of supersonic combustion while retaining the advantages of simulating the flame form.


### 2.2.3. Verification of Compressible Modified FGM Model

In order to verify the stability and integrity of the model, the traditional flamelet model and the modified FGM model are used to numerically simulate the combustion flow field of the typical DLR scramjet combustion chamber. The structure of the combustion chamber is shown in Figure5 [55]. The preheated air expands into the combustion chamber through the nozzle. The conditions of air inflow and hydrogen injection remain the same as the experimental ones, as shown in Table 1. The combustion chamber has a three-dimensional structure and relatively dense nozzle holes are arranged in the span direction. After comparing the two-dimensional and three-dimensional simulations in detail, Luo et al. [56] found that in addition to the three-dimensional effect near the nozzle, the main flow field presents obvious two-dimensional characteristics. Therefore, the RANS method is adopted in this section to simulate the two-dimensional simplified structure, and the nozzle is equivalent to a slit in the same area. For the reaction mechanism of hydrogen and air used to generate laminar flamelet database, a 19-step reaction mechanism of 9 species (H2, O2, N2, H2O, OH, H, O, HO2, H2O2) proposed by Jachimowski [57] is adopted, as shown in Table 2.Figure 5
Schematic diagrams of the DLR scramjet combustion chamber.Table 1
Experiment conditions of the DLR scramjet combustion chamber.
AirH2Ma2.01.0U/(m/s)7301200T/K340250p/MPa0.10.1ρ/(kg/m3)1.0020.097YO20.2320YN20.7360YH2O0.0320YH201k/(m2/s2)102400ε/(1/s)722462963Table 2
H2-O2-N2 9-specie 19-step reaction mechanism [57].
ReactionA(mole/(s∙cm3))E(cal/mol)B (-)H2 + O2 ↔ HO2+H1.0 × 1014560000H + O2 ↔ OH + O2.6 × 1014168000O + H2 ↔ OH + H1.8 × 101089001OH + H2 ↔ H + H2O2.2 × 101351500OH + OH ↔ O + H2O6.3 × 101210900H + OH + M ↔ H2O + M2.2 × 10220−2H + H + M ↔ H2 + M6.4 × 10170−1H + O + M ↔ OH + M6.0 × 10160−0.6H + O2 + M ↔ HO2 + M2.1 × 1015−10000O + O + M ↔ O2 + M6.0 × 1017−18000HO2 + H ↔ OH + OH1.4 × 101410800HO2 + H ↔ H2O + O1.0 × 101310800HO2 + O ↔ O2+OH1.5 × 10139500HO2 + OH ↔ H2O + O28.0 × 101200HO2 + HO2 ↔ H2O2+O22.0 × 101200H + H2O2 ↔ H2+HO21.4 × 101236000O + H2O2 ↔ OH + HO21.4 × 101364000OH + H2O2 ↔ H2O + HO26.1 × 101214300H2O2 + M ↔ OH + OH + M1.2 × 1017455000A structured mesh is adopted, and the height of the first layer of the mesh on the wall is 1×10−6m, which can meet y+≈1 [58, 59]. In order to ensure the precision near the wall surface of the struts, mesh division, as shown in Figure 6, is adopted, and appropriate refining is carried out behind the struts and near the nozzle.Figure 6
A local magnification of the mesh.Figure7 shows the comparison of the combustion state between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. The simulation captures the reflected shock waves on the upper and lower walls due to the separation of the boundary layer as well as the reflected shock waves on the shear layer in the combustion zone. In addition, the position and shape of the shock waves, the position of the narrowest shear layer, and the width of the turbulent vortex region are basically consistent with the experiment.Figure 7
Comparison of the combustion between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. (a) Experiment [55]. (b) Simulation.
(a)(b)In terms of the combustion of hydrogen, the OH group is an important intermediate product. The area with higher OH group concentration is basically the same as the combustion exothermic zone. Therefore, the distribution of the OH group can be used to represent the distribution of the flame. Figure8 shows the mass fraction distribution of the OH group. It can be seen that the distributions of the OH group obtained by the two models are very different. The flame obtained by the flamelet model (Figure 8(b)) is a typical attached-flame. Its starting point is at the corner point of the trailing edge of the strut, which is closely attached to the wall near the two sides of the nozzle and extends downstream along the shear layer. The concentration of the stripe area is higher, and there is almost no OH group in the range of x = 110∼160 mm behind the strut. The flame of the FGM model is lift-off flame (Figure 8(a)), and its starting point is located near x = 120 mm, which proves that the compressible modified FGM model can better capture the phenomenon of “flame lifting” in the supersonic flow field.Figure 8
Distribution of OH group. (a) Compressible modified FGM model. (b) Flamelet model.
(a)(b)Figure9 shows the comparison between the velocity profile of the two models and experiment at three locations along the x-direction. In general, the velocity distribution at three different locations obtained by the two models is not significantly different. In the center area of the flame near x = 120 mm, the airflow in the shear layer decreases to subsonic velocity due to the exothermic combustion. Outside the shear layer, due to the narrowing of the effective channel of the air inlet, the air is compressed. The simulation results are slightly lower than the experimental values in the area near the nozzle. It is in good agreement with the experiment at x = 167 mm. At x = 249 mm, the velocity profile obtained by the FGM model is consistent with the experiment, while the trend obtained by the flamelet model is opposite, which highlights the advantages of the FGM model.Figure 9
Velocity profile of the two models and experiment. (a)x = 120 mm. (b) x = 167 mm. (c) x = 249 mm.
(a)(b)(c)


## 2.2.1. Establishment of Database

FGM is developed based on the flamelet model. The establishment method of its database is the same as that of the flamelet model. It needs to solve the state parameters of the laminar diffusion flame with one-dimensional structure first. This process can be realized by solving the state equations of the laminar flow collision flame (as shown in Figure1).Figure 1
Schematic of the planar laminar counter flow flame structure.Mixture fractionZ is an important concept in diffusion combustion. The mixture fraction represents the scalar of the fuel to oxidant mixing ratio, which is a function of space and time, and can be defined in many ways [45–49]. Masri gives a definition form based on the conservation of local elements [50], as follows:(16)Z=∑iyi1−yi0∑jyj1−yj0.In equation (16): yi0 and yi1 are the mass fractions of element i at the entrance of oxidant and fuel, respectively.According to the definition ofZ given in equation (16), the oxidant inlet Z = 0; the fuel inlet Z = 1. It can be guaranteed that Z is a conserved scalar without a source term. The transport equation is written as follows:(17)∂ρZ∂t+∇⋅ρUZ=∇⋅ρDZ∇Z.In equation (17): DZ is the diffusion coefficient of Z.The state of collision flame is also described by N-S equations given in Section2.1. Considering the geometric conditions and characteristics of the flame structure, equation (1) can be simplified appropriately. A new coordinate system with the normal direction of the flamelet as x is established, which is perpendicular to the flamelet and the other two directions are x1 and x2. The simplified system of equations is written as follows [48]:(18)∂ρu∂x=−ρK,(19)∂ρYi∂x=1Lei∂∂xλcp∂Yi∂x+ωi−ρKYi,i=1,…,Ns−1,(20)∂ρuh∂x=∂∂xλcp∂h∂x+∑Nshiλcp1Lei−1∂Yi∂x−ρKh.In Equations (18)–(20): x, u, and K are the space coordinate, velocity component, and flame stretch rate perpendicular to the flamelet, respectively. For the one-dimensional situation discussed in this section, K is the rate of change of u along the x direction, that is, K = ∂u/∂x.In addition, the transport equation for the mixture fractionZ is written as follows:(21)∂ρuZ∂x=∂∂xλcp∂Z∂x−ρKZ.Combined with the analysis of equations (16) and (21), Z (x) is a pure monotone function along the axis of the colliding flame (normal direction of the flame surface). Therefore, Z can be used to parameterize the flow field variables.In actual nonpremixed combustion, in order to express the state parameters in the form ofφZ in Z, the governing equation in the physical space (Figure 2) needs to be converted into the Z space (Figure 3). This process can be achieved by coordinate transformation in the type of Crocco [42].Figure 2
Variables in physical space.Figure 3
Variables in Z space.The transport equations of species and temperature are transformed as follows [51]:(22)∂Yi∂t−χ2∂2Yi∂Z2−ωiρ=0,(23)∂T∂t−χ2∂2T∂Z2+1ρcp∑i=1Nshiωi−1ρcp∂p∂t=0.In Equations (22) and (23): The dimension of Scalar Dissipation Rate χ is s−1, which characterizes the effect of chemical nonequilibrium effects. The equation of definition is written as follows:(24)χ=2D∇Z2.Equations (22) and (23) are the one-dimensional flamelet equations in Z space, and there is no convection term. After the initial and boundary conditions are determined, the distribution of chemical thermodynamic parameters in flamelet can be obtained by solving these two equations. The solution of the laminar diffusion flame represented by t, Z, and χ is defined as follows:(25)φ=φt,Z,χ.In the FGM model,χ is defined as Progress-variable C that characterizes the degree of combustion. In the case of steady-state combustion, equation (25) can be transformed into the following equation:(26)φ=φZ,C.In general, equation (26) is not directly solved. Instead, according to certain rules, equation (26) is made into a unified format database, which is known as the laminar FGM database.When the thickness of the flamelet is less than the minimum vortex scale (Kolmogorov scale) of turbulent flow (as shown in Figure4), the flamelet can still maintain a laminar state. In this case, the turbulent flame structure can be represented by an ensemble embedded in the turbulent flow field with a one-dimensional wrinkled laminar flamelet. The thermodynamic state parameters of the average flame structure in the turbulent field are given by Statistical average of the laminar flamelet, as follows:(27)φ¯=∫∫φZ,CPZ,CdZdC.Figure 4
Kolmogorov vortex and flamelet in the wrinkled laminar flamelet.In equation (27): φ¯ is the average thermodynamic state parameter required in turbulent flame; P (Z, C) is the Joint Probability Density Function (JPDF) of Z and C. Since φ has been obtained by the solution of equation (26), P (Z, C) only needs to know to calculate φ¯.It is generally assumed thatZ and C are independent of each other, so JPDF can be represented as the product of the Marginal Density-weighted Filter PDF (FPDF):(28)PZ,C=PZPC.By introducing the mean valueZ¯ and variance Z″2¯ of Z, and substituting equation (28) into equation (27), the chemical thermodynamic state of turbulent combustion can be calculated as follows:(29)φ¯=φ¯Z¯,Z″2¯,C.By integrating the calculated multi-group laminar database, the turbulent database under the condition of differentC (progress variables) can be obtained, and then the turbulent database can be made into a lookup table according to certain rules, which is known as the turbulent FGM database.The turbulent FGM database can be combined with CFD solvers to solve the combustion flow field. After the solver calculates the distribution ofZ and C, it uses the lookup table to obtain the components and other required parameters to complete the solution. Therefore, for the FGM model, in addition to solving the continuity equation, and the momentum equation and transport equation of Z, it is also necessary to solve the scalar transport equation of C:(30)∂ρC∂t+∇⋅ρUC=∇ρDC∇C+ωC.In equation (30): DC is the effective transport coefficient of the schedule variable; ωC is the chemical reaction source term of C.


## 2.2.2. Compressible Modification of FGM Model

As far as the FGM model is concerned, the thermochemical state parameters in the flamelet database include but are not limited toZ, T, ρ, Yi, and ωi. The incompressible combustion gas state parameters T, ρ, and p are generally decoupled (regardless of the buoyancy flow). T and ρ can be calculated in advance and obtained from the lookup table. However, for compressible turbulent combustion in supersonic conditions, U, T, ρ, and p are strongly coupled and limited by the ideal gas state equation. A change in kinetic energy can cause a change in temperature, and the temperature in the database cannot reflect this coupling effect. If the T and ρ in the flow field are the same as in the case of incompressibility, the parameters of the flow field will introduce a large error. In addition, the flamelet database is based on the incompressible flamelet equation. The state points represented in the database are different from the supersonic combustion state. Therefore, in order to apply the traditional FGM model to supersonic combustion, it is necessary to modify the FGM model given above.At present, the modification methods of the FGM model can be summarized into two categories: one is to improve the generation method of the database, add more compressible factors to the flame control equation or database, and strive to build a more accurate database; the other is to adjust the query mode of the database and only use part of the parameters in the database, or to make some transformation of the parameters to match the state in supersonic combustion [52, 53]. Among the numerous correction algorithms, the method proposed by Oevermann [54] in 2000 has attracted great attention due to its simplicity, feasibility, and reasonable results. This method is a shock capture algorithm, which only uses the mass fraction of species in the database without temperature or other parameters. The temperature is obtained by the energy equation. In terms of computational efficiency, this processing method only solves one more energy equation than the original flamelet model, which still retains the advantages of the flamelet model. In fact, the huge amount of calculation and the rigidity of the chemical reaction source term brought by the multi-species equation in the combustion solution are difficult to deal with. Oevermann solved these two problems skillfully with the help of the flamelet model and the database query algorithm.Based on the above analysis, the flamelet model modification algorithm proposed by Oevermann was adopted to generate the compressible modified FGM model and conduct a numerical study on supersonic combustion. The specific implementation process is as follows:(1)
To solve the energy distributionh(T) in the flow field from the energy transport equation;(2)
To query the composition distributionYi (T) in the flow field through the database;(3)
To substituteh(T), Yi (T) and the thermodynamic parameters cp,i (T) of each species (equation (10)) into Equation (11). The temperature field T is obtained by the Newton iteration method.Compressibility of gas should be taken into full consideration when simulating the flow field supersonic combustion. With the compressible modification of the FGM model, it is not necessary to query the temperatureT from the FGM database when solving the temperature field, and the database only provides the parameters of the species. At the same time, a step of solving the energy equation is added. By substituting the solution of the energy equation and the species parameters obtained from database query into Equation (11) for iteration, the error caused by directly querying the database to obtain temperature T can be eliminated. Compressible modified FGM model can achieve accurate simulation of supersonic combustion while retaining the advantages of simulating the flame form.


## 2.2.3. Verification of Compressible Modified FGM Model

In order to verify the stability and integrity of the model, the traditional flamelet model and the modified FGM model are used to numerically simulate the combustion flow field of the typical DLR scramjet combustion chamber. The structure of the combustion chamber is shown in Figure5 [55]. The preheated air expands into the combustion chamber through the nozzle. The conditions of air inflow and hydrogen injection remain the same as the experimental ones, as shown in Table 1. The combustion chamber has a three-dimensional structure and relatively dense nozzle holes are arranged in the span direction. After comparing the two-dimensional and three-dimensional simulations in detail, Luo et al. [56] found that in addition to the three-dimensional effect near the nozzle, the main flow field presents obvious two-dimensional characteristics. Therefore, the RANS method is adopted in this section to simulate the two-dimensional simplified structure, and the nozzle is equivalent to a slit in the same area. For the reaction mechanism of hydrogen and air used to generate laminar flamelet database, a 19-step reaction mechanism of 9 species (H2, O2, N2, H2O, OH, H, O, HO2, H2O2) proposed by Jachimowski [57] is adopted, as shown in Table 2.Figure 5
Schematic diagrams of the DLR scramjet combustion chamber.Table 1
Experiment conditions of the DLR scramjet combustion chamber.
AirH2Ma2.01.0U/(m/s)7301200T/K340250p/MPa0.10.1ρ/(kg/m3)1.0020.097YO20.2320YN20.7360YH2O0.0320YH201k/(m2/s2)102400ε/(1/s)722462963Table 2
H2-O2-N2 9-specie 19-step reaction mechanism [57].
ReactionA(mole/(s∙cm3))E(cal/mol)B (-)H2 + O2 ↔ HO2+H1.0 × 1014560000H + O2 ↔ OH + O2.6 × 1014168000O + H2 ↔ OH + H1.8 × 101089001OH + H2 ↔ H + H2O2.2 × 101351500OH + OH ↔ O + H2O6.3 × 101210900H + OH + M ↔ H2O + M2.2 × 10220−2H + H + M ↔ H2 + M6.4 × 10170−1H + O + M ↔ OH + M6.0 × 10160−0.6H + O2 + M ↔ HO2 + M2.1 × 1015−10000O + O + M ↔ O2 + M6.0 × 1017−18000HO2 + H ↔ OH + OH1.4 × 101410800HO2 + H ↔ H2O + O1.0 × 101310800HO2 + O ↔ O2+OH1.5 × 10139500HO2 + OH ↔ H2O + O28.0 × 101200HO2 + HO2 ↔ H2O2+O22.0 × 101200H + H2O2 ↔ H2+HO21.4 × 101236000O + H2O2 ↔ OH + HO21.4 × 101364000OH + H2O2 ↔ H2O + HO26.1 × 101214300H2O2 + M ↔ OH + OH + M1.2 × 1017455000A structured mesh is adopted, and the height of the first layer of the mesh on the wall is 1×10−6m, which can meet y+≈1 [58, 59]. In order to ensure the precision near the wall surface of the struts, mesh division, as shown in Figure 6, is adopted, and appropriate refining is carried out behind the struts and near the nozzle.Figure 6
A local magnification of the mesh.Figure7 shows the comparison of the combustion state between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. The simulation captures the reflected shock waves on the upper and lower walls due to the separation of the boundary layer as well as the reflected shock waves on the shear layer in the combustion zone. In addition, the position and shape of the shock waves, the position of the narrowest shear layer, and the width of the turbulent vortex region are basically consistent with the experiment.Figure 7
Comparison of the combustion between the experimental schlieren diagram and the time-averaged numerical schlieren diagram. (a) Experiment [55]. (b) Simulation.
(a)(b)In terms of the combustion of hydrogen, the OH group is an important intermediate product. The area with higher OH group concentration is basically the same as the combustion exothermic zone. Therefore, the distribution of the OH group can be used to represent the distribution of the flame. Figure8 shows the mass fraction distribution of the OH group. It can be seen that the distributions of the OH group obtained by the two models are very different. The flame obtained by the flamelet model (Figure 8(b)) is a typical attached-flame. Its starting point is at the corner point of the trailing edge of the strut, which is closely attached to the wall near the two sides of the nozzle and extends downstream along the shear layer. The concentration of the stripe area is higher, and there is almost no OH group in the range of x = 110∼160 mm behind the strut. The flame of the FGM model is lift-off flame (Figure 8(a)), and its starting point is located near x = 120 mm, which proves that the compressible modified FGM model can better capture the phenomenon of “flame lifting” in the supersonic flow field.Figure 8
Distribution of OH group. (a) Compressible modified FGM model. (b) Flamelet model.
(a)(b)Figure9 shows the comparison between the velocity profile of the two models and experiment at three locations along the x-direction. In general, the velocity distribution at three different locations obtained by the two models is not significantly different. In the center area of the flame near x = 120 mm, the airflow in the shear layer decreases to subsonic velocity due to the exothermic combustion. Outside the shear layer, due to the narrowing of the effective channel of the air inlet, the air is compressed. The simulation results are slightly lower than the experimental values in the area near the nozzle. It is in good agreement with the experiment at x = 167 mm. At x = 249 mm, the velocity profile obtained by the FGM model is consistent with the experiment, while the trend obtained by the flamelet model is opposite, which highlights the advantages of the FGM model.Figure 9
Velocity profile of the two models and experiment. (a)x = 120 mm. (b) x = 167 mm. (c) x = 249 mm.
(a)(b)(c)


## 2.3. Numerical Setup

In this paper, a three-dimensional structured mesh is established to simulate the combustion of a cavity-based scramjet combustor. In order to accurately simulate the flow and flame at the nozzle and vicinity of the cavity, the mesh is refined at the entrance of the nozzle, the vicinity of the cavity, and the shear layer of the cavity. In order to ensure the accuracy of the hybrid RANS/LES turbulence model, it is recommended that the dimensionless distance between the first layer of mesh and the wall be kept aty+≈1. Based on the theory of dimensionless distance, a mesh with the total number of approximately 6 million is generated, as shown in Figure 10. X, Y, and Z represent flow direction, longitudinal direction, and span direction, respectively. As the flame close to the lower wall is stronger and the temperature gradient is larger, the mesh needs to be secondary refined, so that the mesh on the lower wall is finer than that on the upper wall.Figure 10
Mesh of the cavity-based scramjet combustor. (a) Front view of cavity. (b) Side view of cavity. (c) Top view of nozzle.
(a)(b)(c)The LES model is superior to the RANS model for some complex flows, such as nonequilibrium, three-dimensional separated flows. However, it is still difficult to apply LES to high Reynolds number near-wall flows, because a higher resolution of mesh is required to capture the small-scale vortex structures. In contrast, RANS is more suitable for the simulation of near-wall areas, because only the mesh in the normal direction of the wall needs to be refined, and the coarse meshes can be used in the parallel direction of the wall. Combining the advantages of the RANS model and LES model, the total number of mesh can be greatly reduced by adopting the hybrid RANS/LES turbulence model.The injection conditions and inflow conditions of the fuel are consistent with the experimental conditions [21], which are given in Tables 3 and 4. The outlet adopts supersonic extrapolation boundary. The solid wall is adiabatic and nonslip boundary. The treatment near the wall is the same as the k-ω SST model. A hybrid RANS/LES turbulence model coupled with k-ω SST model proposed by Gritskevich is used to solve the turbulent flow field [60]. The combustion model adopts nonpremixed compressible modified FGM model. The reaction mechanism between hydrogen and air is based on the 19-step reaction mechanism of 9 species (H2, O2, N2, H2O, OH, H, O, HO2, H2O2) given by Jachimowski [57], as shown in Table 3. The viscosity of each species of the gas is given according to equation (6), and the specific heat is determined according to the JANAF physical properties table.Table 3
H2 inlet conditions of experiment.
Nozzle diameter (mm)Ma (-)Total pressure (MPa)Total temperature (K)Flow (g/s)Equivalence ratio (-)2.01.01.52902.8260.097Table 4
Air inlet conditions of experiment.
MaStatic temperature (K)Static pressure (kPa)YO2YH2OYN22.52753.886.20.2320.0810.687


## 3. Numerical Simulation Results

### 3.1. Three-Dimensional Transient Temperature Distribution

In this paper, a compressible modified FGM model is used to simulate the combustion process in a cavity-based scramjet combustor. The results of the simulation at 0.7 ms (combustion is fully developed) are shown byY-Z direction sections (S1 to S7) and middle section. Figure 11 shows the temperature distribution of S1 to S7. There is an obvious high-temperature zone in the cavity and downstream, with a temperature of more than 2300 K, which proves that the combustion reaction between hydrogen and high Mach number incoming air has successfully occurred. In general, the distribution of temperature field can indirectly reflect the flame distribution in the combustion chamber. The high-temperature zone in the figure is mainly distributed in the middle and lower regions of the combustion chamber. In the Z-direction, the flame diffuses to the side of the combustion chamber and has been fully developed, but the flame height near the wall is relatively low. It can be seen that the flame peak is mainly formed at S4 and S5. Under the effect of turbulence, the flame peak is unstable and there is local disturbance, which shows that the flame entrainment is formed near the flame peak. There is also a high-temperature zone in the middle and downstream of the upper wall, but the temperature is significantly lower than that of the flame in main-flow, and the temperature is about 1400–1600 K. In addition, it is observed from S2 that there is a small-scale high-temperature zone near the wall upstream of the hydrogen nozzle, indicating that a weak flame occurs in this zone. The flame height is low without obvious flame peaks. The weak flame only forms near the nozzle and does not diffuse to S1, indicating that it has no effect on the inlet of the combustion chamber.Figure 11
Temperature field of cavity-based scramjet combustor (S1-S7).Figure12 shows the streamline distribution and T = 2300 K temperature iso-surface based on the flow field in the middle section to demonstrate the diffusion process of combustion. Combustion diffuses from the hydrogen nozzle to the cavity and downstream, and presents an obvious trend of turbulent combustion. Part of the incoming air is mixed with the injected hydrogen and combusts, forming a turbulent flame. The other part is not involved in the combustion and is lifted to the top of the combustion chamber under the coupling effect of shock waves and shear layer. This process does not produce “blockage” in the combustion chamber. The detailed combustion evolution process will be analyzed and discussed in Section 3.2–3.3. Figures 11 and 12 show that the compressible modified FGM model can successfully simulate the three-dimensional turbulent combustion in the cavity-based scramjet combustor.Figure 12
Streamline andT = 2300 K temperature iso-surface.


### 3.2. Simulation Verification

Figure13 shows the results of simulation verification. Figure 13(a) shows the high-speed photography of three-dimensional combustion flow field in reference [21], which is obtained by superposition of flames in Z-direction within a certain exposure time. Figure 13(b) shows the simulation result of temperature for easy comparison, which is the mean temperature distribution of the middle section in 0.7 ms. The distribution of the high-temperature zone and the protrusion structure on the flame surface are consistent, which verifies the reliability of the algorithm in this paper. Figure 14 shows the transient simulation results at 0.7 ms. As can be seen from Figure 14(a), the high-temperature zone mainly occurs in the cavity and downstream. Figure 14(b) is the Mach number distribution in the flow field at 0.7 ms, which is superimposed with schlieren effect. In other words, Z-direction coordinates are set as density and Laplace transform processing is performed. The “Oblique compression wave,” “Arc shock wave,” and “Reflected shock wave” can be clearly observed in the figure. The Mach number of the inlet section of the combustion chamber is about 2.52. When gas flows through the inlet section of the combustion chamber, oblique compression waves are generated on the wall. The oblique compression waves belong to the weak wave, which are different from shock waves, and have little impact on downstream combustion, rather than the false shock wave caused by numerical error. The main reason for this wave is that the coupling effect between the supersonic inflow and the no-slip wall leads to a great change in the velocity gradient near the wall, thus generating the oblique compression wave caused by the inlet effect, which is inevitable, and its impact on the combustion in the cavity can be ignored. When passing through the cavity, the local Mach number is reduced to less than 1.00, and a stable flame structure is formed in the cavity. The cavity structure plays a stabilizing role in the combustion.Figure 13
Simulation verification. (a) High-speed photograpy of the combustion flow field [21]. (b) Mean temperature distribution.
(a)(b)Figure 14
Simulation results in the middle section. (a) Temperature distribution. (b) Mach number distribution. (A: Oblique compression wave; B: Arc shock wave; C: Reflected shock wave).
(a)(b)Combined with the analyses in Figure14, the arc shock wave generated by the injection on the lower wall hits the upper wall, and the adverse pressure gradient caused by the high pressure behind the shock wave induces the separation of the boundary layer on the upper wall and forms the reflected shock wave. The reflected shock wave is superimposed with the arc shock wave, and the pressure and temperature behind shock waves are correspondingly increased, which act on the flame structure on the upper cavity: On the one hand, the static temperature and pressure of the air are increased, which is beneficial for ignition and maintaining combustion; on the other hand, the superimposed shock waves act on the shear layer formed by fuel and air, which promotes the transformation of the shear layer to the turbulent flow. These facilitate the mixing of fuel and air and promote the combustion. With the effect of shock waves and combustion, a flame protrusion structure with high temperature and high pressure is formed in the intersection area between the reflected shock wave on the upper wall and the fuel shear layer.


### 3.3. Transient Flame in the Middle Section

Figure15 shows the temperature distribution of the flow field in the middle section at different moments when the combustion is fully developed. In order to facilitate the analysis of the flame propagation mechanism, multiple groups of temperature iso-surface are superimposed and treated transparently. The flame propagation process is shown as follows:(1)
Due to the low velocity in the cavity, stable combustion is first formed in the cavity, as shown in Figure15(a);(2)
As time goes by, the high-pressure environment formed by combustion pushes the high-temperature gas away from the cavity and gradually into the main flow, while the high-temperature gas adheres to the wall and gradually spreads downstream, as shown in Figure15(b);(3)
The high temperature zone starts from the back wall of the cavity and gradually reaches upstream to occupy the main-flow, as shown in Figure15(c);(4)
The phenomenon of “flame lifting” is observed in Figure15(d). A stable flame structure is formed clearly in the main-flow. The high temperature zone has also expanded from about 1/3 of the channel height shown in Figure 15(a) to about 2/3.Figure 15
Simulation results of transient temperature. (a) 0.5 ms. (b) 0.6 ms. (c) 0.7 ms. (d) 0.8 ms.
(a)(b)(c)(d)


## 3.1. Three-Dimensional Transient Temperature Distribution

In this paper, a compressible modified FGM model is used to simulate the combustion process in a cavity-based scramjet combustor. The results of the simulation at 0.7 ms (combustion is fully developed) are shown byY-Z direction sections (S1 to S7) and middle section. Figure 11 shows the temperature distribution of S1 to S7. There is an obvious high-temperature zone in the cavity and downstream, with a temperature of more than 2300 K, which proves that the combustion reaction between hydrogen and high Mach number incoming air has successfully occurred. In general, the distribution of temperature field can indirectly reflect the flame distribution in the combustion chamber. The high-temperature zone in the figure is mainly distributed in the middle and lower regions of the combustion chamber. In the Z-direction, the flame diffuses to the side of the combustion chamber and has been fully developed, but the flame height near the wall is relatively low. It can be seen that the flame peak is mainly formed at S4 and S5. Under the effect of turbulence, the flame peak is unstable and there is local disturbance, which shows that the flame entrainment is formed near the flame peak. There is also a high-temperature zone in the middle and downstream of the upper wall, but the temperature is significantly lower than that of the flame in main-flow, and the temperature is about 1400–1600 K. In addition, it is observed from S2 that there is a small-scale high-temperature zone near the wall upstream of the hydrogen nozzle, indicating that a weak flame occurs in this zone. The flame height is low without obvious flame peaks. The weak flame only forms near the nozzle and does not diffuse to S1, indicating that it has no effect on the inlet of the combustion chamber.Figure 11
Temperature field of cavity-based scramjet combustor (S1-S7).Figure12 shows the streamline distribution and T = 2300 K temperature iso-surface based on the flow field in the middle section to demonstrate the diffusion process of combustion. Combustion diffuses from the hydrogen nozzle to the cavity and downstream, and presents an obvious trend of turbulent combustion. Part of the incoming air is mixed with the injected hydrogen and combusts, forming a turbulent flame. The other part is not involved in the combustion and is lifted to the top of the combustion chamber under the coupling effect of shock waves and shear layer. This process does not produce “blockage” in the combustion chamber. The detailed combustion evolution process will be analyzed and discussed in Section 3.2–3.3. Figures 11 and 12 show that the compressible modified FGM model can successfully simulate the three-dimensional turbulent combustion in the cavity-based scramjet combustor.Figure 12
Streamline andT = 2300 K temperature iso-surface.


## 3.2. Simulation Verification

Figure13 shows the results of simulation verification. Figure 13(a) shows the high-speed photography of three-dimensional combustion flow field in reference [21], which is obtained by superposition of flames in Z-direction within a certain exposure time. Figure 13(b) shows the simulation result of temperature for easy comparison, which is the mean temperature distribution of the middle section in 0.7 ms. The distribution of the high-temperature zone and the protrusion structure on the flame surface are consistent, which verifies the reliability of the algorithm in this paper. Figure 14 shows the transient simulation results at 0.7 ms. As can be seen from Figure 14(a), the high-temperature zone mainly occurs in the cavity and downstream. Figure 14(b) is the Mach number distribution in the flow field at 0.7 ms, which is superimposed with schlieren effect. In other words, Z-direction coordinates are set as density and Laplace transform processing is performed. The “Oblique compression wave,” “Arc shock wave,” and “Reflected shock wave” can be clearly observed in the figure. The Mach number of the inlet section of the combustion chamber is about 2.52. When gas flows through the inlet section of the combustion chamber, oblique compression waves are generated on the wall. The oblique compression waves belong to the weak wave, which are different from shock waves, and have little impact on downstream combustion, rather than the false shock wave caused by numerical error. The main reason for this wave is that the coupling effect between the supersonic inflow and the no-slip wall leads to a great change in the velocity gradient near the wall, thus generating the oblique compression wave caused by the inlet effect, which is inevitable, and its impact on the combustion in the cavity can be ignored. When passing through the cavity, the local Mach number is reduced to less than 1.00, and a stable flame structure is formed in the cavity. The cavity structure plays a stabilizing role in the combustion.Figure 13
Simulation verification. (a) High-speed photograpy of the combustion flow field [21]. (b) Mean temperature distribution.
(a)(b)Figure 14
Simulation results in the middle section. (a) Temperature distribution. (b) Mach number distribution. (A: Oblique compression wave; B: Arc shock wave; C: Reflected shock wave).
(a)(b)Combined with the analyses in Figure14, the arc shock wave generated by the injection on the lower wall hits the upper wall, and the adverse pressure gradient caused by the high pressure behind the shock wave induces the separation of the boundary layer on the upper wall and forms the reflected shock wave. The reflected shock wave is superimposed with the arc shock wave, and the pressure and temperature behind shock waves are correspondingly increased, which act on the flame structure on the upper cavity: On the one hand, the static temperature and pressure of the air are increased, which is beneficial for ignition and maintaining combustion; on the other hand, the superimposed shock waves act on the shear layer formed by fuel and air, which promotes the transformation of the shear layer to the turbulent flow. These facilitate the mixing of fuel and air and promote the combustion. With the effect of shock waves and combustion, a flame protrusion structure with high temperature and high pressure is formed in the intersection area between the reflected shock wave on the upper wall and the fuel shear layer.


## 3.3. Transient Flame in the Middle Section

Figure15 shows the temperature distribution of the flow field in the middle section at different moments when the combustion is fully developed. In order to facilitate the analysis of the flame propagation mechanism, multiple groups of temperature iso-surface are superimposed and treated transparently. The flame propagation process is shown as follows:(1)
Due to the low velocity in the cavity, stable combustion is first formed in the cavity, as shown in Figure15(a);(2)
As time goes by, the high-pressure environment formed by combustion pushes the high-temperature gas away from the cavity and gradually into the main flow, while the high-temperature gas adheres to the wall and gradually spreads downstream, as shown in Figure15(b);(3)
The high temperature zone starts from the back wall of the cavity and gradually reaches upstream to occupy the main-flow, as shown in Figure15(c);(4)
The phenomenon of “flame lifting” is observed in Figure15(d). A stable flame structure is formed clearly in the main-flow. The high temperature zone has also expanded from about 1/3 of the channel height shown in Figure 15(a) to about 2/3.Figure 15
Simulation results of transient temperature. (a) 0.5 ms. (b) 0.6 ms. (c) 0.7 ms. (d) 0.8 ms.
(a)(b)(c)(d)


## 4. Analyses and Discussions

### 4.1. Transient Flame Structure Based on the OH-PLIF Image

In order to study the transient flame structure in the flow field of combustion, the simulated concentration distribution diagram of the OH group was compared with the OH-PLIF image of the experiment in reference [21]. Figures 16 and 17 show the distribution of the OH group in the middle section obtained by experiment and simulation, respectively. The analysis shows that the flame distribution of the lateral injection on the upstream wall of the cavity has the following characteristics: The flame is mainly located in the upper part of the cavity and the mainflow near the rear edge of the cavity. The simulation also captures this phenomenon and finds the distribution of the OH group near the nozzle. It means that part of the fuel is entrained into the reflux zone upstream of the nozzle and combusted, which is different from the experiment. At the same time, the simulation observes that the combustion inside the cavity is relatively intense, while the experimental observation shows that the combustion on the lower wall of the cavity is weak, for two reasons: Due to the weak spontaneous emission signal of the OH group in the flame, a long exposure time (50 μs in the experiment) is needed to obtain the spontaneous emission image with high signal-to-noise ratio. To sum up, the time-averaged flame image is obtained by spontaneous emission, rather than the transient flame image. In addition, the OH-PLIF image is the projection of a three-dimensional entity on a two-dimensional plane, so it cannot detect the information on a single section. However, the simulation results in Figure 17 are obtained from a two-dimensional slice of the middle section, which is a transient flame that cannot be observed through experiments. To sum up, there are certain differences between the experimental and simulation results. In addition, Fan et al. also emphasized in their work that the combustion has an unsteady effect, which will make it difficult to control the measurement precision of the OH-PLIF technology, and the measurement results may have certain errors.Figure 16
OH-PLIF image of the experiment [21].Figure 17
OH mass concentration of the simulation.


### 4.2. Analysis of the Injection Process Upstream of the Cavity

Combined with the schematic diagram of the shock waves and the vortex structure on the middle section near the cavity given in Figure18, the injection process upstream of the cavity was further analyzed. The process can be decomposed into three typical physical processes: transverse jet flow, cavity flow, and space-restricted flow. The mixing of fuel and main flow in the transverse jet is mainly divided into three stages: The transverse jet enters the flow field vertically from the main-flow and deflects to the main-flow; the formation and development of large scale vortexes; and the breaking of large scale vortexes. The shear layer appears above the cavity, while the low-speed recirculation zone is generated in the cavity. In addition, for space-restricted flow, the shock system induced by the shock/boundary layer at the upper wall will have a certain impact on the main flow.Figure 18
Schematic diagram of the shock waves and the vortex structures in a cavity.Due to the large initial dynamic pressure of the jet flow and the small dynamic pressure of the near-wall gas under the viscous effect of the wall, the initial jet flow is basically perpendicular to the main-flow direction and expands further after entering the main-flow, forming a typical “bucket” shock and Mach disc, and tends to the downstream under the impact of the incoming flow. At the same time, a typical “λ-shock” and recirculation zone are formed upstream of the nozzle. At this stage, the jet is cylindrical and has poor mixing with the main-flow.As the jet moves towards the main-flow direction, under the shear effect of the jet and main-flow, the shear layer begins to stretch and deform. Meanwhile, a large vortex structure is formed, which increases the contact zone between the jet and main-flow and improves the mixing efficiency. Because the deceleration and pressurization effect of the cavity counteracts the expansion effect of the cavity, the shear layer is “lifted” to the main-flow by the cavity. Under the effect of the instability of the shear layer and the small vortex structure in the cavity, the large vortex structure in the shear layer is gradually broken into small vortexes. Although the shear layer is parallel to the main-flow direction in the cavity, the air can still be mixed well with the fuel in the cavity due to the transport of small vortexes.Since the height of the combustion chamber is limited, the “λ-shock” formed upstream of the nozzle can hit the upper wall of the combustion chamber, forming a typical shock/boundary layer interference. At the same time, the boundary layer separation is generated near the incident point of the shock wave, while the reflected shock wave hits the shear layer, increasing the instability of the shear layer. It promotes the breaking of large vortexes and improves the mixing efficiency of fuel and oxidant, which is beneficial to the ignition of the fuel and the maintenance of the flame.


### 4.3. Evolution of Fuel Distribution along the Flow Direction

Figure19 shows the three-dimensional distribution of the mean mixing fraction along the flow direction. In order to facilitate the analysis, eight sections (M1–M8) are given along the nozzle to cavity. It can be seen that when the fuel flows downstream, vortex pairs are formed in the Z-direction. In M1–M3, the fuel is weakly mixed with the incoming flow at the initial stage of injection. However, the distribution of mean mixing fraction in M4–M8 fluctuates greatly, indicating that the flow in this zone is violent and the mixing is complex. Figure 20 shows the distribution of streamline and vortex in the middle section. When the main-flow interacts with the jet, vortex structures of different scales will be formed. In the flow direction, the main flow is blocked by the jet and bypasses the jet from the side, and is “lifted” by the high-pressure zone behind the jet to form a large vortex structure. The large vortex structure of the shear layer is broken into a small vortex structure under the combined effect of the instability of the shear layer and the cavity, which promotes mixing. This is also the reason for the uneven distribution of mean mixing fraction in M4–M8 shown in Figure 19.Figure 19
Mean mixture fraction distribution along theX-direction.Figure 20
Streamline near the nozzle.


### 4.4. Combustion Evolution

Figure21 shows the transient temperature distribution in the middle section. As can be seen from Figure 21(a), combustion is first realized at the bottom of the cavity with good fuel mixing, low speed, and high temperature. The shock wave induced by the jet hits the upper wall surface and is reflected on the shear layer at the rear of the combustion chamber. At the same time, the Kelvin–Helmholtz instability caused by the difference of the velocity between two sides of the shear layer causes the fluctuation of the shear layer. The change caused by fluctuation induces a series of weak shock waves, which further aggravates the fluctuation of the shear layer. As shown in Figures 21(b)–21(f), this coupling effect of the shock waves and the shear layer rapidly destabilizes the shear layer and transforms it into a turbulent flow, which promotes the mixing of air and fuel.Figure 21
Transient temperature in the middle section. (a) 0.2 ms. (b) 0.3 ms. (c) 0.4 ms. (d) 0.5 ms. (e) 0.6 ms. (f) 0.7 ms.
(a)(b)(c)(d)(e)(f)Since the combustion zone and the near-wall zone are both subsonic, the flame spreads upstream along the combustion zone between the shear layer and the wall. As the combustion continues, the shear layer is “lifted” by the high pressure generated in the combustion zone, and the phenomenon of “flame lifting” appears. The downstream airflow channel becomes narrower, which forms stronger shock waves in the flow field. The strong shock waves increase the adverse pressure gradient of the upper wall, expand the separation zone, and induce stronger reflected shock waves. When the pressure rise and oblique shock waves generated by combustion are matched with the supersonic incoming flow, the flow field reaches dynamic equilibrium and appears as a small oscillation in the combustion.


### 4.5. Distribution of Pressure and Velocity

Figure22 shows the three-dimensional pressure distribution. It can be seen from Figure 22(a) that the pressure distribution is basically symmetrical in the Z-direction, but there is an obvious gradient in the Y-direction, especially in S3 and S4. Pressure in the upper part of S3 is higher, while the trends of S4 and S5 are opposite. In contrast, the changes of S6 and S7 are relatively flat. As shown in Figure 22(b), sudden changes in pressure exist on the upper wall. The pulsation with low frequency appears in the middle shear layer, and the pressure in the cavity shows a certain degree of unevenness. The analysis believes that the separation of the boundary layer on the upper wall is caused by the intersection of the arc shock wave and the reflected shock wave, resulting in the formation and development of turbulence near the shear layer and in the cavity. This process enhances the mixing of fuel and air, which is beneficial for combustion.Figure 22
Pressure distribution. (a) S1–S7. (b) Middle section.
(a)(b)Figure23 shows the three-dimensional X-velocity distribution. As shown in Figure 23(a), the inlet velocity reaches 1600 m/s under the boundary condition. The velocity decreases after the shock waves. From S3 to S7, it can be observed that the “thermal throat” formed in the middle part and lower part gradually rises upward. The velocity in the thermal throat is about 800–1200 m/s, but it is significantly higher than the surrounding. This phenomenon is consistent with the actual description of the thermal throat in a cavity-based combustor [61]. Combined with Figure 23(b), it can be seen that the thermal throat is not stable, and the change is mainly reflected in the Y-direction. The air not participating in the reaction mainly flows backward through the thermal throat. The instability of the thermal throat may be caused by the instability of the boundary layer and the shear layer. If the disturbance is enhanced, the thermal throat may further narrow, or even reduce the flow to subsonic flow, resulting in blockage. Therefore, attention should be paid to the control of boundary layer separation in engine design. Proper boundary layer separation will facilitate the mixing of air and fuel, but serious boundary layer separation will produce strong shock waves, which may lead to the blockage of the channel.Figure 23
X-velocity distribution. (a) S1–S7. (b) Middle section.
(a)(b)


## 4.1. Transient Flame Structure Based on the OH-PLIF Image

In order to study the transient flame structure in the flow field of combustion, the simulated concentration distribution diagram of the OH group was compared with the OH-PLIF image of the experiment in reference [21]. Figures 16 and 17 show the distribution of the OH group in the middle section obtained by experiment and simulation, respectively. The analysis shows that the flame distribution of the lateral injection on the upstream wall of the cavity has the following characteristics: The flame is mainly located in the upper part of the cavity and the mainflow near the rear edge of the cavity. The simulation also captures this phenomenon and finds the distribution of the OH group near the nozzle. It means that part of the fuel is entrained into the reflux zone upstream of the nozzle and combusted, which is different from the experiment. At the same time, the simulation observes that the combustion inside the cavity is relatively intense, while the experimental observation shows that the combustion on the lower wall of the cavity is weak, for two reasons: Due to the weak spontaneous emission signal of the OH group in the flame, a long exposure time (50 μs in the experiment) is needed to obtain the spontaneous emission image with high signal-to-noise ratio. To sum up, the time-averaged flame image is obtained by spontaneous emission, rather than the transient flame image. In addition, the OH-PLIF image is the projection of a three-dimensional entity on a two-dimensional plane, so it cannot detect the information on a single section. However, the simulation results in Figure 17 are obtained from a two-dimensional slice of the middle section, which is a transient flame that cannot be observed through experiments. To sum up, there are certain differences between the experimental and simulation results. In addition, Fan et al. also emphasized in their work that the combustion has an unsteady effect, which will make it difficult to control the measurement precision of the OH-PLIF technology, and the measurement results may have certain errors.Figure 16
OH-PLIF image of the experiment [21].Figure 17
OH mass concentration of the simulation.


## 4.2. Analysis of the Injection Process Upstream of the Cavity

Combined with the schematic diagram of the shock waves and the vortex structure on the middle section near the cavity given in Figure18, the injection process upstream of the cavity was further analyzed. The process can be decomposed into three typical physical processes: transverse jet flow, cavity flow, and space-restricted flow. The mixing of fuel and main flow in the transverse jet is mainly divided into three stages: The transverse jet enters the flow field vertically from the main-flow and deflects to the main-flow; the formation and development of large scale vortexes; and the breaking of large scale vortexes. The shear layer appears above the cavity, while the low-speed recirculation zone is generated in the cavity. In addition, for space-restricted flow, the shock system induced by the shock/boundary layer at the upper wall will have a certain impact on the main flow.Figure 18
Schematic diagram of the shock waves and the vortex structures in a cavity.Due to the large initial dynamic pressure of the jet flow and the small dynamic pressure of the near-wall gas under the viscous effect of the wall, the initial jet flow is basically perpendicular to the main-flow direction and expands further after entering the main-flow, forming a typical “bucket” shock and Mach disc, and tends to the downstream under the impact of the incoming flow. At the same time, a typical “λ-shock” and recirculation zone are formed upstream of the nozzle. At this stage, the jet is cylindrical and has poor mixing with the main-flow.As the jet moves towards the main-flow direction, under the shear effect of the jet and main-flow, the shear layer begins to stretch and deform. Meanwhile, a large vortex structure is formed, which increases the contact zone between the jet and main-flow and improves the mixing efficiency. Because the deceleration and pressurization effect of the cavity counteracts the expansion effect of the cavity, the shear layer is “lifted” to the main-flow by the cavity. Under the effect of the instability of the shear layer and the small vortex structure in the cavity, the large vortex structure in the shear layer is gradually broken into small vortexes. Although the shear layer is parallel to the main-flow direction in the cavity, the air can still be mixed well with the fuel in the cavity due to the transport of small vortexes.Since the height of the combustion chamber is limited, the “λ-shock” formed upstream of the nozzle can hit the upper wall of the combustion chamber, forming a typical shock/boundary layer interference. At the same time, the boundary layer separation is generated near the incident point of the shock wave, while the reflected shock wave hits the shear layer, increasing the instability of the shear layer. It promotes the breaking of large vortexes and improves the mixing efficiency of fuel and oxidant, which is beneficial to the ignition of the fuel and the maintenance of the flame.


## 4.3. Evolution of Fuel Distribution along the Flow Direction

Figure19 shows the three-dimensional distribution of the mean mixing fraction along the flow direction. In order to facilitate the analysis, eight sections (M1–M8) are given along the nozzle to cavity. It can be seen that when the fuel flows downstream, vortex pairs are formed in the Z-direction. In M1–M3, the fuel is weakly mixed with the incoming flow at the initial stage of injection. However, the distribution of mean mixing fraction in M4–M8 fluctuates greatly, indicating that the flow in this zone is violent and the mixing is complex. Figure 20 shows the distribution of streamline and vortex in the middle section. When the main-flow interacts with the jet, vortex structures of different scales will be formed. In the flow direction, the main flow is blocked by the jet and bypasses the jet from the side, and is “lifted” by the high-pressure zone behind the jet to form a large vortex structure. The large vortex structure of the shear layer is broken into a small vortex structure under the combined effect of the instability of the shear layer and the cavity, which promotes mixing. This is also the reason for the uneven distribution of mean mixing fraction in M4–M8 shown in Figure 19.Figure 19
Mean mixture fraction distribution along theX-direction.Figure 20
Streamline near the nozzle.


## 4.4. Combustion Evolution

Figure21 shows the transient temperature distribution in the middle section. As can be seen from Figure 21(a), combustion is first realized at the bottom of the cavity with good fuel mixing, low speed, and high temperature. The shock wave induced by the jet hits the upper wall surface and is reflected on the shear layer at the rear of the combustion chamber. At the same time, the Kelvin–Helmholtz instability caused by the difference of the velocity between two sides of the shear layer causes the fluctuation of the shear layer. The change caused by fluctuation induces a series of weak shock waves, which further aggravates the fluctuation of the shear layer. As shown in Figures 21(b)–21(f), this coupling effect of the shock waves and the shear layer rapidly destabilizes the shear layer and transforms it into a turbulent flow, which promotes the mixing of air and fuel.Figure 21
Transient temperature in the middle section. (a) 0.2 ms. (b) 0.3 ms. (c) 0.4 ms. (d) 0.5 ms. (e) 0.6 ms. (f) 0.7 ms.
(a)(b)(c)(d)(e)(f)Since the combustion zone and the near-wall zone are both subsonic, the flame spreads upstream along the combustion zone between the shear layer and the wall. As the combustion continues, the shear layer is “lifted” by the high pressure generated in the combustion zone, and the phenomenon of “flame lifting” appears. The downstream airflow channel becomes narrower, which forms stronger shock waves in the flow field. The strong shock waves increase the adverse pressure gradient of the upper wall, expand the separation zone, and induce stronger reflected shock waves. When the pressure rise and oblique shock waves generated by combustion are matched with the supersonic incoming flow, the flow field reaches dynamic equilibrium and appears as a small oscillation in the combustion.


## 4.5. Distribution of Pressure and Velocity

Figure22 shows the three-dimensional pressure distribution. It can be seen from Figure 22(a) that the pressure distribution is basically symmetrical in the Z-direction, but there is an obvious gradient in the Y-direction, especially in S3 and S4. Pressure in the upper part of S3 is higher, while the trends of S4 and S5 are opposite. In contrast, the changes of S6 and S7 are relatively flat. As shown in Figure 22(b), sudden changes in pressure exist on the upper wall. The pulsation with low frequency appears in the middle shear layer, and the pressure in the cavity shows a certain degree of unevenness. The analysis believes that the separation of the boundary layer on the upper wall is caused by the intersection of the arc shock wave and the reflected shock wave, resulting in the formation and development of turbulence near the shear layer and in the cavity. This process enhances the mixing of fuel and air, which is beneficial for combustion.Figure 22
Pressure distribution. (a) S1–S7. (b) Middle section.
(a)(b)Figure23 shows the three-dimensional X-velocity distribution. As shown in Figure 23(a), the inlet velocity reaches 1600 m/s under the boundary condition. The velocity decreases after the shock waves. From S3 to S7, it can be observed that the “thermal throat” formed in the middle part and lower part gradually rises upward. The velocity in the thermal throat is about 800–1200 m/s, but it is significantly higher than the surrounding. This phenomenon is consistent with the actual description of the thermal throat in a cavity-based combustor [61]. Combined with Figure 23(b), it can be seen that the thermal throat is not stable, and the change is mainly reflected in the Y-direction. The air not participating in the reaction mainly flows backward through the thermal throat. The instability of the thermal throat may be caused by the instability of the boundary layer and the shear layer. If the disturbance is enhanced, the thermal throat may further narrow, or even reduce the flow to subsonic flow, resulting in blockage. Therefore, attention should be paid to the control of boundary layer separation in engine design. Proper boundary layer separation will facilitate the mixing of air and fuel, but serious boundary layer separation will produce strong shock waves, which may lead to the blockage of the channel.Figure 23
X-velocity distribution. (a) S1–S7. (b) Middle section.
(a)(b)


## 5. Conclusions

In this paper, a compressible modified FGM model suitable for supersonic combustion is established by improving the method of solving compressible flow field parameters. The integrity and stability of the model are verified by simulating the combustion flow field of a typical DLR scramjet combustion chamber. Combining with hybrid RANS/LES turbulence model, the unsteady combustion process of hydrogen transverse jet in a cavity-based scramjet was numerically simulated. The conclusions of this study can be summarized as follows:(1)
The compressible modified FGM model can capture the fine shock waves and process of the flame propagation in the flow field, which is consistent with the experimental results. In addition, the phenomenon of “flame lifting” can be observed in the results of the simulation.(2)
There is a strong interaction between the shock wave and the shear layer. The disturbance of the shock wave destabilizes the shear layer and transforms into a turbulent flow, which is conducive for the mixing of air and fuel;(3)
The separation of the boundary layer on the upper wall of the combustion chamber will reduce the stability of the shear layer and promote the development of turbulence;(4)
When a scramjet engine is designed, the control of boundary layer separation should be considered. Proper boundary layer separation will facilitate the mixing of air and fuel, but serious boundary layer separation will produce strong shock waves, which may lead to the blockage of the channel.



---

*Source: 8388827-2022-02-08.xml*