The governing equations for gaseous detonations used in present study are the one-dimensionally reactive Euler equationswhere v is the specific volume (v=1/ρ, ρ the density); U is the velocity of fluid particles in the laboratory coordinate system (X); p is the pressure; e is the specific internal energy; Y_i is the mass fraction of species i; Ω_i is the net production rate of species i; and N represents the number of species contained in the complex reactive system. In order to investigate the influence of the elementary reactions, a detailed chemical reaction mechanism of H₂/O₂/Ar (argon) [15] was adopted in this paper, which includes 9 species and 19 elementary reactions. It is assumed that thermally perfect gas in the reaction system is adopted, the equation of state for ideal gas is then used to describe change of the ther-modynamic state in the system.

To facilitate the linearization of the energy equation, original energy equation (Eq. (3)) needs to be expressed as follows, using the form derived by Kao [13]where a_f is the frozen sound speed; G is the Grüneisen coefficient; Y=[Y₁, Y₂, …, Y_N]^T is the vector of mass fractions of the species (Y_i); and Z(v, p, Y) is determined by the partial derivatives of the specific internal energy with respect to each species.

Using the coordinate attached to shock wave and placing the shock wave at the origin of the coordinate axis, the spatial position x and fluid particle velocity u in the new coordinate system can be expressed as follows:where u is the velocity in the shock coordinate system (x); and D is the speed of the detonation wave. In the new coordinate system, the transformation relation for the differential operator is given by

Thus, the governing Eqs. (1), (2), (4), and (5) can be written in the form of the following matrix.where the subscripts “,t” and “,x” represent the partial derivatives with respect to time and spatial distance, respectively; and the matrix z, c, and A can be represented as

In the present work, a method developed by Kabanov and Kasimov [14] is used to solve the linear stability problem of detonation waves. The method consists of two steps. The first step is to numerically solve the linearized equations of perturbations, in order to obtain the time series of speed perturbations of propagating detonation wave. In the second step, a dynamic mode decomposition (DMD) method [16] is used to process the results from the first step to extract the stability spectra of the perturbation. Unlike the normal mode method, this method does not require any special treatment of the sonic plane within the reaction zone of detonation wave, thus making it easier to extend to more complex chemical reactions and equations of state.

Before proceeding with the computations mentioned above, the form of the linearized governing equation of perturbations should be given first. To achieve this, we can decompose the quantities to be solved in governing Eq. (9) into steady-state and perturbation components, resulting in the following expressionwhere the superscript “*” denotes the steady-state solutions; the superscript “'” denotes the perturbation solutions; z′=[v′, U′, p′, Y₁′, …, Y_N′]^T; and ψ′ is the speed perturbation of the propagating detonation wave, 0< ε << 1.

By introducing Eqs. (11)-(12) into Eq. (9) and omitting the high-order nonlinear term, the following linearized governing equations of perturbations are obtainedwhere A^* denotes the steady-state solution of A; and the expression of B^* is as follows:where the subscripts “,x”, “,v” and “,p” represent the partial derivatives of the variables with respect to spatial distance, specific volume and pressure, while the subscripts “,Y₁”, …,“,Y_N” represent the derivatives with respect to the mass fraction of species. Note that the Eqs. (13) and (14) are the extension of original perturbation linear equations [14] in complex reactive system with arbitrary number of species, thus those partial derivatives associated with species i (i=1,…, N), such as, Z_,Yi, Y_i,x, Ω_i,v, Ω_i,p, Ω_i,Yi, Z_,Yi and Y_i,x in matrix B* of Eq. (14) and ∂e/∂Yi in Eq. (6) are determined or calculated independently in present study.

Since Eq. (13) contains N+3 governing equations but has N+4 unknowns (ψ' lacks evolution equations), the equations are not closed. Therefore, it is necessary to supplement the evolution equation for ψ'. For this reason, we adopt the equation derived in the Ref. [17]where the subscript “s” is the state just behind the shock wave; the subscript “a” is the state before the shock wave; and γ is the adiabatic index.

Using Eqs. (11) and (12), we linearize the evolution equation for the detonation speed D given by Eq. (15) and obtain the evolution equation for speed perturbation ψ'

In the previous section, we have linearized the governing equations and obtained the partial differential equations (PDEs) regarding z' (Eq. (13)) and ψ' (Eq. (17)). To solve the PDEs, we first define the length of the one-dimensional computational domain to be 30 times the induction zone length l_ind, ensuring that the reaction zone falls completely within the domain. Here the l_ind is defined as the distance from the location of the leading shock front of the detonation wave to the position of maximum heat release rate in the reaction zone. Next, we discretize the computational domain using a uniform grid {x_i}, i=0, …, M.

Appropriate boundary conditions are specified at the starting and ending positions of the computational domain. The starting position corresponds to x_M=0, while the ending position corresponds to x₀= -30 l_ind. At x_M=0, the boundary condition just behind the leading shock front of the detonation wave can be described using the Rankine-Hugoniot relations, which provide a conservation relationship for the state variables before and after the shock wave, as follows:

Linearize Eq. (20) to obtain the final form representation at the boundary as follows:

The boundary condition at x₀ = -30 l_ind is set as a zero gradient boundary condition, which avoids the complex radiation boundary conditions used in the normal mode method. Thus, the computation is greatly simplified, which is very favorable for the computation of complex chemical reaction models. We also need to provide initial conditions for z' and ψ'. Firstly, the initial speed perturbation is specified as ψ'|_t=0=0.01. The Rankine-Hugoniot conditions are then used to find v'_s, U'_s, and p'_s from Ψ'. Subsequently, initial conditions for perturbations over the whole grid are computed by usingwhere v^*(x), U^*(x), p^*(x) and are the steady-state solutions of specific volume, the velocity of fluid particles, pressure, and species, respectively, as functions of x.

The partial derivatives of the thermodynamic parameters and kinetic parameters that appeared in above equations and matrices can be determined by using Cantera [18] and SDtoolbox [19].

Given the boundary conditions of the linearized Eqs. (13) and (17), a finite difference scheme [14] is adopted to approximate the spatial derivatives on each grid point at the current moment, thus the linearized PDEs are transformed into the ordinary differential equations (ODEs) only with respect to timewhere and are the numerical approximations of L(z^*, z', z'_,x) = A^*(z^*)z'_,x + B^*(z^*)z'-z^*_,xΨ' from Eq. (13) and s'(z^*, z', z'_x) from Eq. (17). The time series of ψ' controlled by the linearized equation can be obtained by integrating the ODEs with an ODE solver VODE [20].

Once the time series of Ψ' are obtained, a DMD method [16] is employed to extract the eigenvalues α_i for mode i (i=0, …, n) from the time series of dimensionless form of Ψ' = /Ψ' (a_f)_a, here a_f is the frozen sonic speed. It should be noted in particular that the eigenvalues α_i used hereinafter in the present work are all in dimensionless form that characters the perturbation mode in detonation system. α_i is a complex number, where its real part α_re represents the growth rate of the perturbation, and its imaginary part α_im represents the frequency of the perturbation. Each eigen-value represents a mode of perturbation, which is sorted by perturbation frequency from low to high as mode 0, mode 1, and so on. If the growth rate of all modes obtained by the DMD is less than 0, the perturbation is decayed, indicating that the detonation wave is stable in this case, otherwise the detonation wave is unstable. The details of computational procedure of the DMD method can be found in the Ref. [16].

To determine the grid resolution in the finite difference approximations, a grid independence test in the computations of the linear stability analysis is carried out, for the case of detonation of stoichiometric H₂/O₂ mixture with 40 vol% Ar dilution ratio at an initial pressure of 1 atm (1 atm=1.01325 × 10⁵ Pa). Table 1 shows the eigenvalues of the first four perturbation modes at three grid resolutions. It can be found that as the grid number per l_ind increases, the real and imaginary parts of the eigenvalue α_i only show the slight changes, indicating that the value of α_i is independent of the grid resolution under the given three grid resolutions. Considering the accuracy and cost of the computational results, the grid resolution of 160 pts/l_ind (pts/l_ind is the number of grid points in the length of each induction zone) is adopted in the subsequent linear stability computations.

For the global reaction systems, the activation energy, adiabatic index, reaction heat release and detonation overdrive can usually be regarded as the parameters affecting the stability of gaseous detonation wave propagation [22]. However, for the complex reaction systems, most of the above parameters can neither be easily determined nor have practical significance. In general, the apparent activation energy, the thermal effect and the adiabatic index in complex reaction systems depend on the composition of reactive gases. Therefore, the selection of gas composition as a controlling parameter can not only investigate the stability of detonation, but also is significant in selecting the appropriate gas in industrial safety applications. In addition, environmental conditions (such as initial pressure and temperature) also have a significant effect on the stability of gaseous detonation waves. Therefore, for the H₂/O₂/Ar mixture in the present study, we choose the initial pressure (p₀) and Ar dilution ratio (volume percentage) as the controlling parameters to study their influence on the detonation stability.

In Sect. 2.3, we have shown that the linear stability spectrum is extracted from the time series of by the DMD method. Therefore, before exploring the effects of controlling parameters on the stability spectrum, the time histories of are first analyzed. Figure 3 shows that variations of with time at the two typical initial pressures p₀ of 0.3 and 0.8 atm, respectively, for stoichiometric H₂/O₂ mixture diluted by 30% Ar. Figure 3a shows a stable detonation case, in which the speed perturbation of the propagating detonation wave shows a two-mode oscillation pattern: A low frequency mode superposed by a high frequency mode. A dashed curve is plotted as an envelope curve of the speed perturbation pattern. Obviously, the low frequency perturbation exhibits the exponential decay , indicating that the detonation process is immune to external infinitesimal perturbations, especially for the low-frequency component. Such results also suggest that the real part α_re of the eigenvalues extracted from time series of ψ' is equal to the coefficient on the exponential term of the envelope curve, at which the coefficient (-0.089) is negative in this case. On the contrary, the envelope curve for speed perturbation pattern in the case with elevated initial pressure (p₀=0.8 atm), as shown in Fig. 3b, shows the unstable detonations, in which the perturbation exponentially grows with a positive coefficient on the exponential term of the envelope curve (0.177). The result in Fig. 3b indicates that the detonation process cannot resist external infinitesimal perturbations and gradually becomes unstable when initial pressure rises.

Given the time series of speed perturbation of detonation, we can determine the growth rate and the frequency of all perturbation modes at the specified controlling parameters of detonation. Figure 4 gives the stability spectra of detonation at several typical initial pressures p₀ for the stoichiometric H₂/O₂ mixture with 30% Ar dilution ratio. Here the stability spectra are exhibited by perturbation modes with corresponding complex α_i, in which the horizontal axis α_re represents the real part of α_i corresponded to the perturbation growth rate, while the vertical axis α_im represents the imaginary part of α_i corresponded to the perturbation frequency. For each initial pressure, the symbols in the figure represent several modes of linear stability of detonation. The mode with the smallest perturbation frequency is defined as mode 0. With the increase of the perturbation frequency, the corresponding modes can be sequentially referred to as mode 1, mode 2, etc. It can be seen from Fig. 4 that there are two unstable modes (mode 0 and mode 1) in the case of p₀=0.8 atm, at which they are located on the right of the vertical dashed line in Fig. 4 (α_re > 0). The perturbation growth rate (α_re) and the perturbation frequency (α_im) of mode 0 are 5.846 × 10^-2 and 0.591, respectively, while the α_re and α_im of mode 1 are 0.177 and 4.195, respectively. Note that the value of α_re for mode 1 is just the coefficient on the exponential term of the envelope in Fig. 3b for the same case. When the initial pressure is reduced to p₀=0.6 atm, the α_re of mode 0 decreases from 5.846 × 10^-2 at 0.8 atm to 8.881 × 10^-3, and the α_re of mode 1 decreases from 0.177 to 0.109, indicating that unstable detonation is diminished with decreasing initial pressure. When the initial pressure is reduced to 0.4 atm, only mode 1 remains unstable, with its growth rate further decreasing to 1.819 × 10^-2, suggesting a further reduction in instability. When the initial pressure is further reduced (p₀=0.2 atm), the unstable mode will no longer appear in the stability spectra and the propagation of detonation wave is linearly stable. In addition, it is also found from the results in Fig. 4 that the perturbation frequency for each mode remains basically unchanged when the initial pressure is changed.

The results of Fig. 4 show that for 30% Ar diluted H₂/O₂ mixture, the decrease of initial pressure is beneficial to stabilizing the detonation. To compare with the nonlinear stability results of detonation wave under the same parameter conditions, we numerically simulate the one-dimensionally nonlinear oscillation characteristics of detonation wave under four initial pressures of 0.2, 0.4, 0.6 and 0.8 atm for the same mixture, as shown in Fig. 5. In the figure, the nonlinear oscillation of the detonation wave is characterized by recording the pressure peak time history of the leading shock front of detonation wave. The results in Fig. 5a show that when the initial pressure is low (p₀=0.2 atm), the pressure peak quickly tends to be constant after an initial large oscillation, indicating that the propagation of the detonation wave under this p₀ is nonlinearly stable, which is consistent with the results of linear stability computed under the same parameter conditions in Fig. 4. In Fig. 5b (p₀=0.4 atm), the pressure peak of the detonation wave firstly shows the pattern of low-frequency, low-amplitude oscillation after experiencing the initial large-amplitude oscillation, and then gradually transforms into an oscillation mode superposed by high-frequency and low-frequency perturbations. The amplitude of the superposed mode decreases slowly with time. The high-frequency, low-amplitude component in this oscillation mode represents a critical state of detonation instability, which will be discussed further in the next section. With the further increase of the initial pressure, the pressure peak time history of the detonation wave finally presents a single oscillation mode. For the case of p₀=0.6 atm, the amplitude of the single oscillation is about 12 p₀ (Fig. 5c), while for the case of p₀= 0.8 atm, the amplitude of the oscillation increases to about 17 p₀, indicating that with the increase of the initial pressure, the pressure peak of the detonation wave becomes more unstable.

Comparing the results in Fig. 4 with those in Fig. 5, it can be found that the results from linear stability analysis are consistent with those from nonlinear stability simulations. In Fig. 4, the linear stability results indicate that at p₀=0.4 atm, there exists only one unstable mode (mode 1) with a relatively small perturbation growth rate (α_re=1.819 × 10^-2), suggesting that the propagation of the detonation wave is approaching a stable state at this initial pressure. Furthermore, the results in Fig. 5b demonstrate that under the same initial pressure, the time history of pressure peak of the detonation wave reaches a critical state characterized by high-frequency, low-amplitude fluctuations, which qualitatively aligns with the findings from the linear stability analysis. As the initial pressure increases, both linear and nonlinear results indicate that the instability of the detonation wave gradually intensifies.

To investigate the impact of Ar dilution ratio on the stability of detonation, Fig. 6 illustrates the linear stability results of the detonation at various Ar dilution ratios and a fixed initial pressure of p₀=1 atm. The results show that at the 20% Ar dilution ratio, there are three unstable modes (mode 0, mode 1, and mode 2), with the most unstable mode1 exhibiting a perturbation growth rate of 0.364. As the Ar dilution ratio increases, the perturbation growth rate of the most unstable mode 1 gradually decreases, indicating that higher Ar dilution ratio leads to more stable detonation wave. Additionally, the results in Fig. 6 reveal that modes with higher perturbation frequencies experience a more rapid decline in growth rate with an increasing Ar dilution ratio. However, the change of Ar dilution ratio does not significantly alter the perturbation frequencies of the same modes, suggesting that within a certain range, the periodicity of the perturbation remains relatively constant. This conclusion is consistent with the experimental results of Mcvey and Toong [23] and Lee's linear stability analysis [5].

Figure 7 presents the results of the one-dimensionally nonlinear oscillation of the detonation wave under the same parameter conditions as in Fig. 6. Figure 7a is the pattern of pressure peak at 20% Ar dilution ratio. One can see that under this Ar dilution ratio, the pressure peak shows a trend of increasing amplitude with time after experiencing a series of irregular oscillations. The raised amplitude makes the pressure at the trough position of oscillation pattern continue to decrease, so that detonation may be extinguished due to the excessive attenuation of shock wave intensity in the later stage, indicating that the detonation is extremely unstable in this case. When the Ar dilution ratio increases to 30%, the time history of pressure peak eventually exhibits the characteristics of dual-mode oscillation with low-frequency, high-amplitude pattern, as shown in Fig. 7b. Although this oscillation is still unstable, it is weaker than that in the case with 20% Ar dilution ratio. When the Ar dilution ratio is further increased to 40%, the pressure peak is finally presented as a single-mode oscillation with a low-frequency, medium-amplitude pattern. The above results clearly show that with the increase of Ar dilution ratio, the detonation of H₂/O₂/Ar mixture is gradually stable, which is qualitatively consistent with the results by linear stability analysis in Fig. 6.

Both the results from linear stability computations and one-dimensionally nonlinear numerical simulations in previous section demonstrate that the initial pressure and the Ar dilution ratio in the complex reaction system significantly influence the stability of detonation waves in H₂/O₂/Ar mixture. To further quantitatively compare the results obtained from the two methods, we plot the critical stability curves for mode 0 and mode 1 on the phase plane of p₀-Ar% based on the results from linear stability analysis. The selected range for p₀ is from 0 to 1.4 atm, while the Ar dilution ratio ranges from 0 to 45%. To determine the stability boundaries for a given mode, we first specify an Ar dilution ratio and then identify the critical value of p₀ at which the growth rate of that mode approaches to zero, allowing for a relative error of 10^-3.

Figure 8 presents two critical stability curves for mode 0 and mode 1 in the p₀-Ar% plane, determined by using the linear stability method. It is evident that the stability curve for mode 1 does not intersect with the curve for mode 0 within the given p₀-Ar% plane, with the stability curve for mode 1 located below that for mode 0. Linear stability results indicate that only below the stability curve for mode 1 do all modes exhibit a negative growth rate, meaning that in the p₀-Ar% plane, the stability of the detonation wave is governed by the critical stability curve of mode 1. Therefore, we call this curve the neutral stability curve of mode 1. On the other hand, linear stability results also show that the perturbations above the critical stability curve of mode 0 exhibit the positive growth rate for low-frequency mode 0, while the perturbations below the curve of mode 0 but above the curve of mode 1 display the negative growth rate for low-frequency mode 0 and the positive growth rate for high-frequency mode.

This suggests that the critical stability curve of mode 0 separates the high-frequency perturbation mode from the low-frequency perturbation mode. Thus, we can also call the critical stability curve of mode 0 the perturbation frequency transition curve. For analytical convenience, we define the stable region below the neutral curve of mode 1 as region I. Above the neutral curve of mode 1, the frequency transition curve of mode 0 further divides this unstable area into two regions: The area above the curve of mode 0 is defined as region III, while the area between both curves of modes 0 and 1 is defined as region II. Clearly, region III represents a linearly unstable area, while region II serves as a transitional area situated between the stable region I and the unstable region III.

To further compare with the nonlinear oscillation behavior of detonation, we selected several suitable parameter pairs (p₀, Ar%) within the p₀-Ar% plane shown in Fig. 8 and conducted numerical simulations of one-dimensionally nonlinear oscillatory detonation. The numerical results show that all cases exhibiting pressure oscillation over time (marked as “×” and “□”) are located above the neutral stability curve of mode 1, while all cases in which the pressure oscillations ultimately decay (marked as “■”) fall below this curve. These results suggest that in complex reactive systems considering detailed elementary chemical mechanisms, the results of linear stability computations well agree with those of one-dimensionally nonlinear oscillation of pressure peak of detonation wave. Ng et al. [24] found that the critical activation energy for one-dimensionally non-linear instability of detonations obtained from numerical simulations of a single-step reaction corresponds well with the linear stability results under the same conditions. Weng and Mével [9] studied the linear stability and one-dimensionally nonlinear stability results using a single-step chemical reaction model that accounts for molecular volume effects, and the good agreement between the linear and nonlinear results was obtained. Sharpe and Falle [25] addressed the stability of pathological detonations in a two-step consecutive reaction model consisting of an exothermic reaction followed by an endothermic reaction, identifying strong correlations between linear stability theory and non-linear oscillation characteristics. In the present study, the results from Fig. 8 also suggest that the linear stability analysis results for complex reaction systems are consistent with those from one-dimensionally nonlinear numerical simulations.

In Fig. 8, region II between the curves of mode 0 and mode 1 is particularly interesting. The parameters selected for the one-dimensionally nonlinear oscillation case shown in Fig. 5b are located in this region. From the results of Fig. 5b, the peak pressure of the detonation wave shows a significant high-frequency characteristic, which is very different from the low-frequency instability in region III (such as cases in Fig. 5c and d). To explore the nonlinear oscillation characteristics of detonation wave in this region, we select five Ar ratios of 15%, 20%, 25%, 30%, and 40% as the Ar dilution parameters. For each Ar dilution ratio, we also select a high initial pressure near the curve of mode 0 and a low initial pressure near the curve of mode 1 as the p₀ parameter, respectively. Such selections produce 10 computational cases that are marked with “□” in region II of Fig. 8.

Figure 9 presents the numerical simulation results for five cases near the lower boundary (the neutral stability curve of mode 1) in region II. When the Ar dilution ratio is 15%, the pattern of pressure peak over time exhibits a clear high-frequency, low-amplitude oscillation pattern (Fig. 9a) after several large initial oscillations and subsequent decay. At 20% and 25% of Ar dilution ratios, the pressure peak first exhibits the significant initial oscillations. These oscillations subsequently transform into a low-frequency, low-amplitude pattern. Ultimately, they shift back to a high-frequency, low-amplitude mode (Fig. 9b and c). As the Ar dilution ratio increases to 30%, the decay of the low-frequency, low-amplitude mode slows down but the growth of the high-frequency, low-amplitude mode becomes more pronounced. Thus, a superposed oscillation mode with low and high frequencies is shown in late stages of time history of pressure peak, as shown in (Fig. 9d). When the Ar dilution ratio is further increased to 40%, the long-term behavior of the pressure peak exhibits a unique oscillation pattern in which the high-frequency perturbation is overlaid by the low-frequency perturbation (Fig. 9e). A locally enlarged pattern of oscillation mode in Fig. 9e is shown in Fig. 9f, the enlarged figure clearly illustrates the pattern from the superposition of these two oscillation modes, which is similar to Fig. 7 in one-dimensional pulsating detonation of two-step chain-branch reaction model by Leung et al. [26]. These results indicate that within the selected range of Ar dilution ratio in present study, region II displays high-frequency, low-amplitude oscillation characteristics near the low boundary of stability. This implies that the results from the linear stability analysis also contain components of high-frequency perturbations near the low boundary of region II. Furthermore, it is found that as Ar dilution ratio increases, the amplitude of high-frequency oscillations gradually decreases, and the transition of pressure peak from low-frequency mode to high-frequency mode requires a longer time. This is attributed to the decreasing difference in growth rates between high-frequency and low-frequency modes when Ar dilution ratio rises (see Fig. 6), leading to a longer period for their growth or decay.

In addition to examining the stability characteristics of detonation near the lower boundary of region II, we also investigate the properties of one-dimensionally nonlinear oscillations near the upper boundary (the perturbation frequency transition curve of mode 0) within this region. Figure 10 shows the time histories of pressure peak of detonation wave under five different Ar dilution ratios near the upper boundary of region II. At Ar dilution ratios of 15% and 20%, the pressure peaks experience significant initial oscillations and then form a superposed pattern of high-frequency and low-frequency perturbations. Over time, this pattern gradually evolves into a low-frequency, high-amplitude oscillation mode (Fig. 10a and b). However, when the Ar dilution ratio increases to 25% or higher (Fig. 10c-e), the initial high-frequency perturbations become less pronounced, and the pressure peaks primarily exhibit low-frequency, high-amplitude oscillation behavior over time.

The results in Fig. 10 also reveal that in all cases, the pressure peaks characterized by low-frequency, high-amplitude pattern gradually decay over time. Especially, this decay becomes more pronounced at higher Ar dilution ratios (as seen in Fig. 10e). At lower Ar dilution ratios, the decay of pressure peaks is weaker, with even slightly initial growth observed (as shown in Fig. 10c). To further investigate these oscillatory behaviors of pressure peak, we perform a fast Fourier transformation on the time series of pressure peaks for the five cases in Fig. 10. From the transformed results, we extract the perturbation frequencies with significant oscillation amplitudes, denoted as f_num. Additionally, we also analyze the perturbation growth rates α_re and perturbation frequencies α_im for modes 0 and 1 obtained from the linear stability results of these five cases. These results are summarized in Table 3. Note that the linear perturbation frequencies α_im and the nonlinear perturbation frequencies f_num in the table are comparable, and their differences can be expressed by using the following relative error εwhere the value of ε represents the difference of perturbation frequencies between the linear and nonlinear results.

The results in Table 3 show that the results of the five cases at the locations near the perturbation frequency transition curve of mode 0 in region II are similar to each other. Therefore, the results of the case 3 with Ar dilution ratio of 25% are taken as an example to illustrate. In this case, the linear stability analysis shows that the perturbation growth rate of mode 0 (α_re= -9.48 × 10^-3) is negative, and its perturbation frequency is small (α_im=0.620), indicating that mode 0 is a decaying low-frequency perturbation. The perturbation growth rate of mode 1 (α_re=0.10752) is positive, and its perturbation frequency is large (α_im=4.221), indicating that mode 1 is a growing high-frequency perturbation. On the other hand, the nonlinear numerical simulation for the same case also produces two perturbation modes. The frequency of the low-frequency perturbation mode is in good agreement with the result of the low-frequency perturbation of linear stability (ε=2.26%), while the high-frequency perturbation result is also in reasonable agreement with the high-frequency perturbation result of linear stability (ε=9.88%). Weng and Mével [9] also compared the perturbation frequencies between the linear and nonlinear results in the detonation stability analysis with single-step reaction considering molecular effect. The range of relative error (ε) of their results is between 0.038% and 5.0% by using the same Eq. (25). This error range is consistent with the error range of the low-frequency perturbations we predicted (1.12-3.74%) and is on the same order of magnitude as that of the high-frequency perturbations we obtained (7.07-10.07%), which can demonstrate the validity of our calculations.

The results of Table 3 demonstrate that oscillation pattern of the pressure peak in Fig. 10 still contains high-frequency perturbation component, and the oscillation pattern is the result of the superposition between high-frequency perturbation and low-frequency perturbation. From the five cases in Table 3, it can be further observed that as Ar dilution ratio increases, both the decay rate of mode 0 and the growth rate of mode 1 show a slower trend. Therefore, the pattern of pressure peak obtained by nonlinear simulations in Fig. 10 is related to the growth and decay of the amplitude of the two modes. For case 1, although the low-frequency, high-amplitude perturbation decays (α_re= -0.0174), the super-imposed perturbation decays slowly, due to the rapid growth of high-frequency, low-amplitude perturbation (α_re= 0.13847) (see Fig. 10a). For case 5, due to the slow growth of high-frequency, low-amplitude perturbation (α_re= 0.04337), the decay of the superposed perturbation becomes more obvious (see Fig. 10e). For case 3, with moderate growth of high-frequency perturbations and the relatively moderate decay of low-frequency perturbations, the overall pressure peak first increases to the maximum value and then gradually decays (see Fig. 10c).

Although the results in Table 3 show that the low-frequency perturbation of mode 0 of these five cases decays, their decay is very slow because these cases are located near the frequency transition curve of mode 0. In the numerical simulation of the longitudinal oscillation of the pressure peak in Fig. 10, the complete decay of the low-frequency oscillation would take a very long time (Fig. 10 does not show this). However, it can be expected that the high-frequency oscillation will eventually appear, if decay time of the low-frequency oscillation is long enough. The above analysis confirms that the curve of mode 0 is the critical curve separating the low-frequency perturbations from the high-frequency perturbations. This conclusion is consistent with the results obtained by Short and Wang [27] using the three-step chain branching reaction model.

Figures 9 and 10 show the difference in longitudinal nonlinear oscillation of the pressure peak of detonation wave between near the lower boundary of region II (the neutral curve of mode 1) and near the upper boundary of region II (the perturbation frequency transition curve of mode 0). For these cases, the former shows a primary high-frequency, low-amplitude oscillation mode, while the latter mainly shows a low-frequency, high-amplitude oscillation mode. To further illustrate the continuous transition from the high-frequency mode to the low-frequency mode in the region, we also add a case with an initial pressure p₀=0.47 atm and a 30% Ar dilution ratio, which is located in the middle of the region II (see the case marked as “Δ” in Fig. 8). Figure 11 shows the time history of pressure peak and its locally enlarged results for this case. This result clearly demonstrates the interaction between high-frequency and low-frequency perturbations, which is similar to that in Fig. 9f. Both results in Figs. 11b and 9f suggest a continue transition of perturbation from low-frequency, high-amplitude oscillation to high-frequency, low-amplitude oscillation in region II, when initial pressure gradually decreases or Ar dilution ratio gradually increases. These results obtained by numerical simulations verify the accuracy and effectiveness of the linear stability theory in predicting the stability of detonation waves in complex chemical reaction systems. By virtue of the low-cost characteristic of linear stability calculations, it can greatly save resources when predicting the stability of detonation waves and selecting appropriate initial conditions for detonation stability research in complex reactive systems. However, our predictions can only cover the one-dimensional perturbation case, and more research is still needed for linear stability analysis on multi-dimensional perturbations in complex reactive systems.