2. Governing equations and non-uniform conditions

Consider a one-dimensional planar detonation wave propagating in a gaseous mixture in the positive $x$ direction (figure 1). The detonation dynamics is governed by the reactive Euler equations (Fickett & Davis Reference Fickett and Davis1979; Kasimov & Goldin Reference Kasimov and Goldin2021):

(2.1)\begin{gather} \rho_{t}+\left(\rho u\right)_{x}=0, \end{gather}

(2.2)\begin{gather}\left(\rho u\right)_{t}+(p+\rho u^{2})_{x}=0, \end{gather}

(2.3)\begin{gather}\left(\rho e\right)_{t}+\left(\rho u\left(e+p/\rho\right)\right)_{x}=0, \end{gather}

(2.4)\begin{gather}\left(\rho\lambda\right)_{t}+\left(\rho u\lambda\right)_{x}=\rho\omega, \end{gather}

that consist of equations of continuity (2.1), momentum (2.2), energy (2.3) and reaction progress (2.4), Here, subscripts $t$ and $x$ denote partial derivatives with respect to time $t$ and space $x$. Variables $\rho$, $u$, $p$, $e$, $\lambda \in [0,1]$ are the density, particle velocity, pressure, total specific energy and reaction progress variable, respectively. The reactive medium is assumed to be an ideal gas described by the equation of state $p=\rho RT$ with gas constant $R$ and temperature $T$. The gas undergoes a single-step irreversible chemical reaction obeying the Arrhenius rate, $\omega =K(1-\lambda )\exp (-E/(RT))$, with rate constant $K$ and activation energy $E$, for reactants $\lambda =0$ and for products $\lambda =1$. The total specific energy $e$ of the gas is given by $e=(\gamma -1)^{-1}p/\rho + u^2/2 -\lambda Q,$ where $\gamma$ is the ratio of specific heats, and $Q$ is the specific heat release.

Figure 1. The schematics of the density profile for detonation in a non-uniform medium.

Consistent with the existing analyses of realistic scenarios in RDE (e.g. Fujii et al. Reference Fujii, Kumazawa, Matsuo, Nakagami, Matsuoka and Kasahara2017; Gaillard, Davidenko & Dupoirieux Reference Gaillard, Davidenko and Dupoirieux2017) wherein injector systems lead to nearly periodic variations of temperature and fuel concentration upstream of the rotating detonation front, we model the upstream conditions (at $x>0$) as

(2.5)\begin{gather} T_{{ {a}}}(x)=T_{{ {a}}}+A(1+\sin kx), \end{gather}

(2.6)\begin{gather}\lambda_{{ {a}}}(x)=A(1+\sin kx). \end{gather}

Here $T_{{ {a}}}$ is the minimal temperature of the fresh mixture, $0\leq A\leq 0.5$ and $k$ are the amplitude and wavenumber of the variations, respectively. We further assume zero fluid velocity (in the laboratory frame) and constant pressure $p_{{ a}}$ throughout the fresh mixture. The variables are rescaled with respect to $p_{{ a}}$ for pressure, $\rho _{{ a}}=p_{{ a}}/(RT_{a})$ for density, and $\sqrt {p_{{ a}}/\rho _{{ a}}}$ for velocity. The length scale is the standard half-reaction length, $l_{1/2}$, i.e. the distance from the shock over which half of the chemical energy is released; the time scale is then $l_{1/2}/\sqrt {p_{{ a}}/\rho _{{ a}}}$.

Governing equations (2.1)–(2.4) are solved numerically by a high-order shock fitting method based on the conservative WENO5M scheme for the spatial discretization and total variation diminishing Runge–Kutta method for time integration as described in Henrick et al. (Reference Henrick, Aslam and Powers2006) and Kasimov & Goldin (Reference Kasimov and Goldin2021). The problem is posed in the moving shock-attached frame of reference. The location of the leading shock wave is always fixed at the origin of the coordinate system. Such an approach allows for high accuracy because the shock is not smeared in contrast to shock capturing methods. Since detonation velocity $D(t)$ explicitly appears in the governing equations, the system requires an additional shock-change equation (see Kasimov & Goldin Reference Kasimov and Goldin2021):

(2.7)\begin{equation} \frac{\mathrm{d}D}{\mathrm{d}t} ={-} \left(\frac{\partial\left( \rho_{{s}} u_{{s}}\right)}{\partial D}\right)^{{-}1} \left. \left(\frac{\partial \left(\rho u \left({u-D}\right)+p\right)}{\partial x} + D \frac{\partial \left(\rho_{{s}} u_{{s}}\right)}{\partial x}\right) \right\vert_{x=0}, \end{equation}

that is integrated numerically together with the reactive Euler equations. Here, subscript ${s}$ denotes the variables immediately after the shock jump. To determine these variables in terms of $D$ and the upstream conditions (2.5)–(2.6), the Rankine–Hugoniot conditions are used:

(2.8)\begin{gather} \rho_{{s}} \left(D-u_{{s}}\right) = \rho_{{a}} \left(D-u_{{a}}\right), \end{gather}

(2.9)\begin{gather}p_{{s}} - p_{{a}} = \left(\rho_{{a}} \left(D - u_{{a}}\right)\right)^{2} \left(\frac 1\rho_{{a}} -\frac 1\rho_{{s}}\right), \end{gather}

(2.10)\begin{gather}e_{{s}} - e_{{a}} = \frac 12\left(p_{{s}} + p_{{a}}\right) \left(\frac 1\rho_{{a}} -\frac 1\rho_{{s}}\right), \end{gather}

(2.11)\begin{gather}\lambda_{{s}} = \lambda_{{a}}. \end{gather}

Numerical integration of $\mathrm {d} x_{{s}}/\mathrm {d}t = D(t)$ gives the shock position $x_{{s}}$ in the laboratory frame of reference. Shock position $x_{{s}}$ is needed to determine the upstream state just ahead of the shock. We solve the governing equations in the domain of length $L=50$ behind the shock wave with $N=1000$ grid points up to time $t=3000$. Then, we extract detonation velocity $D(t)$ and apply analysis of time series and various tools of dynamical systems theory to study the detonation behaviour. The overall third order of accuracy of the numerical method is verified on a sequence of spatial grids with decreasing step size. The accurate numerical solutions allow us to compute the rather subtle effects reported in this paper.

3. Frequency locking

The key property of gaseous detonation that is important for the purpose of the present analysis is that it propagates in a self-sustained oscillatory manner if the activation energy $E$ exceeds a critical value $E_{{ c}}$ with all other parameters fixed (Fickett & Davis Reference Fickett and Davis1979). From the point of view of dynamical systems theory, the detonation can be looked at as a self-sustained nonlinear oscillator. It is therefore reasonable to expect that when subjected to an external periodic influence, the detonation will respond in a manner similar to that of well-known nonlinear oscillators such as, for example, the van der Pol oscillator (Parlitz & Lauterborn Reference Parlitz and Lauterborn1987; Pikovsky, Rosenblum & Kurths Reference Pikovsky, Rosenblum and Kurths2001; Manneville Reference Manneville2004; Strogatz Reference Strogatz2018). In particular, the periodic forcing is expected to lead to resonances and mode locking, the latter closely connected to synchronization. This is indeed what we have found as described next.

The following results are obtained for fixed parameters: $\gamma =1.2$, $Q=50$, $E=26$. Note that $E>E_{{c}}=25.26$. In this case, the steady-state self-sustained CJ detonation velocity is $D_{{ CJ}}=\sqrt {\gamma +(\gamma ^{2}-1)Q/2}+\sqrt {(\gamma ^{2}-1)Q/2}\approx 6.81.$ We investigate the long-time dynamics of the detonation wave at fixed $k=0.1$ and varying $A$. The influence of variations in wavenumber $k$ was previously studied in Kasimov & Goldin (Reference Kasimov and Goldin2021), where resonance and locking effects were discussed. When $A=0$, the time series of detonation velocity $D(t)$ shows a period-1 limit cycle with fundamental frequency $\nu _{{ 0}}=0.085$. These oscillations of $D(t)$ are shown in figure 2(a) with the power spectrum showing its maximum at natural frequency $\nu _0$ in figure 2(d). The addition of upstream variations with amplitude $A=0.025$ leads to the quasiperiodic regime shown in figure 2(b,e) which is seen to contain the linear combinations of two dominant frequencies, $\nu _{{ D}} = 0.071$ as a modified natural frequency and $\nu _{{ f}} = 0.107$ as arising due to the presence of forcing. The further increase of $A$ up to $0.035$ shifts $\nu _{{ D}}$ to $0.054$ which is close to $\nu _{{ f}}/2$ and the oscillations exhibit a period-2 limit cycle with dominant frequency $\nu _{{ D}}$, shown in figure 2(c,f). A more detailed dependence of detonation behaviour on $A$ is shown in figure 3 that depicts the heat map of the power spectra of $D(t)$ for a range of $A$ from $0$ to $0.045$. The frequency $\nu _{{ D}}$ of $D(t)$ is seen to decrease somewhat from $\nu _{{ 0}}$ as the amplitude $A$ increases. The almost vertical line near $\nu _{{ f}}\approx 0.107$ appears in the spectra because it is the time frequency with which the detonation shock propagates through the upstream variations of constant spatial frequency $k/(2{\rm \pi})$. These latter values are related as $\nu _{{ f}}=\bar {D}k/(2{\rm \pi})$, where $\bar {D}$ is the mean detonation velocity.

Figure 2. The normalized detonation speed and its spectra for different amplitudes of the upstream state: (a,d) $A=0$ – period-1 oscillations in a uniform medium; (b,e) $A=0.025$ – a quasiperiodic regime; (c,f) $A\!=\!0.035$ – period-2 synchronization.

Figure 3. The logarithm of the power spectral density for the time series of detonation velocity $D(t)$ for a range of values of amplitude $A$ from $0$ to $0.045$ and fixed forcing wavenumber $k=0.1$. The frequencies $\nu _{0}=0.085$ and $\nu _{f}=0.107$ indicated on the $\nu$ axis are due to the intrinsic and forcing oscillations, respectively.

The ratio $\nu _{{ D}}/\nu _{{ f}}$ determines the possible modes of detonation behaviour. Quasiperiodic regimes develop for incommensurate $\nu _{{ f}}$ and $\nu _{{ D}}$ and produce spectra containing various linear combinations of these two main frequencies (see figure 3 for $A<0.007$, $A\approx 0.01$, $0.015$, $0.026$ and $0.03$). Periodic regimes are observed when $\nu _{{ D}}/\nu _{{ f}}$ becomes rational: $\nu _{{ D}}/\nu _{{ f}}=m/n$, $m$, $n\in \mathbb {N}$, and only one dominating frequency $\nu _{{ D}}/m=\nu _{{ f}}/n$ together with its harmonics are present in the spectrum. Such regions are seen in figure 3 near $A=0.0075$ ($\nu _{{ D}}/\nu _{{ f}}=3/4$), $A=0.0125$ ($\nu _{{ D}}/\nu _{{ f}}=5/7$), $A=0.02$ ($\nu _{{ D}}/\nu _{{ f}}=2/3$), $A=0.028$ ($\nu _{{ D}}/\nu _{{ f}}=3/5$) and $A=0.035$ ($\nu _{{ D}}/\nu _{{ f}}=1/2$) as windows in $A$ with few nearly vertical spectral lines. Rational values of the frequency ratio $\nu _{{ D}}/\nu _{{ f}}$ form the Farey tree structure (Allen Reference Allen1983; Gonzalez & Piro Reference Gonzalez and Piro1985) since there exists a locking region with ratio $(m+m^{\prime })/(n+n^{\prime })$ between two regions with ratios $m/n$ and $m^{\prime }/n^{\prime }$, where $m^{\prime}$, $n^{\prime} \in \mathbb{N}$ and $m^{\prime} \neq m$, $n^{\prime} \neq n$. Moreover, the larger the denominator $n$ of the frequency ratio $\nu _{{ D}}/\nu _{{ f}}$, the thinner the locking interval on the spectrogram. It should also be noted that period-doubling bifurcations occur inside some of the regions (for $\nu _{{ D}}/\nu _{{ f}}=2/3$ and $1/2$) that are observed as the appearance of harmonics with half the dominant frequency. At values $A>0.042$, the spectra are seen to become noisy with a nearly continuous range of frequencies in the signal. This is due to the formation of secondary shock waves in the reaction zone behind the leading detonation shock in which case signal $D(t)$ contains discontinuities that are challenging to characterize by the Fourier spectra and are the cause of the noisy signal. The analysis of this region requires different numerical approaches and is outside the scope of the present paper.

In figure 4, we show the bifurcation diagram displaying normalized local maxima $D_{{ max}}/D_{{ CJ}}$ in $D(t)$ as a function of $A$. Periodic locked signals with $\nu _{{ D}}/\nu _{{ f}}=m/n$ have only a small number $m$ of maxima, while quasiperiodic solutions give a large set of distinct peaks. Period-doubling bifurcations can also be observed in the diagram.

Figure 4. The bifurcation diagram in the plane of the normalized local maxima $D_{{ max}}/D_{{ CJ}}$ of detonation velocity as a function of $A$ for the wavenumber $k=0.1$.

4. Arnold tongues and devil's staircase

Next, we investigate the nature of the locking regions seen in figures 3 and 4. Specifically, we perform two-parameter studies and analyse the period number of $D(t)$ for a range of amplitudes from $A=0$ to $0.02$ and range of wavenumbers from $k=0$ to $0.28$. The period number is calculated as the number of distinct peaks in the signal $D(t)$ for $t$ from $t=2500$ to $3000$. Figure 5 summarizes the results of the calculations. In the wedge-shaped areas in this figure, the period number is small. The areas have a structure similar to that of Arnold tongues for mode-locking transition to chaos in various dynamical systems. In such simple systems as a circle map (Jensen, Bak & Bohr Reference Jensen, Bak and Bohr1984; Schuster & Just Reference Schuster and Just2005), these tongues expand with increasing nonlinearity and eventually fill up the whole horizontal axis at some critical level of nonlinear coupling. For the coupling values greater than this critical one, the tongues overlap, and chaotic behaviour arises because the phase trajectory jumps between various overlapping resonances in an erratic way.

Figure 5. Arnold tongues in the heat map of period number as a function of amplitude $A$ and wavenumber $k$.

To estimate critical amplitude, $A_{{ c}}$, for our case, we investigate the behaviour of the rotation number and estimate the fractal dimension of the set complementary to the mode-locked states. Following the methods of analysis of forced oscillators (e.g. Bohr, Bak & Jensen Reference Bohr, Bak and Jensen1984), we take a large number $N$ of stroboscopic samples from the time series $D(t)$ taken at intervals equal to the period of the upstream oscillations, $T_{{ f}}=1/\nu _{{ f}}$. Then, we calculate the sample phases $\phi _{j}$, $j=1,\dots,N$, as complex arguments of points on the phase trajectory in the plane $(D-\bar {D},\dot {D})$, where $\dot {D}={\rm d}D/{\rm d}t$. The rotation number is defined as a limit of the ratio between the total phase increase and the number of samples: $W=\lim _{N\to \infty }(\phi _{N}-\phi _{1})/N$. This is in essence another form of the frequency ratio $\nu _{{ D}}/\nu _{{ f}}$ discussed above in relation to figure 3. The dependence of the rotation number on forcing frequency represents the Cantor function (often called the devil's staircase, see Bak Reference Bak1986) that has constant rational values within synchronization regions. Analytical and experimental studies show that this function arises in various nonlinear systems as part of the transition to chaos via quasiperiodicity or mode locking (Bohr et al. Reference Bohr, Bak and Jensen1984; Jensen et al. Reference Jensen, Bak and Bohr1984; Schuster & Just Reference Schuster and Just2005). Amazingly enough, certain universal scaling laws have been found near the critical parameters, the main one being the fractal dimension of the set complementary to the steps of the critical staircase. At transition to chaos, resonance regions fill up the whole frequency axis while quasiperiodic trajectories form the Cantor set of zero Lebesgue measure with non-integer dimension $d=0.87$ (Bak Reference Bak1986).

In figure 6, we show how the fractal dimension varies with the increasing amplitude $A$. The numerical results have been obtained with steps in $k$ equal to $0.00075$, which allows us to estimate accurately only the locations of the widest tongues in the analysed regions. They correspond to rotation numbers $W=1$, $2/3$, $1/2$, $2/5$ and $1/3$ (see also figure 7). Simply put, the values of dimension are measures of complexity for different devil's staircases. They indicate how the number of points scales in a neighbourhood of an arbitrary point of the Cantor set, if one changes the radius of this neighbourhood. The fractal dimension does not depend on the choice of the mode-locked regions because of the self-similar nature of the devil's staircase, and we use two approaches to estimate this value. The first one follows Jensen et al. (Reference Jensen, Bak and Bohr1983); Bak (Reference Bak1986) and uses a range of scales $r$ to find the total width $S(r)$ of all steps wider than a given scale $r$. The space between the steps is found as $N(r) = S - S(r)$, where $S$ is the length of the $k$ range under consideration. When the tongues expand, the synchronization regions grow, and the space between them $N(r)$ eventually shrinks to a Cantor set and, if measured on scale $r$, indicates a power law with specific exponent $N(r)/r \sim r^{-d}$. One can find the dimension $d$ from the relation $\log N(r) \sim (1-d) \log r$ via linear regression. This approach works well when one has a sufficiently large number of the detected tongues up to very narrow lengths, which is not our case. The second approach to estimating the dimension is to solve equation $\sum _{i}(s_{i}/\bar {s})^{d}=1$ for $d$, where $\bar {s}$ is the gap between two steps and $s_{i}$ are smaller gaps found inside the first larger gap (Cvitanovic et al. Reference Cvitanovic, Jensen, Kadanoff and Procaccia1985; Haucke & Ecke Reference Haucke and Ecke1987; Sagdeev, Usikov & Zaslavsky Reference Sagdeev, Usikov and Zaslavsky1988). One can choose any two locking regions corresponding to ratios $m/n$ and $m^{\prime }/n^{\prime }$ to calculate the large gap $\bar {s}$, but the Farey tree structure between rational $m/n$ and $m^{\prime }/n^{\prime }$ should be used to identify the smaller gaps $s_{i}$. This option gives good results for experimental data with a moderate number of detected steps, hence it is also a preferable method for our case. Based on this analysis, we have found an estimate for the critical amplitude as $A_{{ c}}=0.0121$. Here, we should emphasize that mathematically one observes an abrupt transition of the dimension value from $0.87$ to $1$ with the decreasing amplitude, but numerical estimates show a smooth cross-over as in figure 6 even for simpler dynamical systems such as the circle map (Jensen et al. Reference Jensen, Bak and Bohr1984, § V).

Figure 6. The fractal dimensions of the set complementary to the mode-locked states computed by methods from Jensen, Bak & Bohr (Reference Jensen, Bak and Bohr1983) (red dots) and Cvitanovic et al. (Reference Cvitanovic, Jensen, Kadanoff and Procaccia1985) (blue dots). The black dashed line shows the universal value $d=0.87$.

Figure 7. The critical devil's staircase for the rotation number $W$ as a function of wavenumber $k$ at $A_{{ c}}=0.0121$.

In figure 7, we construct the dependence of $W$ on $k$ at $A=A_{{ c}}$. It clearly shows the presence of plateaus for $k$ in the locking regions. We note certain local universal features of the mode locking (Jensen et al. Reference Jensen, Bak and Bohr1984) such as the slope discontinuity near rational rotation numbers in figure 7 and several period-doubling bifurcations when amplitude $A$ increases for fixed $k$ in figure 4. Just above the critical amplitude $A_{{ c}}$, where the synchronization regions begin to intersect, chaotic solutions may arise. However, they are difficult to find in practice since their measure is only slightly above zero near $A_{{ c}}$ (Schuster & Just Reference Schuster and Just2005). Precisely how the chaotic solutions arise at $A>A_{{ c}}$ is an interesting problem that we believe should be investigated in the future. The nature of the Arnold tongues and the path to chaos for the detonation wave are likely to be more complex than in the simple case of the one-dimensional circle map. For example, the critical line may not necessarily be horizontal and various scaling relations may change compared with those of the circle map if the observed dynamics belong to a different universality class (Jensen et al. Reference Jensen, Bak and Bohr1983; Bohr et al. Reference Bohr, Bak and Jensen1984). Our present results should be looked at as first estimates in this context.