2. The main governing equations of reactive flow

We assume that the medium is described by the following system of equations expressing respectively the laws of conservation of mass, momentum and energy, and the chemical heat release:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}{\it\rho}}{\text{D}t}+{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=\text{div}(\unicode[STIX]{x1D64F}), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}e}{\text{D}t}=\unicode[STIX]{x1D64F}\boldsymbol{ : }\unicode[STIX]{x1D63F}-\text{div}(\boldsymbol{q}_{e}), & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}{\it\Lambda}}{\text{D}t}={\it\rho}\tilde{W}-\text{div}(\boldsymbol{q}_{{\it\Lambda}}). & \displaystyle\end{eqnarray}$$

$\text{D}/\text{D}t=\partial /\partial t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$

$\mathit{v}=1/{\it\rho}$

${\it\rho}$

$\boldsymbol{u}$

$\unicode[STIX]{x1D64F}$

$\unicode[STIX]{x1D63F}=(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})/2$

$\unicode[STIX]{x1D64F}\boldsymbol{ : }\unicode[STIX]{x1D63F}=\sum _{i,j}\unicode[STIX]{x1D61B}_{ij}\unicode[STIX]{x1D60B}_{ij}$

$e=e_{i}-\tilde{Q}{\it\Lambda}$

$\tilde{Q}$

$\boldsymbol{q}_{e}$

$\boldsymbol{q}_{{\it\Lambda}}$

$\tilde{W}({\it\Lambda},T)$

${\it\Lambda}$

${\it\Lambda}=0$

${\it\Lambda}=1$

We make the following standard modelling assumptions (e.g. Williams Reference Williams1985).

1. (1) The fluid is Newtonian, with the Stokes assumption on the bulk viscosity, so that $\unicode[STIX]{x1D64F}=-(p+(2/3){\it\mu}\text{div}(\boldsymbol{u}))\unicode[STIX]{x1D644}+2{\it\mu}\unicode[STIX]{x1D63F}$ , where ${\it\mu}$ is the dynamic viscosity, $p$ is the pressure and $\unicode[STIX]{x1D644}$ is the identity tensor.

2. (2) The species flux is given by Fick’s law, $\boldsymbol{q}_{{\it\Lambda}}=-{\it\rho}d\boldsymbol{{\rm\nabla}}{\it\Lambda}$ , with $d$ denoting the diffusion coefficient.

3. (3) The energy flux has contributions from both the heat conduction (given by Fourier’s law) and the species diffusion, so that $\boldsymbol{q}_{e}=-{\it\kappa}\boldsymbol{{\rm\nabla}}T+\tilde{Q}{\it\rho}d\boldsymbol{{\rm\nabla}}{\it\Lambda}$ , where ${\it\kappa}$ is the heat diffusion coefficient.

4. (4) The medium is a perfect gas described by the ideal-gas equation of state, $p={\it\rho}RT$ , with the internal energy given by $e_{i}=p\mathit{v}/({\it\gamma}-1)$ , where $R$ is the universal gas constant divided by the molecular weight and ${\it\gamma}$ is the ratio of specific heats, assumed to be constant.

5. (5) For simplicity, we take the rate of reaction to be $\tilde{W}=\tilde{k}(1-{\it\Lambda})\exp (-\tilde{E}/RT)$ , with the added ignition-temperature assumption that $\tilde{W}=0$ for $T<T_{i}$ for some temperature, $T_{i}$ . With this assumption, used throughout the rest of the paper, reactions are switched off ahead of the propagating detonation wave. Here, $\tilde{k}$ is the rate constant and $\tilde{E}$ is the activation energy. More general rate functions can be considered, in principle, as long as appropriate sensitivity to temperature is preserved.

We can then rewrite (2.1)–(2.4) as

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}{\it\rho}}{\text{D}t}+{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(-\left(p+\frac{2}{3}{\it\mu}\text{div}(\boldsymbol{u})\right)\unicode[STIX]{x1D644}+2{\it\mu}\unicode[STIX]{x1D63F}\right), & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \hspace{-4.79993pt}{\it\rho}\frac{\text{D}e}{\text{D}t}=-p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}-\frac{2}{3}{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})^{2}+{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u})+{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}\boldsymbol{{\rm\nabla}}T-\tilde{Q}{\it\rho}d\boldsymbol{{\rm\nabla}}{\it\Lambda}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}{\it\Lambda}}{\text{D}t}={\it\rho}\tilde{W}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\text{}{\it\rho}\text{}d\boldsymbol{{\rm\nabla}}{\it\Lambda}). & \displaystyle\end{eqnarray}$$

$e=e_{i}-\tilde{Q}{\it\lambda}=RT/({\it\gamma}-1)-\tilde{Q}{\it\Lambda}$

(2.9) $$\begin{eqnarray}\displaystyle & & \displaystyle {\it\rho}\frac{\text{D}T}{\text{D}t}-\frac{{\it\gamma}-1}{R{\it\gamma}}\frac{\text{D}p}{\text{D}t}=\frac{{\it\gamma}-1}{R{\it\gamma}}\left(\tilde{Q}{\it\rho}\tilde{W}-\frac{2}{3}{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})^{2}+{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u})\right.\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left.(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}\boldsymbol{{\rm\nabla}}T)\right).\qquad\end{eqnarray}$$

We focus now on the two-dimensional case for simplicity. Consider a localized wave moving into an equilibrium quiescent state and let ${\it\rho}_{a}$ , $p_{a}$ , $T_{a}$ and $u_{a}=\sqrt{p_{a}/{\it\rho}_{a}}$ denote the density, pressure, temperature and Newtonian sound speed in the fresh mixture ahead of the wave respectively. We rescale the dependent variables with respect to this reference state. The independent variables are rescaled as follows:

(2.10a−c ) $$\begin{eqnarray}x=\frac{X-D_{0}t}{x_{0}},\quad y=\frac{Y}{y_{0}},\quad {\it\tau}=\frac{t}{t_{0}},\end{eqnarray}$$

where $X$ and $Y$ are the original spatial variables and $D_{0}$ is a typical wave speed, which is to be determined in the process of deriving the asymptotic model by requiring non-triviality of the leading-order corrections. The length scales, $x_{0}$ , $y_{0}$ , and the time scale, $t_{0}$ , are chosen to reflect the appropriate physics of weakly nonlinear waves. We assume that ${\it\epsilon}=x_{0}/(u_{a}t_{0})$ is small, which means that the spatial scale of interest in the $x$ direction, which is related to the size of the reaction zone (the half-reaction length defined below), is small compared with the typical distance covered by acoustic waves in time $t_{0}$ . For the transverse dimension, we assume the scaling $y_{0}=x_{0}/\sqrt{{\it\epsilon}}$ . This follows from the fact that, along a weakly curved front, a distance ${\it\epsilon}$ in the normal direction corresponds to a distance $O(\sqrt{{\it\epsilon}})$ in the transverse direction.

Several dimensionless groups appear upon rescaling of the governing equations. We define the Reynolds, Prandtl and Lewis numbers respectively as follows:

(2.11a−c ) $$\begin{eqnarray}\mathit{Re}=\frac{{\it\rho}_{a}u_{a}x_{0}}{{\it\mu}},\quad \mathit{Pr}=\frac{c_{p}{\it\mu}}{{\it\kappa}},\quad \mathit{Le}=\frac{{\it\kappa}}{{\it\rho}_{a}c_{p}d},\end{eqnarray}$$

where $c_{p}={\it\gamma}R/({\it\gamma}-1)$ . We also introduce the following dimensionless parameters:

(2.12a−c ) $$\begin{eqnarray}Q=\frac{\tilde{Q}}{RT_{a}},\quad E=\frac{\tilde{E}}{RT_{a}},\quad K=t_{0}\tilde{k}\exp (-E).\end{eqnarray}$$

From now on, all variables are dimensionless even though we retain the same notation as for the dimensional quantities above.

Writing $\boldsymbol{u}=(U,V)^{\text{T}}$ , it is convenient to introduce the differential operator

(2.13) $$\begin{eqnarray}L=\partial _{{\it\tau}}+\frac{1}{{\it\epsilon}}(U-D_{0})\partial _{x}+\frac{1}{\sqrt{{\it\epsilon}}}V\partial _{y}.\end{eqnarray}$$

Then, the dimensionless compressible reactive Navier–Stokes equations become

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle L[{\it\rho}]+{\it\rho}\left(\frac{1}{{\it\epsilon}}U_{x}+\frac{1}{\sqrt{{\it\epsilon}}}V_{y}\right)=0, & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}L[U]+\frac{1}{{\it\epsilon}}p_{x}=\frac{1}{3{\it\epsilon}\mathit{Re}}(U_{xx}+\sqrt{{\it\epsilon}}V_{xy})+\frac{1}{{\it\epsilon}\mathit{Re}}(U_{xx}+{\it\epsilon}U_{yy}), & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}L[V]+\frac{1}{\sqrt{{\it\epsilon}}}p_{y}=\frac{1}{3{\it\epsilon}\mathit{Re}}(\sqrt{{\it\epsilon}}U_{xy}+{\it\epsilon}V_{yy})+\frac{1}{{\it\epsilon}\mathit{Re}}(V_{xx}+{\it\epsilon}V_{yy}), & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle {\it\rho}L[T]-\frac{({\it\gamma}-1)}{{\it\gamma}}L[p] & = & \displaystyle \frac{{\it\gamma}-1}{{\it\gamma}}\left(Q{\it\rho}W-\frac{2}{3{\it\epsilon}\mathit{Re}}(U_{x}+\sqrt{{\it\epsilon}}V_{y})^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{1}{{\it\epsilon}\mathit{Re}}(U_{x}^{2}+{\it\epsilon}U_{y}^{2}+V_{x}^{2}+{\it\epsilon}V_{y}^{2})\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{\it\gamma}-1}{{\it\gamma}}\frac{1}{{\it\epsilon}\mathit{Re}}(U_{x}^{2}+\sqrt{{\it\epsilon}}U_{y}V_{x}+\sqrt{{\it\epsilon}}V_{x}U_{y}+{\it\epsilon}V_{y}^{2})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{{\it\epsilon}\mathit{Re}\,\mathit{Pr}}(T_{xx}+{\it\epsilon}T_{yy}),\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle {\it\rho}L[{\it\Lambda}]={\it\rho}W+\frac{1}{{\it\epsilon}\mathit{Re}\,\mathit{Pr}\,\mathit{Le}}(({\it\rho}{\it\Lambda}_{x})_{x}+{\it\epsilon}({\it\rho}{\it\Lambda}_{y})_{y}), & & \displaystyle\end{eqnarray}$$

where $W$ is defined as

(2.19) $$\begin{eqnarray}W=K(1-{\it\Lambda})\exp \left[E\left(1-\frac{1}{T}\right)\right].\end{eqnarray}$$

3. Weakly nonlinear approximation

We develop an asymptotic simplification of the above general formulation by considering a weakly nonlinear detonation wave, for which we assume that the heat release is small, the activation energy is large and the Newtonian limit, $0<{\it\gamma}-1\ll 1$ , holds. We also identify the conditions on $\mathit{Le}$ , $\mathit{Pr}$ and $\mathit{Re}$ that allow us to neglect the dissipative effects in the main reaction region. To be precise, we start from (2.14)–(2.18) and make the following assumptions.

1. (1) $K=k/{\it\epsilon}$ , $k=O(1)$ . This assumption is chosen to ensure that the reaction rate affects the leading-order expansion of ${\it\Lambda}$ . Since $K=t_{0}\tilde{k}\exp (-E)$ , this assumption implies that the characteristic time scale, $t_{0}$ , of weakly nonlinear detonations is large compared with the collision time, $\tilde{t}_{c}=1/\tilde{k}$ , namely, $t_{0}=O[\tilde{t}_{c}/({\it\epsilon}\exp (-E))]$ .

2. (2) $({\it\gamma}-1)Q/{\it\gamma}={\it\epsilon}^{2}q$ , $q=O(1)$ . This assumption implies that the heat release does not play a role at the linear level. It does not mean that the chemistry is unimportant, but that the heat release must have the appropriate size to balance the nonlinear effects. The extra factor of $({\it\gamma}-1)/{\it\gamma}$ in front of $Q$ arises naturally in the governing equations (see (2.17)) and is retained in the definition of $q$ . With the further assumption below of small ${\it\gamma}-1$ , this implies that $Q$ is $O({\it\epsilon})$ .

3. (3) $E={\it\theta}/{\it\epsilon}^{2},$ ${\it\theta}=O(1)$ . This ensures that small temperature deviations – which are $O({\it\epsilon}^{2})$ for weak shocks in the Newtonian limit – have an $O(1)$ influence on the reaction rate.

4. (4) ${\it\gamma}-1={\it\gamma}_{1}{\it\epsilon}$ , ${\it\gamma}_{1}=O(1)$ . This assumption is needed to balance the temperature fluctuations with both the acoustics and the chemistry at the same order. Without this assumption, the leading-order corrections for density, velocity and temperature all behave in the same way, as in a weakly nonlinear inert shock (Rosales Reference Rosales1991; Hunter Reference Hunter1995). As we show later, having a temperature profile that is different from the density/velocity profiles is crucial, as it allows the model derived here to incorporate the dynamical instabilities of detonation waves.

5. (5) $\mathit{Le}^{-1}$ and $\mathit{Pr}^{-1}$ are $O(1)$ , while $\mathit{Re}^{-1}$ is $o({\it\epsilon})$ . These assumptions follow from the kinetic theory of gases which yields $\mathit{Re}\sim {\it\epsilon}^{-2}\exp ({\it\theta}/{\it\epsilon}^{2})$ . Therefore, all of the dissipative terms in the governing equations have coefficients that are exponentially small as ${\it\epsilon}\rightarrow 0$ . This justifies their neglect to some extent. It is clear, however, that the role of the dissipative terms is non-trivial because of their singular perturbation nature. We discuss this point further in the conclusions.

To elucidate the intuition behind the asymptotic ordering above, we recall the following well-known fact for weak shocks (Whitham Reference Whitham1974). If the shock strength is measured by the relative jump in pressure across the shock, ${\rm\Delta}p=(p_{s}-p_{a})/p_{a}$ (with subscripts $s$ and $a$ denoting post- and pre-shock states respectively), then the shock Mach number is $\mathit{Ma}=1+[({\it\gamma}+1)/(4{\it\gamma})]\,{\rm\Delta}p+O(({\rm\Delta}p)^{2})$ and the shock temperature is $T_{s}/T_{a}=1+[({\it\gamma}-1)/{\it\gamma}]\,{\rm\Delta}p+O(({\rm\Delta}p)^{2})$ . Therefore, for weak shocks, with $\mathit{Ma}-1=O({\it\epsilon})$ , in the Newtonian limit, ${\it\gamma}-1=O({\it\epsilon})$ , the leading-order temperature correction is $O({\it\epsilon}^{2})$ . That is, all variables have an $O({\it\epsilon})$ jump across the shock, but the temperature jump is smaller, only $O({\it\epsilon}^{2})$ . In the chosen asymptotic approximation, we therefore expect similar temperature behaviour in the reaction zone as well, at least for non-dissipative detonations. This also means that the leading-order temperature perturbation is decoupled from the other variables, which is crucial in terms of instability, as explained later in the paper.

Now, we assume the following expansions in the reaction zone:

(3.1) $$\begin{eqnarray}\displaystyle & {\it\rho}=1+{\it\rho}_{1}{\it\epsilon}+{\it\rho}_{2}{\it\epsilon}^{3/2}+{\it\rho}_{3}{\it\epsilon}^{2}+o({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & T=1+T_{1}{\it\epsilon}+T_{2}{\it\epsilon}^{3/2}+T_{3}{\it\epsilon}^{2}+o({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & p=1+p_{1}{\it\epsilon}+p_{2}{\it\epsilon}^{3/2}+p_{3}{\it\epsilon}^{2}+o({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}=\boldsymbol{u}_{1}{\it\epsilon}+\boldsymbol{u}_{2}{\it\epsilon}^{3/2}+\boldsymbol{u}_{3}{\it\epsilon}^{2}+o({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & {\it\Lambda}={\it\lambda}+o({\it\epsilon}). & \displaystyle\end{eqnarray}$$

$O(\sqrt{{\it\epsilon}})$

Expanding $p={\it\rho}T$ , we find that $p_{1}={\it\rho}_{1}+T_{1}$ , $p_{2}={\it\rho}_{2}+T_{2}$ and $p_{3}={\it\rho}_{3}+T_{3}+{\it\rho}_{1}T_{1}$ . Using these relations to eliminate pressure perturbations, (2.14)–(2.18) yield

(3.6) $$\begin{eqnarray}\displaystyle & & \displaystyle (-D_{0}{\it\rho}_{1x}+U_{1x})+\sqrt{{\it\epsilon}}(-D_{0}{\it\rho}_{2x}+U_{2x}+V_{1y})\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\epsilon}({\it\rho}_{1{\it\tau}}-D_{0}{\it\rho}_{3x}+U_{1}{\it\rho}_{1x}+U_{3x}+{\it\rho}_{1}U_{1x}+V_{2y})=o({\it\epsilon}),\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & & \displaystyle (-D_{0}U_{1x}+T_{1x}+{\it\rho}_{1x})+\sqrt{{\it\epsilon}}(-D_{0}U_{2x}+T_{2x}+{\it\rho}_{2x})\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\epsilon}(U_{1{\it\tau}}-D_{0}U_{3x}+U_{1}U_{1x}+T_{1}{\it\rho}_{1x}+T_{3x}+{\it\rho}_{3x}-D_{0}{\it\rho}_{1}U_{1x})=o({\it\epsilon}),\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & & \displaystyle (-D_{0}V_{1x})+\sqrt{{\it\epsilon}}(-D_{0}V_{2x}+T_{1y}+{\it\rho}_{1y})\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\epsilon}(V_{1{\it\tau}}-D_{0}V_{3x}+U_{1}V_{1x}+T_{2y}+{\it\rho}_{2y})=o({\it\epsilon}),\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & & \displaystyle (-D_{0}T_{1x})+\sqrt{{\it\epsilon}}(-D_{0}T_{2x})\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\epsilon}\left(T_{1{\it\tau}}-D_{0}T_{3x}+U_{1}T_{1x}+\frac{{\it\gamma}_{1}}{{\it\gamma}}D_{0}({\it\rho}_{1x}+T_{1x})\right)={\it\epsilon}q{\it\omega}+o({\it\epsilon}),\end{eqnarray}$$

(3.10) $$\begin{eqnarray}-D_{0}{\it\lambda}_{x}={\it\omega}+o(1).\end{eqnarray}$$

Because of the weak heat-release assumption, we need to expand only the reaction-progress variable and the reaction rate to the leading order. As we shall see later, the leading-order corrections to temperature are of order ${\it\epsilon}^{2}$ , i.e. $T_{1}=T_{2}=0$ , and therefore the leading-order reaction rate is given by

(3.11) $$\begin{eqnarray}{\it\omega}=k(1-{\it\lambda})\exp ({\it\theta}T_{3}).\end{eqnarray}$$

Balancing $O(1)$ terms, we find that

(3.12) $$\begin{eqnarray}\left[\!\begin{array}{@{}cccc@{}}-D_{0} & 1 & 0 & 0\\ 1 & -D_{0} & 0 & 1\\ 0 & 0 & -D_{0} & 0\\ 0 & 0 & 0 & -D_{0}\end{array}\!\right]\left[\!\begin{array}{@{}c@{}}{\it\rho}\\ U\\ V\\ T\end{array}\!\right]_{1x}=\left[\!\begin{array}{@{}c@{}}0\\ 0\\ 0\\ 0\!\end{array}\right].\end{eqnarray}$$

This homogeneous system has non-trivial solutions if and only if the coefficient matrix (denoted by $\unicode[STIX]{x1D63C}$ from now on) is singular. Therefore, we must have $D_{0}=\pm 1$ . We focus on a right-going wave for which $D_{0}=1$ .

After substituting $D_{0}=1,$ we integrate the first equation in (3.12) to obtain $U_{1}=u(x,y,{\it\tau})+\bar{U}(y,{\it\tau}),{\it\rho}_{1}=u(x,y,{\it\tau})+\bar{{\it\rho}}(y,{\it\tau})$ for some, so far arbitrary, functions $\bar{U}(y,{\it\tau})$ and $\bar{{\it\rho}}(y,{\it\tau})$ . In order for our asymptotic expansions to hold as we approach the upstream state, which is assumed to be quiescent and uniform at all times, we require ${\it\rho}_{1}=U_{1}\equiv 0$ as $x\rightarrow \infty$ (upstream of the wave). As a consequence, it follows that $\bar{U}(y,{\it\tau})=\bar{{\it\rho}}(y,{\it\tau})\equiv 0$ . Had we considered the case where the wave moves into a non-uniform background, the functions $\bar{{\it\rho}}$ and $\bar{U}$ would provide the freedom needed to match the asymptotic expansion in the wave zone to the flow ahead (e.g. Rosales & Majda Reference Rosales and Majda1983). A similar reasoning can be used to deduce from the third and fourth equations in (3.12) that $T_{1}=V_{1}\equiv 0$ .

At $O({\it\epsilon}^{1/2})$ , we obtain

(3.13) $$\begin{eqnarray}\left[\begin{array}{@{}cccc@{}}-1 & 1 & 0 & 0\\ 1 & -1 & 0 & 1\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -1\end{array}\right]\left[\begin{array}{@{}c@{}}{\it\rho}\\ U\\ V\\ T\end{array}\right]_{2x}=\left[\begin{array}{@{}c@{}}0\\ 0\\ -u_{y}\\ 0\end{array}\right],\end{eqnarray}$$

with the same matrix of coefficients, $\unicode[STIX]{x1D63C}$ , as before. Since $\unicode[STIX]{x1D63C}$ is singular, this system has solutions if and only if the right-hand side is orthogonal to $\boldsymbol{l}_{0}=[1\,1\,0\,1]$ , which spans the left null space of $\unicode[STIX]{x1D63C}$ . Clearly, this condition is always satisfied here. It should be noted that from the third and fourth equations in (3.13), we obtain that $V_{2x}=u_{y}$ and $T_{2}=0$ , i.e. letting

(3.14) $$\begin{eqnarray}v=V_{2},\end{eqnarray}$$

we find

(3.15) $$\begin{eqnarray}v_{x}=u_{y}.\end{eqnarray}$$

Finally, at $O({\it\epsilon})$ , we obtain

(3.16) $$\begin{eqnarray}\left[\begin{array}{@{}cccc@{}}-1 & 1 & 0 & 0\\ 1 & -1 & 0 & 1\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -1\end{array}\right]\left[\begin{array}{@{}c@{}}{\it\rho}\\ U\\ V\\ T\end{array}\right]_{3x}=\left[\begin{array}{@{}c@{}}-u_{{\it\tau}}-2uu_{x}-v_{y}\\ -u_{{\it\tau}}\\ -{\it\rho}_{2y}\\ q{\it\omega}-{\displaystyle \frac{{\it\gamma}_{1}}{{\it\gamma}}}u_{x}\end{array}\right].\end{eqnarray}$$

It should be noticed that all terms in (3.16) are $O(1)$ . The solvability condition for this system is that the right-hand side be orthogonal to $\boldsymbol{l}_{0}=[1\,1\,0\,1]$ , which yields

(3.17) $$\begin{eqnarray}\displaystyle 2u_{{\it\tau}}+2uu_{x}+v_{y}=q{\it\omega}-\frac{{\it\gamma}_{1}}{{\it\gamma}}u_{x}. & & \displaystyle\end{eqnarray}$$

This equation together with

(3.18) $$\begin{eqnarray}\displaystyle & v_{x}=u_{y}, & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{x}=-{\it\omega}, & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}T_{3}=\frac{{\it\gamma}_{1}}{{\it\gamma}}u+q{\it\lambda},\end{eqnarray}$$

giving

(3.21) $$\begin{eqnarray}{\it\omega}=k(1-{\it\lambda})\exp \left({\it\theta}\left(\frac{{\it\gamma}_{1}}{{\it\gamma}}u+q{\it\lambda}\right)\right).\end{eqnarray}$$

Therefore, we obtain

(3.22) $$\begin{eqnarray}\displaystyle & 2u_{{\it\tau}}+2uu_{x}+{\displaystyle \frac{{\it\gamma}_{1}}{{\it\gamma}}}u_{x}+v_{y}=-q{\it\lambda}_{x}, & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & v_{x}=u_{y}, & \displaystyle\end{eqnarray}$$

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\lambda}_{x}=-k(1-{\it\lambda})\exp \left({\it\theta}\left(\frac{{\it\gamma}_{1}}{{\it\gamma}}u+q{\it\lambda}\right)\right). & \displaystyle\end{eqnarray}$$

It is convenient to further rescale the variables as

(3.25a−e ) $$\begin{eqnarray}x\mapsto x-{\it\gamma}_{1}/(2{\it\gamma}){\it\tau},\quad y\mapsto 2^{-1/2}q^{-1/4}y,\quad {\it\tau}\mapsto q^{-1/2}{\it\tau},\quad u\mapsto q^{1/2}u,\quad v\mapsto 2^{1/2}q^{3/4}v,\end{eqnarray}$$

where the scale for $u$ is chosen so that the travelling wave solution found later in § 4.1 has speed $1$ . We also choose ${\it\epsilon}=({\it\gamma}-1)/{\it\gamma}$ , which means that ${\it\gamma}_{1}={\it\gamma}$ . Other choices of ${\it\epsilon}$ are possible as long as ${\it\epsilon}$ and ${\it\gamma}-1$ are of the same asymptotic order when ${\it\epsilon}\rightarrow 0$ . These choices, although equivalent in the limit ${\it\epsilon}\rightarrow 0$ , will have some effect on the quantitative numerical predictions for finite ${\it\epsilon}$ . With this rescaling and the choice of ${\it\epsilon}$ , we obtain the asymptotic system of equations for weakly nonlinear detonations in its final form:

(3.26) $$\begin{eqnarray}\displaystyle & u_{{\it\tau}}+uu_{x}+v_{y}=-{\textstyle \frac{1}{2}}{\it\lambda}_{x}, & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & v_{x}=u_{y}, & \displaystyle\end{eqnarray}$$

(3.28) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{x}=-k(1-{\it\lambda})\exp [{\it\theta}(\sqrt{q}u+q{\it\lambda})]. & \displaystyle\end{eqnarray}$$

Equations (3.26)–(3.28) must be supplemented by the appropriate jump conditions across shocks. If the shock locus is defined by ${\it\phi}(x,y,{\it\tau})=x-s(y,{\it\tau})=0$ and $[z]$ denotes the jump of $z$ across the shock, the Rankine–Hugoniot conditions are

(3.29) $$\begin{eqnarray}\displaystyle & s_{{\it\tau}}[u]-{\textstyle \frac{1}{2}}[u^{2}]+s_{y}[v]=0, & \displaystyle\end{eqnarray}$$

(3.30) $$\begin{eqnarray}\displaystyle & s_{y}[u]+[v]=0, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\displaystyle & [{\it\lambda}]=0. & \displaystyle\end{eqnarray}$$

These equations follow from the conservation form of (3.26)–(3.28).

We recall that, implicit in the definition, (2.19), of the rate of reaction, $W$ , lies the ignition-temperature assumption, such that $W\equiv 0$ for $T<T_{i}/T_{a}$ . Thus, the leading-order reaction rate, ${\it\omega}$ , given by (3.11) also satisfies ${\it\omega}\equiv 0$ for $T_{3}<(T_{i}/T_{a}-1)/{\it\epsilon}^{2}$ or, equivalently, for $u+q{\it\lambda}<(T_{i}/T_{a}-1)/{\it\epsilon}^{2}$ . In order to prevent reactions from occurring in the ambient state, a reasonable constraint on the ignition temperature is that it be larger than the ambient temperature, $T_{a}$ . Furthermore, for the ZND solution to exist, it is necessary that the ignition temperature, $T_{i}$ , be smaller than the temperature at the von Neumann state of the ZND solution. In the weakly nonlinear regime considered in this paper, we must therefore have $0<T_{i}-T_{a}=O({\it\epsilon}^{2})$ . This should be viewed as a modelling assumption about the chemical kinetics. It should be noticed that the ignition temperature, if taken to be $T_{i}\gtrapprox T_{a}$ , has no effect on the ZND profiles and their stability properties. Simply, it states that the reactions occur only after the shock, and therefore ${\it\omega}=0$ ahead of the wave. When considering viscous effects, however, the travelling wave profiles and their dynamical evolution will depend on the ignition temperature, and $T_{i}$ should therefore be considered as another parameter that affects the properties of the solutions.

For the convenience of the reader, in table 1, we collect the relations between various dimensionless and dimensional quantities.

Table 1.

Summary of scaling relationships.

The relationships between the asymptotic variables $u$ , $v$ and ${\it\lambda}$ and the physical (dimensionless, rescaled with the upstream state, except for ${\it\Lambda}$ ) quantities ${\it\rho},U,V,T$ and ${\it\Lambda}$ are given by

(3.32) $$\begin{eqnarray}\displaystyle & {\it\rho}=1+{\it\epsilon}\sqrt{q}u+O({\it\epsilon}^{3/2}), & \displaystyle\end{eqnarray}$$

(3.33) $$\begin{eqnarray}\displaystyle & U={\it\epsilon}\sqrt{q}u+O({\it\epsilon}^{3/2}), & \displaystyle\end{eqnarray}$$

(3.34) $$\begin{eqnarray}\displaystyle & V={\it\epsilon}^{3/2}\sqrt{2}q^{3/4}v+O({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$

(3.35) $$\begin{eqnarray}\displaystyle & T=1+{\it\epsilon}^{2}(\sqrt{q}u+q{\it\lambda})+O({\it\epsilon}^{5/2}), & \displaystyle\end{eqnarray}$$

(3.36) $$\begin{eqnarray}\displaystyle & {\it\Lambda}={\it\lambda}+O({\it\epsilon}). & \displaystyle\end{eqnarray}$$

It is worth noting here that, after the last transformation given in (3.25), the spatial variable $x$ corresponds to a frame of reference moving with speed $(1+{\it\epsilon}/2)u_{a}$ (see table 1). Thus, the speed of the detonation in the asymptotic model is to be interpreted relative to this moving frame. This means that the detonation speed, as predicted by the asymptotic model, is given by

(3.37) $$\begin{eqnarray}\displaystyle D=\left(1+\frac{{\it\epsilon}}{2}\right)u_{a}+{\it\epsilon}\sqrt{q}\frac{u_{s}}{2}+O({\it\epsilon}^{2})=\frac{1}{\sqrt{{\it\gamma}}}\left(1+\frac{{\it\epsilon}}{2}+\frac{{\it\epsilon}}{2}\sqrt{q}u_{s}\right)c_{a}+O({\it\epsilon}^{2}), & & \displaystyle\end{eqnarray}$$

where $u_{s}(t)$ denotes the post-shock state (and therefore $u_{s}/2$ is the shock speed in the model) and $c_{a}$ is the sound speed in the fresh mixture. The $\sqrt{q}$ factors appear because of the additional velocity rescaling in (3.25). Upon substitution of ${\it\gamma}=(1-{\it\epsilon})^{-1}$ and expansion in powers of ${\it\epsilon}$ in (3.37), the detonation Mach number is given by

(3.38) $$\begin{eqnarray}\mathit{Ma}=1+{\it\epsilon}\frac{\sqrt{q}u_{s}}{2}+O({\it\epsilon}^{2}).\end{eqnarray}$$

Therefore, as long as $u_{s}$ remains positive (which it does), the wave moves supersonically, i.e. $\mathit{Ma}>1$ .

We note here that without the chemical reaction, (3.26)–(3.28) reduce to

(3.39) $$\begin{eqnarray}\displaystyle & u_{{\it\tau}}+\left({\displaystyle \frac{u^{2}}{2}}\right)_{x}+v_{y}=0, & \displaystyle\end{eqnarray}$$

(3.40) $$\begin{eqnarray}\displaystyle & v_{x}=u_{y}, & \displaystyle\end{eqnarray}$$

${\it\phi}$

$\boldsymbol{{\rm\nabla}}{\it\phi}=(u,v)$

(3.41) $$\begin{eqnarray}{\it\phi}_{x{\it\tau}}+\left(\frac{{\it\phi}_{x}}{2}\right)_{x}^{2}+{\it\phi}_{yy}=0,\end{eqnarray}$$

which is the well-known unsteady transonic small-disturbance equation (Lin et al. Reference Lin, Reissner and Tsien1948). In nonlinear acoustics, this equation is known as the Zabolotskaya–Khokhlov equation (Zabolotskaya & Khokhlov Reference Zabolotskaya and Khokhlov1969). In one dimension, when $v\equiv 0$ , (3.26), (3.28) are the equations in Clavin & Williams (Reference Clavin and Williams2002). With a different reaction rate ( $q{\it\lambda}$ missing in the exponential), the model (3.26)–(3.28) appears in Rosales (Reference Rosales1989). When linearized, equations (3.26)–(3.28) reduce to those obtained in Clavin & Williams (Reference Clavin and Williams2009).

It may also be of interest to note that the nonlinearities in (3.26)–(3.28) are present only in the reaction rate in (3.28) and the advection term in (3.26), which are common with the corresponding one-dimensional model (Clavin & Williams Reference Clavin and Williams2002). As elaborated in the introduction, the advective nonlinearity arises quite generally from the nature of the hyperbolic system. The reaction-rate nonlinearity is of course particular to detonations and is explained by the choice of the distinguished limit under consideration. It should be kept in mind, however, that the evolution in both directions is nonlinearly coupled through not only (3.26)–(3.28), but also the shock conditions (3.29)–(3.31), in which the terms $s_{y}[u]$ and $s_{y}[v]$ represent direct nonlinear coupling of the longitudinal and transverse variations of the solution. A question arises as to the implications of these facts for the origin of cellular structures in gaseous detonations. Since triple points are very natural in the context of non-planar shock dynamics (even without chemical reactions, see, e.g., Sturtevant & Kulkarny Reference Sturtevant and Kulkarny1976), it is tempting to view cellular detonations as the result of a ‘combination’ of the one-dimensional longitudinal oscillatory dynamics (which leads to pulsating detonations) and the nonlinear multidimensional dynamics of the lead shock. However, the true mechanism behind the formation and evolution of detonation cells is likely to require subtler reasoning than this, precisely because of the nonlinear couplings just described – for example, recall that cellular detonations can arise in situations where purely one-dimensional perturbations are stable. In our opinion, much further analysis is still necessary to elucidate the physical mechanisms of cellular detonations.

In the following sections, we investigate the steady-state solutions of (3.26)–(3.28), their spectral stability and nonlinear dynamics, and demonstrate that the asymptotic theory contains the essential features of not only the steady-state one-dimensional solution, but also the unsteady and multidimensional travelling waves of the reactive Euler equations.

4. Travelling wave solutions and their linear stability

In this section, we analyse the travelling wave solutions of the asymptotic equations and their spectral stability, showing that the asymptotic solutions are analogous to their ZND counterparts. We start by presenting the one-dimensional travelling wave solutions in § 4.1, which are the basis for the one- and multidimensional stability analyses presented in § 4.2.

4.1. Travelling wave solutions of the asymptotic model

Seeking the one-dimensional travelling wave solutions of the inviscid asymptotic model (3.26)–(3.28) of the form $\bar{u}=\bar{u}(x-\bar{D}{\it\tau})$ , we obtain

(4.1) $$\begin{eqnarray}\bar{u}=\bar{D}+\sqrt{\bar{D}^{2}-\bar{{\it\lambda}}},\end{eqnarray}$$

where $\bar{D}$ is the speed of the wave and the bar denotes the steady state. It is seen that, for the solution to remain real at the end of the reaction zone, we must choose $\bar{D}\geqslant 1$ . The precise choice of the value is related to the degree of overdrive of the wave. Even though overdriven detonations can be included in the analysis, we focus on the important case of a self-sustained detonation in which the steady state has a sonic point at the end of the reaction zone. Then, $\bar{D}=1$ and

(4.2) $$\begin{eqnarray}\bar{u}=1+\sqrt{1-\bar{{\it\lambda}}},\end{eqnarray}$$

where $\bar{{\it\lambda}}$ solves

(4.3) $$\begin{eqnarray}\bar{{\it\lambda}}_{{\it\xi}}=-k(1-\bar{{\it\lambda}})\exp ({\it\theta}(\sqrt{q}\bar{u}+q\bar{{\it\lambda}}))\end{eqnarray}$$

with ${\it\xi}=x-{\it\tau}$ and boundary condition $\bar{{\it\lambda}}(0)=0$ (this is analogous to the ZND theory, see Fickett & Davis Reference Fickett and Davis2011). This solution is used in the analysis that follows.

The steady-state structure of the asymptotic model contains substantial information about the underlying modelling assumptions. A crucial qualitative feature of the structure is that, because of the behaviour of the leading-order correction to temperature, $T_{3}=u+q{\it\lambda}$ , it is possible to have a maximum of the reaction-rate function, (4.3), inside the reaction zone, as we show in figure 1. The presence of the internal maximum in the reaction rate appears to be a key factor responsible for the observed unsteady dynamics of real detonations, consistent with our previous work on a closely related model equation (Kasimov et al. Reference Kasimov, Faria and Rosales2013; Faria et al. Reference Faria, Kasimov and Rosales2014). In fact, we have verified that even the qualitative models of Fickett,

(4.4) $$\begin{eqnarray}\displaystyle u_{t}+\left(\frac{u^{2}}{2}+q{\it\lambda}\right)_{x}=0,\quad {\it\lambda}_{t}={\it\omega}({\it\lambda},u), & & \displaystyle\end{eqnarray}$$

and Majda (in the inviscid limit),

(4.5) $$\begin{eqnarray}\displaystyle (u+q{\it\lambda})_{t}+\left(\frac{u^{2}}{2}\right)_{x}=0,\quad {\it\lambda}_{t}={\it\omega}({\it\lambda},u), & & \displaystyle\end{eqnarray}$$

are capable of capturing the one-dimensional instabilities of detonation waves if the reaction rate is taken as (3.21). We recall that, in these models, the reaction rate usually takes the form

(4.6) $$\begin{eqnarray}{\it\omega}=(1-{\it\lambda}){\it\varphi}(u),\end{eqnarray}$$

where the ignition function, ${\it\varphi}(u)$ , is

(4.7) $$\begin{eqnarray}{\it\varphi}=\left\{\begin{array}{@{}ll@{}}{\it\varphi}_{0}(u),\quad & u>u_{i},\\ 0,\quad & u<u_{i},\end{array}\right.\end{eqnarray}$$

and $u_{i}$ is the ‘ignition-temperature’ parameter (e.g. Majda Reference Majda1980; Humpherys et al. Reference Humpherys, Lyng and Zumbrun2013). The new rate function, given by (3.21), reflects the important feature of the heat release in unstable detonations, which can exhibit a maximum inside the reaction zone. Thus, both Fickett’s and Majda’s models possess the necessary complexity needed to capture the qualitative dynamics of one-dimensional unstable detonations, provided that the reaction-rate function is chosen appropriately.

Figure 1.

Steady-state profiles and the rate function of the asymptotic solution for $q=1.7$ and ${\it\theta}=1.7$ .

The steady-state solution also provides a first quantitative test of the accuracy of the asymptotic approximation. In figure 2, we show a comparison between the ZND solutions of the reactive Euler equations and their asymptotic counterparts as predicted by the present theory. We see that the asymptotic approximation performs rather well when the heat release is small and the activation energy is large, i.e. when $Q\sim {\it\epsilon}$ and $E\sim 1/{\it\epsilon}^{2}$ . For example, for the realistic value of ${\it\gamma}=1.2$ $({\it\epsilon}=1/6)$ , the relative error is only a small percentage (figure 2 c). As expected, for the smaller value of ${\it\gamma}=1.1$ ( ${\it\epsilon}=1/11$ ), the agreement is seen to improve (figure 2 b,d), with the maximum relative error approximately 2 % of the flow velocity. As also expected, the approximation worsens as the values of ${\it\gamma}-1$ and $Q$ increase or the value of $E$ decreases.

Figure 2.

Comparison between the exact and asymptotic steady-state ZND profiles. The asymptotic solutions are calculated using (3.32)–(3.36). (a) Profiles for ${\it\gamma}=1.2,~Q=0.4,~E=50$ ; ${\it\epsilon}=1/6$ . (b) Profiles for ${\it\gamma}=1.1,~Q=0.4,~E=50$ ; ${\it\epsilon}=1/11$ . (c) Error for ${\it\gamma}=1.2,~Q=0.4,~E=50$ ; ${\it\epsilon}=1/6$ . (d) Error for ${\it\gamma}=1.1,~Q=0.4,~E=50$ ; ${\it\epsilon}=1/11$ .

4.2. Linear stability theory for the asymptotic model

With the steady-state solutions in reasonable agreement with the reactive Euler equations, we next investigate their stability properties. It is well known that the reactive Euler equations for the ideal-gas equation of state and simple-depletion Arrhenius kinetics predict that the steady-state detonations are unstable when the activation energy is large enough at a fixed heat release (Erpenbeck Reference Erpenbeck1964; Lee & Stewart Reference Lee and Stewart1990; Short & Stewart Reference Short and Stewart1998). In this section, we analyse the linear stability of the travelling wave solutions obtained in the previous section to see whether the asymptotic theory agrees with the Euler equations in this regard as well. We show that, indeed, the asymptotic steady-state solutions become unstable if either the heat release or the activation energy crosses a certain threshold. It is demonstrated that multidimensional effects play a non-trivial role in the onset of instabilities. We also note that the linearized equations below, (4.15)–(4.17), are the same as the linear equations analysed in Clavin & Williams (Reference Clavin and Williams2009).

To proceed with the analysis, let $\bar{u}$ and $\bar{{\it\lambda}}$ be the steady-state solution obtained in § 4.1. Rewriting (3.26)–(3.28) in a shock-attached frame, ${\it\chi}=x-s(y,{\it\tau})$ , where $s(y,{\it\tau})$ is the shock position, we obtain

(4.8) $$\begin{eqnarray}\displaystyle & u_{{\it\tau}}+(u-s_{{\it\tau}})u_{{\it\chi}}+{\textstyle \frac{1}{2}}{\it\lambda}_{{\it\chi}}+v_{y}-s_{y}v_{{\it\chi}}=0, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & u_{y}-s_{y}u_{{\it\chi}}-v_{{\it\chi}}=0, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{{\it\chi}}=-k(1-{\it\lambda})\exp ({\it\theta}(\sqrt{q}u+q{\it\lambda})). & \displaystyle\end{eqnarray}$$

Next, we expand the solution in normal modes,

(4.11) $$\begin{eqnarray}\displaystyle & u=\bar{u}({\it\chi})+{\it\delta}u_{1}({\it\chi})\exp ({\it\sigma}{\it\tau}+\text{i}ly)+O({\it\delta}^{2}), & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & v={\it\delta}v_{1}({\it\chi})\exp ({\it\sigma}{\it\tau}+\text{i}ly)+O({\it\delta}^{2}), & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & {\it\lambda}=\bar{{\it\lambda}}({\it\chi})+{\it\delta}{\it\lambda}_{1}({\it\chi})\exp ({\it\sigma}{\it\tau}+\text{i}ly)+O({\it\delta}^{2}), & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & s=\bar{D}{\it\tau}+{\it\delta}\exp ({\it\sigma}{\it\tau}+\text{i}ly), & \displaystyle\end{eqnarray}$$

and let ${\it\delta}\rightarrow 0$ . The linearized equations are then

(4.15) $$\begin{eqnarray}\displaystyle & (\bar{u}-\bar{D})u_{1}^{\prime }=-({\it\sigma}+\bar{u}^{\prime })u_{1}+{\it\sigma}\bar{u}^{\prime }-\text{i}lv_{1}-{\textstyle \frac{1}{2}}(g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1}), & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & v_{1}^{\prime }=\text{i}lu_{1}-\text{i}l\bar{u}^{\prime }, & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{1}^{\prime }=g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1}, & \displaystyle\end{eqnarray}$$

where the prime denotes a differentiation with respect to ${\it\chi}$ and

(4.18) $$\begin{eqnarray}\displaystyle & \displaystyle g({\it\chi})=-\frac{\partial {\it\omega}}{\partial u}(\bar{u},\bar{{\it\lambda}})=-k{\it\theta}\sqrt{q}(1-\bar{{\it\lambda}})\exp [{\it\theta}(\sqrt{q}\bar{u}+q\bar{{\it\lambda}})], & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle h({\it\chi})=-\frac{\partial {\it\omega}}{\partial {\it\lambda}}(\bar{u},\bar{{\it\lambda}})=-k[{\it\theta}q(1-\bar{{\it\lambda}})-1]\exp [{\it\theta}(\sqrt{q}\bar{u}+q\bar{{\it\lambda}})]. & \displaystyle\end{eqnarray}$$

The boundary conditions for (4.15)–(4.17) are obtained from linearizing (3.29)–(3.31):

(4.20a−c ) $$\begin{eqnarray}u_{1}(0)=2{\it\sigma},\quad v_{1}(0)=-2\text{i}l,\quad {\it\lambda}_{1}(0)=0.\end{eqnarray}$$

Noticing that for self-sustained detonation, $\bar{u}-\bar{D}\rightarrow 0$ as ${\it\chi}\rightarrow -\infty$ , we require that the right-hand side of (4.15) vanish in the limit as well, i.e.

(4.21) $$\begin{eqnarray}\displaystyle H({\it\sigma},l)=-({\it\sigma}+\bar{u}^{\prime })u_{1}+{\it\sigma}\bar{u}^{\prime }-\text{i}lv_{1}-{\textstyle \frac{1}{2}}(g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1})\rightarrow 0\quad \text{as }{\it\chi}\rightarrow -\infty .\quad & & \displaystyle\end{eqnarray}$$

Because $\bar{u^{\prime }}\rightarrow 0,~g({\it\chi})\rightarrow 0$ as ${\it\chi}\rightarrow -\infty$ , this solvability condition (alternatively called the ‘boundedness’ or the ‘radiation’ condition (Lee & Stewart Reference Lee and Stewart1990)) simplifies to

(4.22) $$\begin{eqnarray}-{\it\sigma}u_{1}-\text{i}lv_{1}-{\textstyle \frac{1}{2}}h(-\infty ){\it\lambda}_{1}\rightarrow 0\quad \text{as }{\it\chi}\rightarrow -\infty .\end{eqnarray}$$

To eliminate the numerical inconvenience of ${\it\sigma}=0,~l=0$ always being an eigenvalue – a consequence of the translation invariance of the travelling wave – we rescale $u_{1}$ and ${\it\lambda}_{1}$ by ${\it\sigma}$ and $v_{1}$ by $il$ . The stability problem is then posed as follows.

Solve

(4.23) $$\begin{eqnarray}\displaystyle & \displaystyle (\bar{u}-\bar{D})u_{1}^{\prime }=-({\it\sigma}+\bar{u}^{\prime })u_{1}+\bar{u}^{\prime }+\frac{l^{2}}{{\it\sigma}}v_{1}-\frac{1}{2}(g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1}), & \displaystyle\end{eqnarray}$$

(4.24) $$\begin{eqnarray}\displaystyle & \displaystyle v_{1}^{\prime }={\it\sigma}u_{1}-\bar{u}^{\prime }, & \displaystyle\end{eqnarray}$$

(4.25) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\lambda}_{1}^{\prime }=g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1}, & \displaystyle\end{eqnarray}$$

subject to $u_{1}(0)=2$ , $v_{1}(0)=-2$ and ${\it\lambda}_{1}(0)=0$ at the shock and the boundedness condition, (4.22), at negative infinity.

The preceding eigenvalue problem is solved numerically using the shooting method of Lee & Stewart (Reference Lee and Stewart1990). We solve the problem for different values of ${\it\theta}$ , $q$ and $l$ , which are the only remaining parameters. Consistent with the behaviour of the stability spectrum of detonation waves in the reactive Euler equations (Lee & Stewart Reference Lee and Stewart1990), we find that unstable modes do exist either for large enough $q$ or for large enough ${\it\theta}$ . We also find that the transverse modes, where $l\neq 0$ , tend to be more unstable than purely longitudinal disturbances (Short & Stewart Reference Short and Stewart1998).

Figure 3.

Contour plots of the spectral function (4.26) for ${\it\theta}=1.8$ and increasing $q$ : (a)  $q=1.8$ , (b) $q=2.2$ , (c) $q=2.6$ , (d) $q=2.75$ , (e) $q=3.5$ , (f) $q=4.5$ . It should be noted that the real part of the non-oscillatory root increases with $q$ . The red dots in the figures represent the eigenvalues. In (e) and (c), the dominant non-oscillatory root is not shown.

First, we consider the purely one-dimensional problem, i.e. with $l=0$ . In figure 3, we show the contour plot of the absolute value of the stability function,

(4.26) $$\begin{eqnarray}|H({\it\sigma},0)|=|-({\it\sigma}+\bar{u}^{\prime })u_{1}+{\it\sigma}\bar{u}^{\prime }-{\textstyle \frac{1}{2}}(g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1})|,\end{eqnarray}$$

as a function of the real and imaginary parts of ${\it\sigma}$ . The valleys in the plot of $|H({\it\sigma},0)|$ , which correspond to the darker regions in figure 3, provide an initial guess for the location of eigenvalues. A root solver is then used wherein the complex function, $H({\it\sigma},0)$ , is set to zero in order to accurately locate the eigenvalues. An increasing number of unstable eigenvalues is seen as the neutral boundary is crossed by increasing the heat release, $q$ , for a given value of the activation energy, ${\it\theta}$ . The qualitative behaviour of the spectrum is in agreement with that known for Euler detonations (Lee & Stewart Reference Lee and Stewart1990). Furthermore, figure 3(a–d) shows the migration of the oscillatory complex conjugates, with ${\it\sigma}_{i}\neq 0$ , into non-oscillatory unstable modes that are also observed in the Euler equations (Short Reference Short1997). Finally, in figure 3(e,f), we see that far into the unstable regime, many eigenvalues are found, indicating the increasing complexity of the linear spectrum.

We also obtain the neutral stability curves for the first two unstable eigenvalues in the asymptotic model, as shown in figure 4(a). It is seen that the lowest frequency mode 1 is more unstable than mode 2 for a wide range of $q$ and ${\it\theta}$ . Far away from the neutral boundary, the non-oscillatory root may become dominant, as is seen in the example shown in figure 3(c,d).

In figure 4(b), we plot the neutral curve for the lowest-frequency eigenvalue, $Q$ , indicated by the solid line in figure 4(a), in the plane of the heat release, $Q$ , and the activation energy, $E$ . The result is compared with the neutral curve computed directly from the reactive Euler equations (Lee & Stewart Reference Lee and Stewart1990). We observe that, as expected, the agreement improves with smaller $Q$ and larger $E$ .

Figure 4.

Neutral stability curves. (a) Neutral stability boundaries for the first two low-frequency modes in the asymptotic model. The modes are unstable to the right of the curves. (b) Comparison of the asymptotic and the exact (Lee & Stewart Reference Lee and Stewart1990) neutral stability boundaries for ${\it\gamma}=1.2$ , in the heat release, $Q$ , versus activation energy, $E$ , plane.

If $l\neq 0$ , then there is the possibility that transverse waves will trigger instabilities. This occurs in the Euler equations, where it is found that multidimensional instability prevails over purely longitudinal instability. Again, using the shooting method, we solve (4.23)–(4.25) numerically for various values of $l$ . Solving for the roots of the radiation function,

(4.27) $$\begin{eqnarray}H({\it\sigma},l)=-({\it\sigma}+\bar{u}^{\prime })u_{1}+{\it\sigma}\bar{u}^{\prime }-\text{i}lv_{1}-{\textstyle \frac{1}{2}}(g({\it\chi})u_{1}+h({\it\chi}){\it\lambda}_{1}),\end{eqnarray}$$

at ${\it\chi}=-\infty$ , we obtain the two-dimensional stability spectrum. We first fix $q=1.7$ and vary ${\it\theta}=1.65,\,1.60,\,1.55$ . In figure 5, we show the real and imaginary parts of the unstable modes for relatively small values of $l$ . It is seen that the model predicts purely two-dimensional instabilities for a certain range of the transverse wavenumbers, and that increasing the activation-energy parameter, ${\it\theta}$ , has a destabilizing effect on the steady-state solutions.

Figure 5.

Growth rate and frequency for the most unstable mode versus the wavenumber,  $l$ , for several values of the activation energy, ${\it\theta}$ .

We also perform a quantitative comparison between the two-dimensional stability of the asymptotic model and the known stability diagram for the Euler equations. We choose the values of ${\it\gamma}=1.2,~Q=0.4$ and $E=50$ . Then, after performing the appropriate conversion between dimensionless variables, we compare the asymptotic results with those obtained in Taylor, Kasimov & Stewart (Reference Taylor, Kasimov and Stewart2009) (see figure 6). We observe a fair agreement. There are, however, some differences. We see that for the parameters chosen in figure 6, the asymptotic model predicts one-dimensional instabilities, while the reactive Euler equations do not. Moreover, we observe that the disagreement between the imaginary parts of the eigenvalues increases with increasing wavenumber. This is expected because short transverse wavelengths cannot be represented accurately in the weak-curvature limit assumed for this model.

In the next section, we investigate the long-time nonlinear dynamics of the asymptotic solutions in regimes corresponding to linearly unstable steady-state one-dimensional solutions. The calculations are performed in both one and two spatial dimensions.

Figure 6.

Comparisons of the growth rate and frequency for the most unstable mode versus the wavenumber, for $Q=0.4$ and $E=50$ . The dashed curve corresponds to the model in this paper. The solid line corresponds to the reactive Euler equations as computed in Taylor et al. (Reference Taylor, Kasimov and Stewart2009).

5. Nonlinear dynamics of the asymptotic model

In the previous section, we showed that the asymptotic model exhibits the same linear stability behaviour as the reactive Euler equations. The question of what happens after the onset of instabilities is investigated through direct numerical simulations of the model equations. We show below that the travelling wave solutions of the asymptotic model reproduce, in both one and two spatial dimensions, the complexity observed in solutions of the reactive Euler equations.

5.1. Pulsating detonations

We focus here on the ability of the asymptotic model to predict the complex nonlinear dynamics of pulsating detonations. We perform a detailed numerical investigation of the large-time asymptotic behaviour of oscillatory solutions of the model. In the one-dimensional case, (3.26)–(3.28) reduce to

(5.1) $$\begin{eqnarray}\displaystyle & u_{{\it\tau}}+uu_{x}=-{\textstyle \frac{1}{2}}{\it\lambda}_{x}, & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{x}=-k(1-{\it\lambda})\exp ({\it\theta}(\sqrt{q}u+q{\it\lambda})). & \displaystyle\end{eqnarray}$$

The system (5.1)–(5.2) is equivalent to the one derived, using a somewhat different approach, in Clavin & Williams (Reference Clavin and Williams2002). No numerical calculations of either linear stability or nonlinear dynamics are reported in Clavin & Williams (Reference Clavin and Williams2002); however, theoretical arguments in that work suggest that the linear instability is present in the model. Our calculations in this section demonstrate that the system (5.1), (5.2) is faithful to the physics of pulsating detonations not only in terms of the linear instability, but also in terms of the dynamics of large-amplitude pulsating solutions. This system closely resembles the one derived in Rosales & Majda (Reference Rosales and Majda1983), with one crucial difference: as a consequence of the ${\it\gamma}-1=O({\it\epsilon})$ assumption, the reaction-rate function in (5.2) has a more complicated ${\it\lambda}$ dependence, which is in fact at the heart of the complexity of the solutions obtained here.

We solve (5.1), (5.2) numerically in a shock-attached frame (Kasimov & Stewart Reference Kasimov and Stewart2004), using a second-order finite volume scheme with a second-order total variation diminishing Runge–Kutta temporal discretization (LeVeque Reference LeVeque2002). Because no differentiation across the shock is performed, true second-order convergence is obtained for the cases tested, which include the convergence to various stable steady-state ZND solutions. We further verified the numerical algorithm by performing a cross-validation between the linear stability solver and the direct numerical solver. That is, we made sure that for several values of $q$ , the linear stability prediction of the neutral boundary agrees with the neutral stability boundary of the numerical scheme for the nonlinear model. For example, when $q=1.7$ , the linear stability curve in figure 4(a) indicates that ${\it\theta}_{c}=1.710$ is the neutral value of the activation energy, i.e. the ZND wave is unstable for ${\it\theta}>{\it\theta}_{c}$ and stable for ${\it\theta}<{\it\theta}_{c}$ . Numerical solutions of (5.1), (5.2) with ${\it\theta}$ slightly above and slightly below ${\it\theta}_{c}$ confirm this prediction, as seen in figure 7, where the shock state, $u_{s}$ , is plotted as a function of time.

Figure 7.

Nonlinear dynamics of (5.1), (5.2) near the neutral boundary, ${\it\theta}_{c}=1.710$ , as predicted by the linear stability theory for $q=1.7$ . (a) Stability of the ZND solution at ${\it\theta}=1.705<{\it\theta}_{c},~q=1.7$ . (b) Instability of the ZND solution and the limit-cycle attractor at ${\it\theta}=1.715>{\it\theta}_{c},~q=1.7$ . It should be noticed that the amplitude of the limit cycle scales roughly as the square root of the distance to the bifurcation point, as in a supercritical Hopf bifurcation.

We also compute solutions of (5.1), (5.2) further away from the neutral boundary in order to check whether the model captures a sequence of bifurcations leading to chaos, as occurs in the reactive Euler equations (Ng et al. Reference Ng, Higgins, Kiyanda, Radulescu, Lee, Bates and Nikiforakis2005; Henrick, Aslam & Powers Reference Henrick, Aslam and Powers2006). Such a sequence of bifurcations is indeed present in the model. Long-time simulations show that the solutions tend to either a fixed point, a limit cycle or (what appears to be) a chaotic attractor. We run the simulations at $q=5$ and plot the post-shock state, $u_{s}$ (by the Rankine–Hugoniot conditions, $u_{s}=2D$ , where $D$ is the shock speed), as a function of time for several different types of solution, as shown in figure 8. Beyond the stability boundary, the shock velocity becomes oscillatory. Near the neutral boundary, the oscillations have a small amplitude and are periodic (figure 8 a), but the structure of each period becomes more complex as we move away from the neutral boundary by increasing the activation energy (figure 8 b,c). Eventually, a value of ${\it\theta}$ is reached at which no obvious period is evident, as seen in figure 8(d).

Figure 8.

The shock state as a function of time for increasing values of ${\it\theta}$ and fixed $q=5$ showing pulsations of different complexity: (a) ${\it\theta}=0.75$ , (b) ${\it\theta}=0.8$ , (c) ${\it\theta}=0.83$ , (d) ${\it\theta}=0.85$ .

The behaviour described in the previous paragraph can be understood in terms of a period-doubling sequence of bifurcations where the period becomes longer and more complex with each bifurcation. In order to construct a bifurcation diagram illustrating this process, we proceed as follows. For each value of ${\it\theta}$ , we follow the evolution of the solution until it settles on the attractor. Then, we extract the set of local minima for $u_{s}({\it\tau})$ – a finite set for any periodic solution. The values of these minima are then plotted versus the bifurcation parameter, ${\it\theta}$ . The result is shown in figure 9, which is reminiscent of the standard Feigenbaum period-doubling cascade leading to chaos (Strogatz Reference Strogatz1994).

It is important to note that the further we move into the unstable region, the harder it is to numerically capture the wave dynamics with good accuracy. That is, in the highly unstable regime ( ${\it\theta}\gtrapprox 0.85$ in figure 9), the quantitative details of the bifurcation diagram are very sensitive to the grid resolution. In a truly chaotic regime, such sensitivity is intrinsic and reflects the nature of the system. However, another reason, which is at play here even before the apparently chaotic regime sets in, is that the wave dynamics can involve spatial scales that undergo large changes (by orders of magnitude) during the wave evolution. This is a direct consequence of the Arrhenius exponential dependence of the reaction rate, which can trigger large variations in the reaction rate from moderate changes in the temperature when the activation energy is large. This issue of stiffness associated with high-activation-energy detonations in unstable regimes is discussed in appendix B in more detail.

Because of the difficulties outlined here and in appendix B, without resorting to adaptive mesh refinement, 500–1000 points per half-reaction length are needed to obtain a bifurcation diagram with features that are grid-independent away from the chaotic regimes (see figure 9). While recognizing that such sensitivity to initial conditions or discretization errors is natural for chaotic dynamics, caution is still required when interpreting the results of numerical simulations of unstable detonations with high activation energies.

The appearance of fast temporal scales raises the question of whether or not the underlying assumption of slow evolution is still valid. Certainly, as the time scales become smaller, we should expect the asymptotic approximation to worsen. One can only hope, by comparison with solutions of the full system a posteriori, that the asymptotic theory still provides a good approximation. In the cases shown in figure 8, we observe the appearance of fast transients that correspond to the sharp spikes in the plots. These transients arise as a result of the wave steepening to a near-breaking state in the reaction zone. As a consequence, short spatial and temporal transients do occur, especially when the activation energy ${\it\theta}$ is large. Nevertheless, the observed complex dynamics are the properties of the reduced system, whether they are in precise quantitative agreement with the full system or not. Even in the presence of fast transients, at least a qualitative agreement with the full system is evident.

Figure 9.

The bifurcation diagram at $q=5$ showing the local minima of the attractor solution shock value, $u_{s}({\it\tau})$ , as a function of the activation energy, ${\it\theta}$ . The simulations were carried out at $N=15\,000$ (red) and $N=30\,000$ (black) grid points on the computational domain of length $L=30$ . Away from the chaotic regimes, the predictions are seen to be nearly identical.

5.2. Cellular detonation

Multidimensional instability is very important in gaseous detonations and results in cellular structures involving triple-point interactions on the detonation lead shock (Oran & Boris Reference Oran and Boris2001; Lee Reference Lee2008; Fickett & Davis Reference Fickett and Davis2011). It is therefore crucial to check whether the asymptotic model, (3.26)–(3.28), can reproduce the dynamics of not only one-dimensional, but also multidimensional, detonations. As we have seen in § 4.2, multidimensional instabilities can be dominant in the asymptotic model, with the dispersion relation showing a maximum growth rate for some non-zero transverse wavenumber $l$ (see figure 5). In this section, we calculate the long-term dynamics of detonation waves when two-dimensional effects are present. In particular, we show that the asymptotic model retains the essential complexity required to reproduce multidimensional cellular patterns. Solving (3.26)–(3.28) numerically turns out to be a non-trivial task and requires special care. Therefore, in the following subsection, we present a discussion of the algorithm employed in this work.

5.2.1. Numerical algorithm for the two-dimensional asymptotic system

In order to appreciate the subtlety associated with (3.26)–(3.28), we note that the equations comprise a nonlinear hyperbolic system with one of its characteristic planes orthogonal to time. This means the following.

1. (1) The initial data are given on a characteristic surface, the $x$ – $y$ plane. The absence of a time derivative in (3.27), (3.28) requires the initial conditions to satisfy $v_{x}=u_{y}$ and ${\it\lambda}_{x}={\it\omega}({\it\lambda},u)$ .

2. (2) Evolving in ${\it\tau}$ is a non-local procedure and, in the presence of a shock, care has to be taken to avoid spurious numerical oscillations.

Many existing numerical methods for solving (3.39), (3.40) are based on a formal rewriting of the system as a single equation,

(5.3) $$\begin{eqnarray}u_{x{\it\tau}}+\left(\frac{u^{2}}{2}\right)_{xx}+u_{yy}=0,\end{eqnarray}$$

by cross-differentiation and substitution. Two concerns arise with this approach. First, the validity of such a transformation is not obvious when $u$ and $v$ are discontinuous functions. Second, in the presence of chemical reactions, differentiation of (3.26) with respect to $x$ produces a delta forcing at the ignition-temperature locus due to the discontinuous nature of the reaction rate. Thus, techniques based on solving (5.3) are inadequate for our purposes.

Another common technique is to solve for the potential, ${\it\phi}$ , that satisfies

(5.4) $$\begin{eqnarray}{\it\phi}_{x{\it\tau}}+\left(\frac{{\it\phi}_{x}}{2}\right)_{x}^{2}+{\it\phi}_{yy}=0.\end{eqnarray}$$

Again, when chemical reactions are present, complications arise because the reaction rate now becomes an exponential function of ${\it\phi}_{x}$ , and the discretization errors in ${\it\phi}_{x}$ therefore exponentially amplify, requiring very-high-order methods for good accuracy.

The method we employ here is based on a direct semi-implicit discretization of (3.26), (3.27) following some of the ideas found in Hunter & Brio (Reference Hunter and Brio2000). In this method, all terms except for $v_{y}$ are treated explicitly. By choosing second-order spatial and temporal discretizations, we obtain an algorithm that is formally second order in time and space. The general procedure is outlined here to explain our reasoning for the choice of the algorithm.

Assuming that the solution at time ${\it\tau}={\it\tau}^{n}$ is known, we evolve it to ${\it\tau}={\it\tau}^{n+1}$ as follows.

1. (1) First, we employ a semi-implicit time discretization, where the $v_{y}$ term is treated implicitly. The motivation for this comes from the fact that some waves propagate infinitely fast in the $x$ – $y$ plane. In the simple case of a forward Euler time discretization, we obtain

   (5.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{u^{n+1}-u^{n}}{{\rm\Delta}{\it\tau}}+(F^{n}(u))_{x}+v_{y}^{n+1}={\it\lambda}_{x}^{n}, & \displaystyle\end{eqnarray}$$

   (5.6) $$\begin{eqnarray}\displaystyle & \displaystyle v_{x}^{n+1}=u_{y}^{n+1}, & \displaystyle\end{eqnarray}$$

   (5.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\lambda}_{x}^{n+1}={\it\omega}({\it\lambda}^{n+1},u^{n+1}), & \displaystyle\end{eqnarray}$$
   where $F(u)=u^{2}/2$ . A more quantitative reason for treating $v_{y}$ implicitly can be seen from the von Neumann stability analysis of the linearized system without chemical reactions, wherein a fully explicit scheme can be unstable (see appendix A for details).

2. (2) We approximate the explicit terms ${\it\lambda}_{x}^{n}$ and $F_{x}^{n}$ using a shock-capturing scheme, e.g. finite volume or essentially non-oscillatory/weighted essentially non-oscillatory algorithms (LeVeque Reference LeVeque2002).

3. (3) Using (5.6), we write a forward difference representation of $v_{i,j}^{n+1}$ in terms of $u_{i,j}^{n+1}$ and $u_{i^{\prime },j}^{n+1}$ for $i^{\prime }>i$ . For instance, a first-order forward difference scheme can be used, i.e.

   (5.8) $$\begin{eqnarray}v_{i,j}^{n+1}=v_{i+1,j}^{n+1}-{\rm\Delta}x(u_{i,j}^{n+1})_{y},\end{eqnarray}$$
   where the $y$ derivative approximation is postponed until the next step. The use of the forward difference here is a consequence of up-winding the infinitely fast waves propagating from right to left.

4. (4) We approximate the $y$ derivatives, e.g. by centred differences, to obtain the fully discrete scheme:

   (5.9) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i,j}^{n+1}=u_{i,j}^{n}+{\rm\Delta}{\it\tau}\left[-\frac{1}{2}({\it\lambda}_{x})_{i,j}^{n}-(F_{x})_{i,j}^{n}-\frac{1}{2{\rm\Delta}y}(v_{i,j+1}^{n+1}-2v_{i,j}^{n+1}+v_{i,j-1}^{n+1})\right], & \displaystyle\end{eqnarray}$$

   (5.10) $$\begin{eqnarray}\displaystyle & \displaystyle v_{i,j}^{n+1}=v_{i+1}^{n+1}-\frac{{\rm\Delta}x}{2{\rm\Delta}y}(u_{i,j+1}^{n+1}-u_{i,j}^{n+1}+u_{i,j-1}^{n+1}). & \displaystyle\end{eqnarray}$$

5. (5) In order to solve (5.9), (5.10), we sweep from right to left, assuming that the right boundary values of $u$ and $v$ are known at all times, as follows.

   1. (i) For a fixed $i$ , we insert (5.10) into (5.9) and solve the linear system for the vector $u_{i,j}^{n+1}$ .

   2. (ii) With $u_{i,j}^{n+1}$ , we find $v_{i,j}^{n+1}$ using (5.9).

   3. (iii) We repeat (i) and (ii) with $i=i-1$ until the left boundary is reached.

6. (6) Finally, we compute ${\it\lambda}^{n+1}$ by solving (5.7) with a boundary condition on the right end of the domain, which is given by

   (5.11) $$\begin{eqnarray}{\it\lambda}_{x}^{n+1}(x_{right},y)=0.\end{eqnarray}$$
   It should be noticed that (5.7) is actually an initial-value problem for ${\it\lambda}^{n+1}$ for fixed $y$ wherein $-x$ is a time-like direction. It can be solved with any initial-value solver (e.g. a Runge–Kutta method) if desired.

It may seem counterintuitive at first that the simplified model requires a semi-implicit method while the reactive Euler equations can be solved explicitly. The reason is that the asymptotic approximation is performed in a limit where the reactive Euler equations themselves would have to be treated implicitly. In order to understand this, we look at the three waves present in the Euler equations. Since the weak heat-release approximation implies that the detonation velocity is nearly acoustic, and an acoustic wave induces no flow behind it, we see that the speed of the forward acoustic characteristic, given by $u+c-D$ in a frame moving with the wave, is actually an $O({\it\epsilon})$ quantity ( $D\approx c$ and $u\approx {\it\epsilon}$ ). The entropy and backward acoustic characteristics, on the other hand, have $O(1)$ speeds. In the asymptotic model, a slow time, ${\it\tau}$ , is chosen so that the dynamics occurs on $O(1)$ time scales, and therefore some characteristics have speeds of $O(1/{\it\epsilon})$ in the slow time variable. When an explicit method is used to solve the Euler equations in this limit, a typical CFL (Courant–Friedrichs–Lewy) condition would require a time step, ${\rm\Delta}{\it\tau}\sim {\it\epsilon}{\rm\Delta}x$ , and as ${\it\epsilon}\rightarrow 0$ , it becomes clear that the time step restriction becomes unattainable and an implicit method is needed.

Figure 10.

Dynamics of the asymptotic model in channels of different widths: (a)  $L_{y}=2$ , (b) $L_{y}=10$ , (c)  $L_{y}=50$ , (d) $L_{y}=100$ . The plots show the asymptotic variable $u$ . The parameters are $q=1.7$ , ${\it\theta}=1.65$ , as in figure 5. The white region corresponds to the ambient state ahead of the wave, with $u=0$ .

The infinite characteristic velocity arises because the $x,y$ plane is a characteristic surface for the equations. A simple rotation in space–time ‘resolves’, to some extent, this issue such that in the ‘rotated’ space–time the system behaves as a standard system of conservation laws. The tradeoff for recovering a classical hyperbolic system with finite wave speeds is that one must now consider a grid with moving boundaries. This idea can be exploited to produce fully explicit schemes, as shown in Tabak & Rosales (Reference Tabak and Rosales1994) for the case without reactions.

5.2.2. Two-dimensional cellular detonations

Using the algorithm described in the previous section, we solve (3.26)–(3.28) in a frame moving with constant speed, $D_{0}=1$ . All simulations are initialized with the one-dimensional ZND solutions obtained in § 4.1. The computational domain in the $x$ direction in all of the simulations below is $[-35,5]$ , with inflow and outflow conditions used at the right and left ends of the domain respectively. We investigate, in the context of the asymptotic equations, (1) the effect of the width of the channel on the nonlinear stability properties of the travelling wave solutions, (2) the effect of the periodic boundary conditions and (3) the effect of increasing heat release on the wave dynamics.

In figure 10, we show the profile of $u$ at ${\it\tau}\approx 500$ for varying widths of the domain, $L_{y}$ , in the $y$ direction. The $y$ boundaries are modelled as rigid walls, appropriate for a detonation in a two-dimensional channel (i.e. a channel with negligible depth). The parameters are fixed at $q=1.7$ and ${\it\theta}=1.65$ , for which the linear stability calculation shows that the ZND wave is unstable only to two-dimensional perturbations. In a channel of width $L_{y}=2$ , the ZND wave remains stable, as seen in figure 10(a), as the spacing is too narrow for the two-dimensional instability to develop. Increasing the width of the channel, however, allows for a number of transverse modes to be excited, leading to the formation of a multidimensional cellular pattern, as shown in figure 10(c,d). The cell size appears to remain around 10 in all panels of figure 10(b–d). It should be noted that, for a channel with finite width, $L_{y}$ , the allowed transverse wavenumbers, $l$ , are given by $l={\rm\pi}n/L_{y}$ , where $n\in \mathbb{Z}$ . Therefore, only a discrete set of modes can be excited, and larger values of $L_{y}$ typically allow for more unstable modes to appear.

When we replace the solid-wall boundary conditions with periodic conditions, we observe a two-dimensional version of a spinning detonation. There is only one family of transverse shock waves which all propagate in the same direction. Such a wave can be imagined to form in a narrow gap between two concentric cylindrical tubes, as in a rotating detonation engine (Voitsekhovskii et al. Reference Voitsekhovskii, Mitrofanov and Topchian1966). A snapshot of the solution field, $u$ , is shown in figure 11, where relatively strong transversely propagating shocks can be seen.

Figure 11.

Detonation in a channel of width $L=100$ with periodic boundary conditions. The parameters are $q=1.7$ and ${\it\theta}=1.65$ .

Figure 12.

Dynamics of the asymptotic model for varying heat release at a fixed width $L=100$ and activation energy ${\it\theta}=1.65$ : (a) $q=1.7$ , (b)  $q=1.8$ , (c) $q=1.9$ , (d)  $q=2.5$ .

In figure 12, we show the effect of the heat release on the solution structure. We fix ${\it\theta}=1.7$ and increase the value of $q$ from 1.7 to 2.5. Consistent with the behaviour of detonations in the reactive Euler equations, regular cells are observed at small $q$ (figure 12 a,b), but with increasing $q$ , the structure of the detonation front becomes more complex with the formation of irregular cells (figure 12 d).

All of the previous results show that, at least at a qualitative level, the asymptotic model captures many important characteristics of multidimensional detonations. Next, we investigate how quantitatively close the predictions are to the solutions of the reactive Euler equations. The numerical simulations of the reactive Euler equations were carried out using PyCLAW (Ketcheson et al. Reference Ketcheson, Mandli, Ahmadia, Alghamdi, Quezada de Luna, Parsani, Knepley and Emmett2012), which provides, among other things, a Python wrapper for the classic routines in CLAWPACK. The software requires the user to provide a Riemann solver. We used a Roe-linearized Riemann solver with a Harten–Hyman entropy fix (LeVeque Reference LeVeque2002). The classic package of PyCLAW employs a second-order finite volume algorithm with a fractional step method for the time evolution (LeVeque Reference LeVeque2002). For the numerical simulations, the reactive Euler equations were non-dimensionalized in the conventional form by the ambient state with velocities scaled by $u_{a}=\sqrt{p_{a}/{\it\rho}_{a}}$ , spatial variables scaled by the half-reaction length, $x_{1/2}$ , and time by $u_{a}/x_{1/2}$ .

In figure 13, we show a comparison between the asymptotic solutions and the solutions of the reactive Euler equations. The parameters chosen are the same as in figure 6(b), i.e. ${\it\gamma}=1.2$ , $Q=0.4$ and $E=50$ , which in the asymptotic variables correspond to $q=2.4$ and ${\it\theta}=50/36\approx 1.389$ respectively. We start both simulations with the ZND solution and solve for a time interval large enough such that instabilities have already fully developed. It should be noticed that since the asymptotic model uses the slow time variable, ${\it\tau}$ , we only need to solve it for a relatively short time interval, ${\it\tau}\approx 100$ . The reactive Euler system, on the other hand, was non-dimensionalized in the conventional way using the regular time variable, $t$ , and thus we must solve the system up to $t\approx {\it\tau}/{\it\epsilon}=600$ . We perform all necessary scaling conversions between the variables so that the lengths and amplitudes shown in figure 13 are consistent. The length of the $y$ domain is fixed in the asymptotic model to be 20, which corresponds to a length of 27.83 in the regular dimensionless $y$ , while the $x$ domain is of length 40.

We see from figure 13 that the asymptotic model captures the salient features of multidimensional detonations with good quantitative agreement. That is, the characteristic scales of the detonation cells and the range of variation of the variable are seen to be close in both cases. A small difference in the transverse propagation velocity of the triple points is apparent, the speed being higher in the asymptotic case. The transverse shocks in the asymptotic solutions far from the lead shock are seen to be smoother than in the full solutions, which is due to a larger numerical diffusion in the algorithm for the asymptotic model.

Finally, we remark that weakly nonlinear detonations operate in the same fluid dynamic regime as weak shock wave focusing or weak shock reflection at near-grazing incidence. As pointed out by von Neumann (Reference von Neumann1963), in this limit the ‘classical’ Mach triple-point shock structure is impossible. Yet, both experiments and numerical calculations (Glaz et al. Reference Glaz, Colella, Glass and Deschambault1985) consistently exhibit triple shock structures in this regime. This apparent contradiction is known as the ‘von Neumann paradox’. Although the structures behind the lead shock in our calculations resemble the triple points observed in detonations (figure 14), a simple algebraic argument can show that the asymptotic equations do not admit ‘classical’ triple points, where three shocks separated by smooth flows meet (Tabak & Rosales Reference Tabak and Rosales1994; Hunter & Brio Reference Hunter and Brio2000). Thus, much like in the problem of weak shock reflection, the observed cells in the weak heat-release detonations considered in this work present yet another instance of the ‘von Neumann paradox’. Here, we make no attempt to address this problem and simply note that triple-point-like structures appear in our numerical simulations behind the lead shock; whether they are true discontinuities, sharp waves or some singularities remains an open problem.

Figure 13.

Comparison between the asymptotic and full solutions of the reactive Euler equations for a cellular detonation in a channel. (a) Density field, ${\it\rho}=1+{\it\epsilon}\sqrt{q}u$ , as predicted by the asymptotic model. (b) Density field computed from the reactive Euler equations.

6. Discussion and conclusions

In this work, we develop an asymptotic theory of multidimensional gaseous detonations within the framework of the compressible Navier–Stokes equations for a perfect gas reacting with a one-step heat-release law. The main outcome is a reduced model that consists, in two spatial dimensions, of a forced Burgers-like equation coupled with a heat-release equation and an equation that enforces zero vorticity. We have analysed it in detail and have shown that it captures many of the dynamical features of detonations in the reactive Euler equations. These include (1) the existence of steadily travelling waves (i.e. ZND solutions), (2) linear instability of the ZND solutions, (3) existence of a period-doubling bifurcation route to chaos in the one-dimensional case and (4) onset of cellular detonations in the two-dimensional case. Our theory considers weakly nonlinear detonations in the distinguished limit of small heat release, large activation energy and small ${\it\gamma}-1$ . In this limit, the dynamics occurs on time scales that are long compared with the scale of a characteristic chemical reaction time. Our reduced evolution system builds on the weakly nonlinear theory in Rosales & Majda (Reference Rosales and Majda1983) and Rosales (Reference Rosales1989) by adding the Newtonian approximation to it. In the one-dimensional case, our theory agrees with the model in Clavin & Williams (Reference Clavin and Williams2002).

Figure 14.

Triple-point-like structures in the asymptotic solution of the weakly nonlinear detonation, as shown by the magnitude of the density gradient.

As a consequence of the ${\it\gamma}-1=O({\it\epsilon})$ assumption, the temperature expansion in our analysis starts with an $O({\it\epsilon}^{2})$ correction to the leading-order term, as opposed to the $O({\it\epsilon})$ correction in Rosales & Majda (Reference Rosales and Majda1983) and Rosales (Reference Rosales1989). This is precisely what allows us to escape the fact that in Rosales & Majda (Reference Rosales and Majda1983) and Rosales (Reference Rosales1989) there must be leading-order corrections to temperature, velocity and density that behave in exactly the same way. The present scaling highlights an important difference between a reactive and a non-reactive shock, namely that the temperature in a reactive shock can increase at some distance from the shock front because of heat release. This increase in the post-shock temperature means that the energy release can possess a maximum inside the reaction zone. As we previously elucidated with a reactive Burgers model (Kasimov et al. Reference Kasimov, Faria and Rosales2013; Faria et al. Reference Faria, Kasimov and Rosales2014), the presence of such an internal maximum appears to be an important factor for the mechanism that generates resonances and subsequent amplification of the waves in the reaction zone.

It should be noted that all one-dimensional weakly nonlinear asymptotic models as well as the ad hoc models devised to capture detonation dynamics (Fickett Reference Fickett1979; Majda Reference Majda1980; Radulescu & Tang Reference Radulescu and Tang2011; Kasimov et al. Reference Kasimov, Faria and Rosales2013) are extensions of the Burgers equation modified by adding a reaction forcing term, with an extra equation for the reaction-progress variable. The reason is that the dynamics of any one-dimensional genuinely nonlinear hyperbolic system reduces to that of an inviscid Burgers equation in the weakly nonlinear limit (Rosales Reference Rosales1991).

While the reduced weakly nonlinear asymptotic theory has clear predictive power, as demonstrated in this work, recognizing its limitations is important and helps to highlight open problems that require further investigation. For example, at parameter values (especially of the heat release) corresponding to strong cellular and pulsating detonations typical in experiments and numerical simulations, the asymptotic predictions are not sufficiently accurate when compared with the reactive Euler equations. While it is not clear how to extend the present approach directly to strongly nonlinear detonations, the ability of the reduced system obtained in this work to capture essential dynamical characteristics of unsteady and multidimensional detonations serves as a strong indication that a theory of comparable simplicity may be possible for strongly nonlinear detonations.

The role of dissipative effects in weakly nonlinear combustion deserves some comment. Provided that $\mathit{Re}$ , $\mathit{Le}$ and $\mathit{Pr}$ are treated as parameters independent of ${\it\epsilon}$ and appropriate scaling assumptions are employed on their sizes, dissipative effects can be included without any essential modification into the asymptotic procedure outlined in § 3. In this case, terms proportional to $u_{xx}$ , $T_{xx}$ and ${\it\lambda}_{xx}$ must then be added to (3.26)–(3.28) (Faria Reference Faria2014). For example, if viscous effects only are considered (this requires appropriate conditions on the size of $\mathit{Le}$ , $\mathit{Pr}$ and $\mathit{Re}$ ), the result is a reduced model that differs from (3.26)–(3.28) only by an extra term on the right-hand side of (3.26). This term is ${\it\nu}u_{xx}$ , where ${\it\nu}=2/(3\sqrt{q}{\it\epsilon}\mathit{Re})$ is the dimensionless viscosity coefficient, which must be assumed to be $O(1)$ (Faria Reference Faria2014). This then endows the lead detonation shock with structure, just as in non-reactive shocks (Whitham Reference Whitham1974). However, for ideal gases, the assumptions used in this paper imply that $\mathit{Re}=\mathit{O}[\exp ({\it\theta}/{\it\epsilon}^{2})/{\it\epsilon}^{2}]$ (see item (5) at the beginning of § 3). Thus, $\mathit{Re}$ is exponentially large, so that the dissipative processes may be important only within exponentially thin transition layers. Whether or not such thin layers have an appreciable effect on the overall dynamics of detonations is not obvious, but recent studies seem to indicate this to be the case (Ziegler et al. Reference Ziegler, Deiterding, Shepherd and Pullin2011; Mazaheri, Mahmoudi & Radulescu Reference Mazaheri, Mahmoudi and Radulescu2012; Romick, Aslam & Powers Reference Romick, Aslam and Powers2012, Reference Romick, Aslam and Powers2015; Barker et al. Reference Barker, Humpherys, Lyng and Zumbrun2015). Although it is still possible to incorporate viscous effects into the model, the process is significantly more technical than the procedure presented in this paper, as it involves a singular perturbation problem in the presence of exponentially thin layers, and requires the use of matched asymptotic expansions.

Another problem that merits exploration, in our opinion, is that of the evolution of a small-amplitude localized initial perturbation, considered in the same weakly nonlinear asymptotic regime as in the present work. All of the underlying assumptions of the theory remain valid for this initial-value problem, which is relevant to the problem of detonation initiation (i.e. under what conditions the perturbation triggers a detonation wave).

From a mathematical point of view, the reduced system obtained in this work and its connection with the issue of the non-existence of triple-point shock structures for the Zabolotskaya–Khokhlov equation pose some puzzling and challenging problems.

Finally, we recall that the qualitative models introduced by Fickett and Majda do not contain instabilities with the state-dependent rate functions used in prior work with these models (e.g. ${\it\omega}=k(1-{\it\lambda}){\it\phi}(u)$ , Humpherys et al. (Reference Humpherys, Lyng and Zumbrun2013)). However, using the rate function in (3.28) (or, for that matter, the two-step rate function used in Radulescu & Tang (Reference Radulescu and Tang2011)), it can be shown to lead to the same complexity of the solutions as occurs in the reactive Euler equations. A modified Fickett’s analogue can be written as

(6.1) $$\begin{eqnarray}\displaystyle & u_{t}+uu_{x}=-{\textstyle \frac{1}{2}}{\it\lambda}_{x}, & \displaystyle\end{eqnarray}$$

(6.2) $$\begin{eqnarray}\displaystyle & {\it\lambda}_{t}=k(1-{\it\lambda})\exp ({\it\alpha}u+{\it\beta}{\it\lambda}), & \displaystyle\end{eqnarray}$$

${\it\alpha}$

${\it\beta}$

${\it\omega}={\it\psi}({\it\lambda}){\it\phi}(u)$

${\it\psi}$

${\it\phi}$