Large Eddy Simulation and Dynamic Mode Decomposition of Supersonic Combustion Instability in a Strut-Based Scramjet Combustor

Abstract

Supersonic combustion instability studies are crucial for the future maturation of scramjet engines. In the present paper, the supersonic combustion instability in a strut-based scramjet combustor is investigated through large eddy simulation and dynamic mode decomposition. The results show significant pressure oscillation in the strut-based scramjet combustor when the air parameters at the combustor inlet and the fuel parameters at the injector outlet are under certain conditions, and these pressure oscillation situations correspond to supersonic combustion instability. The oscillations have multiple dominant frequencies, including relatively low frequency of 2984 Hz, high frequency of 62,180 Hz, and very high frequency of 110,562 Hz. Large pressure oscillations in the strut-based scramjet combustor are closely related to wake instability, shear layer instability, shear layer and wave interactions, and combustion. Reducing the air total temperature at the combustor inlet can attenuate the pressure oscillations, and reducing the fuel flow rate at the injector outlet can also attenuate the pressure oscillations.

1. Introduction

Supersonic combustion is an important phenomenon in hypersonic air-breathing propulsion systems, and corresponding research has made great progress in the past decades [1,2]. Supersonic combustion in cavity-based combustors [3,4,5,6,7,8,9,10], strut-based combustors [11,12,13,14,15,16,17,18], and mixing layers [19,20] have all been investigated in detail. At present, unsteady supersonic combustion has become an issue of great concern in the field of supersonic combustion. This is especially true for supersonic combustion instability, which has become a topical subject of research [21,22]. The causes of supersonic combustion instability are complex, and the study of supersonic combustion instability is important for all types of supersonic combustors with different configurations [22].

The cavity configuration has the advantages of simple structure, high flame stability, low total pressure loss and drag [23,24], and is one of the most widely used supersonic combustor configurations. Initially, it was believed that acoustic waves would not propagate upstream in a supersonic combustor, and thus, a closed-loop feedback loop between acoustic waves and combustion heat release could not be formed. Consequently, the problem of supersonic combustion instability has been ignored. However, studies on cavity configuration supersonic combustors have shown that there are multiple subsonic flow regions inside the combustor, which allow the disturbance-induced acoustic waves to propagate upstream. Therefore, combustion stability can no longer be ignored [25]. Choi et al. [26,27,28] performed numerical simulations of flow oscillations in a scramjet with and without a cavity and investigated the effect of the interaction between the fuel jet, the shock wave structure, and the boundary layer on the unsteady flow characteristics. They found that the coupling between flow oscillations and unsteady heat release in the reaction zone leads to combustion instabilities. Using ethylene as fuel, Tian et al. [29] investigated the effect of inlet Mach number and transverse injection on combustion instability and found the presence of a high frequency oscillation peak frequency of about 4 kHz at the trailing edge of the cavity. Similarly, Gao et al. [30] investigated the effect on combustion instability when the equivalent ratio was in the range of 0.2 to 0.63. Wang et al. [31] reported the effect of the equivalence ratio on flame oscillation in a hydrogen-fueled cavity engine. The authors found that an increase in equivalence ratio leads to a shift in the combustion stabilization mode from cavity to jet-trail mode, with less flame front fluctuation in the mode and enhanced flame front oscillation in the 1000–1500 Hz range from the cavity to the jet-trail stabilization mode. Peng et al. [32] further examined the effect of different equivalence ratios on combustion instability in the spanwise direction. They found that the peak oscillation frequency in the super-combustion mode lies between 150–200 Hz and the low-frequency oscillation becomes more pronounced with increasing equivalence ratio. Ouyang et al. [33] studied the effects of the number and position of cavities, length-to-depth ratio, and rear inclination angle on combustion instability. It was found that, compared with other parameters, the cavity position has a greater effect on the flame, and the premixing distance of ethylene and oxygen varies with different cavity positions. The flame will show two forms of small high-frequency oscillations and large low-frequency oscillations with the change of premixing distance. It is difficult to change the main frequency of combustion oscillations by increasing the number of cavities provided that the cavity position is sufficient to make the fuel and oxidizer premix. Zhao et al. [34] evaluated the effects of length-to-depth ratio, rear inclination angle, and air throttling downstream of the cavity on combustion oscillations. The authors found that quasi-periodic combustion oscillations occurred when the length-to-depth ratio and rear inclination angle were large and the air throttling downstream of the cavity was close. Ma et al. [35] studied the thermoacoustic instability of the cavity configuration combustor using a quasi one-dimensional model. They concluded that the instability is related to the acoustic–convective interaction between the injector and the flame zone. Lin et al. [36] experimentally investigated the mechanism of thermoacoustic instability in a cavity configuration supersonic combustor with different injection schemes. It was shown that the main frequency increases with increasing the equivalent ratio in the range of 100–400 Hz. This phenomenon may be due to the acoustic feedback loop between the excitation and flame regions and the acoustic–convective feedback loop between the fuel injection and the flame region. A similar feedback loop was found by Allison et al. [37]. They performed a fast Fourier transform (FFT) analysis of CH chemiluminescence images and found a characteristic frequency of 340 Hz. They inferred that the frequency is related to the instability of the reflection and convection of acoustic waves between the shock wave trains and the flame front. Nakaya et al. [38] used fast Fourier and dynamic mode decomposition (DMD) techniques to find two different characteristic frequencies of 128 Hz and 1600 Hz. The occurrence of the 128 Hz peak frequency was associated with the whole flame oscillation in the steady combustion of the cavity shear layer, while the other oscillation frequency was observed in the oscillation between the steady jet wake and ramming combustion. However, no characteristic frequency associated with acoustic feedback and convective acoustic feedback was found.

Combustion instability is also critical for supersonic combustors with strut configurations. Huang et al. [15] studied the pressure oscillations and flame front movement in the combustor of a DLR scramjet using a three-dimensional large eddy simulation. They suggested that the mechanisms causing pressure oscillations are different in different locations within the combustor. The flow instability and combustion instability dominated in the wake region near the strut and in the turbulent combustion region, respectively, while in the transition region, both were of equal status and jointly affected the pressure oscillations. In addition, the flame front moved backward and forward with time, indicating that the flame was inherently unstable. In the same group, Qin et al. [39] studied the combustion oscillations in the same combustor under different air inflow conditions, and found that the oscillation characteristics were quite different under these conditions. Li and Wang [40] further investigated the reaction flow dynamics within the DLR scramjet in detail using the DMD method. They provided the dominant modes of OH, HO₂, pressure, and velocity fields, and based on the results of the DMD analysis, they concluded that the oscillation of the pressure was related to the vortex motions in the near wake region of the strut and the interactions between the shock waves and the mixing layers. Yuan et al. [41] studied the combustor flame and thermoacoustic coupling characteristics using a thick flame model and the proper orthogonal decomposition (POD) technique. They found the presence of thermoacoustic coupling oscillations at a frequency of 4997 Hz, with the source of the thermoacoustic coupling disturbance coming from vortex shedding and the resulting downstream shock wave oscillations. Zhang et al. [42] investigated the flame oscillation characteristics in a dual-mode combustor with a strut plate configuration fueled by liquid kerosene. The results revealed that the flame oscillation phenomenon was caused by the variation of the pilot flame intensity and the inverse pressure gradient near the combustion zone.

In summary, the supersonic combustion instability of the cavity configuration scramjet combustor has gained much research attention as mentioned above. However, there are only limited studies on the supersonic combustion instability of the strut configuration scramjet combustor. This paper contributes to research on the supersonic combustion instability in a strut configuration scramjet combustor. The remainder of this paper is organized as follows. Section 2 outlines the numerical simulation method, Section 3 discusses the pressure oscillations and flow characteristics in supersonic combustion instability and the effect of air and fuel conditions on supersonic combustion instability, and Section 4 provides the conclusions.

2. Numerical Method

2.1. LES Equations and Chemical Kinetics

In the current study, the Large Eddy Simulation (LES) method is used to simulate the complex flow and combustion in the combustor of the DLR scramjet. This is a numerical simulation method of turbulence between Direct Numerical Simulation (DNS) and Reynolds Averaging (RANS). With the rapid improvement in computing hardware conditions, LES has become the main method for studying supersonic combustion in recent years [43,44,45,46]. The LES method performs direct numerical simulations for large scale motions that are strongly influenced by boundary conditions, while small scale vortices with more commonalities are simulated by constructing subgrid models. The transient Navier–Stokes equations are based on the mesh size and Favre averaging to obtain the following LES equations by applying a spatial filter:

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

(1)

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

(2)

$$\begin{array}{l}\frac{\partial (\overline{\rho }\tilde{e}+\overline{\rho }{\tilde{u}}_i{\tilde{u}}_i/2+k^{sgs})}{\partial t}+\frac{\partial [(\overline{\rho }\tilde{e}+\overline{\rho }{\tilde{u}}_i{\tilde{u}}_i/2+\overline{p}+k^{sgs}){\tilde{u}}_j]}{\partial x_j}=+\frac{\partial }{\partial x_j}\times \\ [{\tilde{u}}_j{\overline{\tau }}_{ij}+\overline{\rho }\nu \frac{\partial k^{\mathit{sgs}}}{\partial x_j}+{\overline{q}}_j-H_j^{sgs}+{\sigma }_{ij}^{sgs}]\end{array}$$

(3)

$$\frac{\partial (\overline{\rho }{\tilde{Y}}_m)}{\partial t}+\frac{\partial (\overline{\rho }{\tilde{Y}}_m{u}_j)}{\partial x_j}=\frac{\partial }{\partial x_j}(\overline{\rho }D\frac{\partial {\tilde{Y}}_m}{\partial x_j}-\overline{\rho }(\widetilde{u_j{Y}_m}-{\tilde{u}}_j{\tilde{Y}}_m))+{\overline{\dot{\omega }}}_m (m=1,\dots ,N)$$

(4)

where variables with overbar and tilde denote spatially filtered and Favre-averaged quantities, respectively. ρ is the density, (u_i)_i = 1,2,3 are the velocity components, p is the pressure, Y_m is the mass fraction of the mth species, N is the total number of species, $\lambda$ is the heat transfer coefficient, the molecular diffusion coefficient is evaluated from $D=\mu /\rho S_c$, S_c is the constant number, μ is the molecular dynamic viscosity, $H_j^{sgs}$ is the subgrid scale stress flux energy, ${\sigma }_{ij}^{sgs}$ is the subgrid scale stress viscous dissipation, k_sgs is the subgrid scale turbulent kinetic energy, and the Lewis number is assumed to be 1. The filtered heat flux ${\overline{q}}_i$ is given as follows:

$${\overline{q}}_j=\overline{\rho }\frac{\mu }{\mathrm{Pr}}\frac{\partial \tilde{h}}{\partial x_j}$$

(5)

The filtered subgrid scale stress (SGS) ${\overline{\tau }}_{ij}$ is given as:

$${\overline{\tau }}_{ij}=2{\mu }_t{\overline{S}}_{ij}-\frac{2}{3}{\mu }_t{\overline{S}}_{kk}{\delta }_{ij}$$

(6)

where $S_{ij}=\frac{1}{2}(\frac{\partial ui}{\partial xj}+\frac{\partial ui}{\partial xi})$ is the strain rate tensor, ${\delta }_{ij}$ is the Kronecker function, and ${\mu }_t$ is given as:

$${\mu }_t=C_k\rho k_{sgs}^{1/2}{\Delta }_f$$

(7)

In this paper, a one-equation eddy dynamic model is used for the closure of subgrid scale stresses [47,48,49]:

$$\rho \frac{\partial {\overline{k}}_{sgs}}{\partial t}+\rho \frac{\partial ({\overline{u}}_j{\overline{k}}_{sgs})}{\partial x_j}=-{\tau }_{ij}\frac{\partial {\overline{u}}_i}{\partial x_j}-C_{\epsilon }\rho \frac{k_{sgs}^{3/2}}{{\Delta }_f}+\frac{\partial }{\partial x_j}\left(\frac{{\mu }_t}{{\sigma }_k}\frac{\partial k_{sgs}}{\partial x_j}\right)$$

(8)

where the constants $C_k$ and $C_{\epsilon }$ are dynamically determined and ${\sigma }_k$ are fixed at 1.0.

Due to the time and money costs involved in calculating the supersonic combustion of hydrogen-air using a multi-step reaction chemistry, this paper uses a one-step total package reaction mechanism. The combustion model used is the Finite-Rate/Eddy-Dissipation model, which is a combination of the Arrhenius equation and vortex dissipation, taking into account the effects of both chemical reaction kinetics and turbulent fluctuations. This model provides more accurate calculations for hypersonic combustion processes.

2.2. Dynamic Mode Decomposition

The supersonic combustion processes within the DLR scramjet are complex, exhibiting complex dynamics over a wide range of spatial and temporal scales. To further improve understanding of the underlying mechanisms of supersonic combustion instabilities, the dynamic mode decomposition method (DMD) is used to quantify the flow field. DMD is a purely data-driven algorithm that extracts dynamically coherent features as spatial modes, and each mode has a fixed oscillation frequency associated with it. A detailed report on the DMD algorithm can be found in the literature [50,51].

An nth-order matrix snapshot X can be obtained by arranging the n moments of flow field data obtained from the numerical simulation into a column in time order, and the time interval between two adjacent snapshots is ∆t. The nth-order matrix snapshots are denoted as matrices X₁ and X₂ respectively:

$$X_1=\left[x_1,x_2,\cdots ,x_{n-1}\right],X_2=\left[x_2,x_3,\cdots ,x_n\right]$$

(9)

It is assumed that the flow field x_i+1 and x_i can be mapped by the following linear relationship:

$$x_{i+1}=A{x}_i$$

(10)

where A is a linear operator whose eigenvalues and eigenvectors describe the dynamic properties of the flow field. For non-linear dynamic systems, Equation (16) is a linear estimation mapping relation. Thus, the relationship between the matrices X₁ and X₂ can be written as:

$$X_2=\left[A{x}_1,A{x}_2,\cdots ,A{x}_{n-1}\right]=A{X}_1$$

(11)

Due to the large amount of data obtained from the numerical simulation, it is necessary to find a low-dimensional matrix to replace matrix A. The singular value decomposition method is used to reconstruct X₁ as follows:

$$X_1=UE{V}^T$$

(12)

$$A=U\tilde{A}{U}^T$$

(13)

where the matrix U is a left singular vector, V is a right singular vector, E is a diagonal matrix containing r singular values, and the accompanying matrix $\tilde{A}$ is given as:

$$\tilde{A}=U^H{X}_{2}V{E}^{-1}\approx A$$

(14)

According to $\tilde{A}W=W\Lambda$, the eigenvalues λ_j and eigenvectors w_j of $\tilde{A}$ are obtained, so that the jth DMD mode is given as:

$${\Phi }_j=U{w}_j$$

(15)

The real part of the logarithmic form of the eigenvalue λ_j represents the growth/decay rate g_j of the DMD mode, and the imaginary part represents the frequency ω_j of the DMD mode:

$$g_j=\mathrm{Re}(lg({\lambda }_j)/\Delta t)$$

(16)

$${\omega }_j=Im(lg({\lambda }_j)/\Delta t)$$

(17)

Mode stability can be judged by the distribution of eigenvalues and the value of the modal growth rate. If the growth rate is positive, this means that the corresponding mode is unstable and the eigenvalues corresponding to that mode fall outside the unit circle; if the growth rate is negative, the corresponding mode is stable and the eigenvalues corresponding to that mode lie inside the unit circle; if the growth rate is zero, the corresponding mode is periodic and the eigenvalues corresponding to that mode lie on the unit circle.

W is a matrix whose column vectors are eigenvectors w_j. The eigen-decomposition of A can be expressed as:

$$\tilde{A}=WN{W}^{-1} N=diag({\lambda }_1,\cdots ,{\lambda }_r)$$

(18)

Thus, the snapshot at any moment X is estimated by giving:

$$x_i=A{x}_{i-1}=U\tilde{A}{U}^T{x}_{i-1}=UWN{W}^{-1}{U}^T{x}_{i-1}=UW{N}^{i-1}{W}^{-1}{U}^T{x}_1$$

(19)

According to Equation (15), each column of Φ is defined as a DMD mode:

$$\Phi =UW$$

(20)

The mode amplitude α is given as a criterion for judging the level of mode contribution to the flow field:

$$\alpha =W^{-1}{U}^T{x}_1,\alpha =\left[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_{\mathrm{r}}\right]$$

(21)

The real flow field can be reconstructed from DMD modes, mode amplitudes, and eigenvalues, and even predicted for the development of the flow field. Simply taking Equations (20) and (21) into Equation (19) gives:

$$x_i=\Phi {\Lambda }^{i-1}\alpha ={\displaystyle \sum _{j=1}^r{\Phi }_j{({\lambda }_j)}^{i-1}{\alpha }_j}$$

(22)

The parametric number of each DMD mode indicates the magnitude of the energy of the mode. The magnitude of the modal energy can also be judged by the magnitude of the contribution of the mode to the flow field, using the following equation for mode energy:

$$\Vert {\Phi }_j\Vert =\sqrt{{\displaystyle \sum _{i=1}^n{{\displaystyle \left|{\Phi }_j\right|}}^2}}$$

(23)

2.3. Computational Configuration and Flow Conditions

2.3.1. Geometry Details and Computational Domain

The geometric model of the DLR scramjet combustor in this paper is in agreement with Waidmann et al. [52,53,54]. Figure 1 illustrates a two-dimensional schematic of the combustor with a total length of 340 mm, an air inlet height of 50 mm, and a spanwise length of 2.4 mm. The upper wall surface is deflected outward at a 3° divergence angle from a position 100 mm from the air inlet. The length of the wedge-shaped strut on the centerline of the combustor is 32 mm, and its apex angle is 12°. The coordinates’ origin O is located at the bottom center of the strut at a distance of 109 mm from the air inlet. In the experiments of Waidmann et al. [52,53,54], the bottom of the strut has a row of 15 jet holes with a 1 mm diameter and a 1.4 mm pitch. In this study, only one hydrogen jet port is selected due to the computational cost.

2.3.2. Computational Mesh

Figure 2 shows a multi-block structured mesh of the entire 3D combustor, which is referenced from the setup of Génin and Menon [55]. The combustor is decomposed into 318 and 32 nodes in the x-streamwise direction and the z-spanwise direction, respectively, 56 nodes between the upper wall surface of the strut and the upper wall surface of the combustor, 58 nodes between the lower wall surface of the strut and the lower wall surface of the combustor, and 68 nodes in the y-transverse direction at the bottom of the strut. In order to capture the complex flow and combustion phenomena near the strut, the mesh around the strut, hydrogen injection, and the centerline of the combustor is locally refined, as shown in Figure 2b. The minimum distance in the direction perpendicular to the wall of the strut is set to 0.008 mm, so that the y⁺ of the first layer grid is less than 1. The number of computational meshes used for solving the supersonic combustion reaction process is approximately 1.5 million. This grid is used for the computations of the seven cases listed in

Table 1. Air and fuel parameters of different cases.

	Inlet	Ma	Total Temperature K	Velocity m/s	Static Pressure Pa	Mass Flux kg/(m2/s)
1 (baseline case)	air	2.0	568	730	100,000	-
	H2	1.0	300	1200	-	116
2	air	2.0	460	657	90,000	-
	H2	1.0	300	1200	-	116
3	air	2.0	960	949	130,000	-
	H2	1.0	300	1200	-	116
4	air	2.0	568	730	100,000	-
	H2	1.0	300	1200	-	232
5	air	2.0	568	730	100,000	-
	H2	1.0	300	1200	-	58
6	air	2.0	960	949	130,000	-
	H2	1.0	300	1200	-	232
7	air	2.0	960	949	130,000	-
	H2	1.0	300	1200	-	58

to provide preliminary results and comparisons between cases for the supersonic combustion instability community, because the amount of grid points is not very large by current standards.

2.3.3. Boundary Conditions

The pressure inlet and mass flow boundary conditions are used for the combustor air inlet and the strut hydrogen injection, respectively, and the pressure outlet boundary condition is used for the combustor outlet. Seven flow conditions are simulated in this paper, and the air and fuel parameters for the baseline case are kept consistent with the literature [56], with Ma = 2 for the air inlet and Ma = 1 for the hydrogen injection. Detailed inlet and outlet condition parameters are provided in Table 1. Without considering the effect of the upper and lower wall boundary layers, wall adiabatic slip wall boundary conditions are used on the upper and lower walls. Periodic boundary conditions are used on both side walls. No-slip adiabatic wall boundary conditions are used for the strut walls.

2.3.4. Numerical Methods

Numerical computations are performed using the density-based and double precision solver. The large eddy simulation equations are solved using the finite volume method and the second-order implicit algorithm. The convective terms are discretized using the AUSM scheme, which has a strong ability to capture shock waves and discontinuities. The spatial gradient is solved using the Least Squares Cell-Based method. In order to accelerate the process of getting usable LES unsteady data, the LES unsteady simulations are performed using a convergent flow field (obtained from the RANS steady method) as an initial field. In total, 2.4 ms (approximately 8 flow-through times) of the flow are computed. The maximum CFL (Courant–Friedrichs–Lewy) number of the reactive flow is set to 0.1 to ensure the stability of the simulation, corresponding to a calculation time step of 3 × 10⁻⁸ s. Pressure data at sampling points are collected at every time step (i.e., with a time interval of 3 × 10⁻⁸ s). A total of 200 snapshots in the seventh and the eighth flow-through times with a time interval of 3 × 10⁻⁶ s are collected from the pressure field, velocity field, and temperature field for conducting DMD.

2.4. Numerical Validation

Before studying the supersonic combustion instability characteristics of the DLR scramjet using large eddy simulation, we first verify the numerical reactive flow results of the baseline case compared with the available experimental [53] and numerical [55] results. The mean velocity and temperature at six flow positions are examined.

Figure 3 compares the x-streamwise direction mean velocity U and temperature T at six different flow locations: x = 11 mm, 58 mm, 90 mm, 115 mm, 140 mm, and 166 mm. As shown in Figure 3a, the mean velocity peaks at the flow positions x = 90 mm and x = 115 mm deviate considerably from those of Waidmann et al. [53] and Génin et al. [55]. The most likely reason for this result may come from the fact that a one-step reaction mechanism is used in this paper. The mean velocities at the remaining flow positions in general agree with the calculations of Génin and Menon [55], but there are small discrepancies with the experiments. For the baseline case, it can be seen from Figure 3b that the mean temperature of the six flow positions is in good agreement with the results of Génin and Menon [55] and the peaks are close, which verifies the accuracy of the numerical model. The comparison of the above data verifies the correctness and accuracy of the calculation method in this study.

3. Results and Discussion

The combustion and flow characteristics are obtained by conducting large eddy simulations and dynamic mode decomposition for the seven cases in Table 1. In the following section, the typical combustion instability case, case 3, and the case without combustion instability, case 5, are discussed in Section 3.1. Then, Section 3.2 discusses the effect of air conditions on combustion instability by comparing cases 3 and 1, cases 6 and 4, and cases 7 and 5, while the effect of fuel conditions on combustion instability is discussed by comparing cases 4, 1, and 5, and cases 6, 3, and 7.

3.1. Pressure Oscillation and Flow Characteristics in Supersonic Combustion Instability

Pressure oscillations are an important feature of supersonic combustion oscillations. In this paper, the parameters including pressure are monitored at 18 points in the flow field of the supersonic combustor for each case in Table 1. By comparing the pressure evolution of the monitoring points for each case, it is found that there are significant differences in the pressure evolution of these cases. The pressure evolution profiles of the monitoring points for case 3, in which typical pressure oscillations exist, and the pressure evolution profiles of the monitoring points for case 5, in which no pressure oscillations exist, are shown in Figure 4 and Figure 5, respectively. The 1, 2…18 indicate the monitoring points, and the locations of these monitoring points are shown in Figure 6 and Figure 7 (the background of the figure is a transient pressure contour).

As illustrated in Figure 4 and Figure 6, for case 3, the pressure at the monitoring points p1, p3, p7, and p8 remain essentially unchanged. The pressure at the monitoring points p4, p5, and p6 experience small fluctuations. The pressure at the monitoring points p2, p9, p10, p11, p12, p13, p14, p15, p16, p17, and p18 have dramatic fluctuations. The physics behind these phenomena are explained as follows: the p1 is located in the region undisturbed by the strut, p3 is located in the region after the shock wave generated by the head of the strut, and p7 and p8 are located in the region behind the reflected shock wave generated by the head of the strut and the reflection from the upper and lower walls, and the flow in these regions is stable. The p4 is located in the region of the expansion wave fan generated at the bottom of the strut, p5 and p6 are located in the near wake at the bottom of the strut, and the flow in these regions is mainly influenced by wake instability. The p2, p9, p10, p12, and p13 are located in the region of the combustion shear layer in the strut wake, p14 and p17 are located in the region of the unburned shear layer downstream of the combustion shear layer in the strut wake, and p11, p15, p16, and p18 are located in the regions after the “shock system originating from the combustion shear layer in the strut wake”, where the flow is affected by wake instability, shear layer instability, and combustion.

For case 5, as shown in Figure 5 and Figure 7, the pressure at monitoring points p1, p3, p7, and p8 is maintained unchanged, while the pressure fluctuates slightly at monitoring points p2, p4, p5, p6, p9, p10, p11, p12, p13, p14, p15, p16, p17, and p18, among which the pressure characteristics at monitoring points p2, p9, p10, p11, p12, p13, p14, p15, p16, p17, and p18 are significantly different from those of case 3. This indicates that different combustion characteristics match different pressure characteristics.

For case 3, points p5 and p6 with small pressure fluctuations and points p9 and p13 with strong pressure fluctuations are selected for further analysis using the FFT (Fast Fourier Transform) method. As can be observed from the frequency-amplitude profiles in Figure 8, the peak frequencies correspond to relatively large amplitude at points p5 and p6 are mainly below 20,000 Hz, while the peak frequencies correspond to relatively large amplitude at points p9 and p13 include frequencies below 20,000 Hz and frequencies around 50,000 Hz. The frequency of the maximum peak at point p9 is 48,333 Hz, which is exactly close to the frequency of the quasi-periodic occurrence of peaks in the pressure evolution profiles at point p9 (see Figure 4a). Therefore, it is necessary to discuss in detail the evolution of the pressure at point p9 between the two adjacent peaks.

The period near the moment 1.816 ms where the highest pressure peak is located in Figure 4a and the periods near the moments 2.166 ms and 2.187 ms where the two secondary pressure peaks are located are selected to observe the pressure evolution near the p9 point. Figure 9 shows the pressure evolution near the p9 point from 1.807 to 1.828 ms and from 2.166 to 2.187 ms. As seen in Figure 9, the lower pressure region upstream of the p9 point periodically forms a local high-pressure region behind the shock wave. This local high-pressure region gradually convects downstream with time (as shown by the black circles at each moment in Figure 9) while the pressure gradually increases. The pressure in this high-pressure region already has a high value when convection reaches the p9 point position, forming the pressure peak shown in Figure 4.

Further, Figure 10 shows the evolution of H₂ mass fraction near the p9 point from 1.807 to 1.828 ms and from 2.166 to 2.187 ms (the position of the black circle in Figure 10 is the same as the position of the black circle in the pressure contour at the same moment in Figure 9). It can be seen that the H₂ mass fraction inside the black circle gradually decreases with time, which indicates that the mass fraction of H₂ gradually decreases during the process of convection from the high-pressure region shown in the black circle to the downstream while the pressure gradually increases in Figure 9, i.e., the process of pressure increase gradually reaches the peak accompanied by the combustion process.

In addition to the pressure oscillation analysis and local flow field characterization at the monitoring point, the dynamic mode decomposition of the spatial and temporal evolution fields of the whole combustor is also carried out in this paper. A total of 200 snapshots of the flow field is obtained by sampling the flow field data at 3 μs intervals from 1.8 to 2.4 ms with a sampling frequency of 333 kHz, and then the dynamic mode decomposition is carried out.

Figure 11 shows the DMD spectrum of the pressure field. The top 3 frequencies except 0 Hz are 2984 Hz, 62,180 Hz, and 110,562 Hz, respectively, and the corresponding DMD mode contours are provided in Figure 12. The phenomena and the physics behind the phenomena are explained below. The mode with a frequency of 2984 Hz (see Figure 12a) exhibits a significant wave-like structure, indicating that the oscillations of the shocks and expansion waves are more significant at this frequency. By comparing Figure 6a with Figure 12a, it can be seen that the shock waves generated at the leading edge of the strut and its reflected shock waves generated by the reflection from the upper and lower walls, the expansion wave fan generated at the bottom of the strut, and the expansion wave fan generated at the divergence angle of the upper wall do not appear in Figure 12a. These waves are more stable and do not have significant unsteady properties; The wave structure presented in Figure 12a is a wave system formed by waves originating from the wake of the strut and combustion shear layer and their reflected waves at the wall and shear layer, which are significantly unsteady due to the combined effects of wake instability, shear layer instability, and combustion. The modes at a frequency of 62,180 Hz (see Figure 12b) and 110,562 Hz (see Figure 12c) demonstrate prominent shear layer flow structures and wave system flow structures originating from the shear layer, indicating that shear layer instability and shear layer and wave interactions lead to high frequency oscillations in the flow.

The maximum peak frequency in the pressure oscillation spectrum at p9 is 48,333 Hz, and there is a frequency with a larger value in the pressure DMD spectrum (48,622 Hz) that is closer to this frequency. The modes with close frequencies in the DMD modes are morphologically closer, so the flow structure corresponding to the maximum peak frequency of the pressure oscillation at p9 can be understood by observing the DMD mode with a frequency of 48,622 Hz. Figure 13 shows this flow structure, which can be seen in the prominent shear layer structure and the wave system structure originating from this shear layer near the p9 point. As mentioned before, the maximum peak frequency of pressure oscillation at point p9 is exactly close to the frequency of the quasi-periodic peak of the pressure evolution curve at point p9 (see Figure 4a). The process of pressure increase gradually reaches the peak accompanied by the combustion process, so the large pressure oscillation at point p9 is related to shear layer instability, shear layer and wave interaction, and combustion.

3.2. Effect of Air and Fuel Conditions on Supersonic Combustion Instability

From the above analysis, it is clear that a typical significant pressure oscillation occurs in case 3. Comparing the numerical simulation results of all cases in Table 1, it is found that if the air parameter at the combustor inlet or the fuel parameter at the injector outlet is changed, it will be able to change the pressure oscillation characteristics.

The pressure evolution profiles at the monitoring points for case 1 are given in Figure 14. It can be seen that the pressure oscillation at most of the monitoring points for case 1 is slowed down compared with the pressure evolution at the monitoring points for case 3 in Figure 4. The pressure DMD spectrum can reflect the characteristics of the whole-field pressure evolution with time. The comparison of the pressure DMD spectrum of case 3 and case 1 is given in Figure 15. Although it is not shown in Figure 15 that the amplitude of each frequency in one case is lower than that in the other case, the black dots are higher than the red dots overall in this figure. That is, the pressure oscillation of case 1 is less than that of case 3 as a whole. Considering that the total air temperature at the combustor inlet in case 1 is lower than that at the combustor inlet in case 3, this may indicate that reducing the total air temperature at the combustor inlet can attenuate the pressure oscillations.

To further observe the effect of the total air temperature at the combustor inlet on the pressure oscillation, we compare the pressure evolution profiles (shown in Figure 16 and Figure 17) and pressure DMD spectra (shown in Figure 18) at the monitoring point for case 6 and case 4, and the pressure evolution profiles (shown in Figure 19 and Figure 5) and pressure DMD spectra (shown in Figure 20) at the monitoring point for case 7 and case 5. It can be seen that the level of pressure oscillation is less in case 4 than in case 6 and in case 5 than in case 7, which further demonstrates the effect of reducing the total temperature of the combustor inlet air to attenuate the pressure oscillation.

Comparing the pressure evolution profiles of the monitoring points for case 4, case 1, and case 5 (Figure 17, Figure 14, and Figure 5), it is observed that for most of the monitoring points, case 4 has the highest pressure oscillations, case 1 has the second highest pressure oscillations, and case 5 has the weakest pressure oscillations. The pressure DMD spectra represent the time evolution of the whole-field pressure. A comparison of the pressure DMD spectra for case 4, case 1, and case 5 is given in Figure 21. Considering that the fuel flow rate at the injector outlet decreases from case 4, to case 1, to case 5, this may indicate that reducing the fuel flow rate at the injector outlet can attenuate the pressure oscillations.

To further observe the effect of the fuel flow rate at the injector outlet on the pressure oscillation, the pressure evolution profiles (shown in Figure 16, Figure 4, and Figure 19) and the pressure DMD spectra (shown in Figure 22) at the monitoring points for case 6, case 3, and case 7 are compared. It can be seen that the pressure oscillation is the strongest in case 6, followed by case 3, and it is the weakest in case 7, which further demonstrates the effect of reducing the fuel flow rate at the injector outlet to attenuate the pressure oscillation.

4. Conclusions

In this paper, we conducted large eddy simulation and dynamic mode decomposition of the unsteady combustion process under different operating conditions in a strut-based scramjet combustor. The conclusions are as follows:

(1)

The pressure in the strut-based scramjet combustor shows significant oscillation characteristics when the air parameters at the combustor inlet and the fuel parameters at the injector outlet are under certain conditions, and these pressure oscillation situations correspond to supersonic combustion instability.

(2)

The pressure oscillations of the sampling points analyzed using FFT have broadband and two dominant frequency ranges. One dominant frequency range is below 20,000 Hz while the other is around 50,000 Hz.

(3)

The pressure oscillations of the whole combustor analyzed using DMD have multiple dominant frequencies, including relatively low frequency of 2984 Hz, high frequency of 62,180 Hz, and very high frequency of 110,562 Hz.

(4)

Large pressure oscillations in the strut-based scramjet combustor supersonic combustion instability are closely related to wake instability, shear layer instability, shear layer and wave interactions, and combustion.

(5)

Reducing the air total temperature at the combustor inlet can attenuate the pressure oscillations in supersonic combustion instability, and reducing the fuel flow rate at the injector outlet can also attenuate the pressure oscillations in supersonic combustion instability.