1 Introduction
The stability of two-dimensional shocks has been continually investigated since the pioneering studies of D’Yakov (Reference D’Yakov1954) and Kontorovich (Reference Kontorovich1957). They found that under certain conditions related to the Rankine–Hugoniot (RH) slope, the shock oscillates permanently, generating corresponding entropic and rotational perturbations downstream and emission of constant-amplitude sonic waves; the latter is commonly referred to as spontaneous acoustic emission (SAE). This regime is commonly referred to in the literature as the D’Yakov–Kontorovich (DK) instability.
The literature on this topic mostly focuses on adiabatic shocks moving into media that cannot be modelled as perfect gases (Freeman Reference Freeman1955; Zaidel Reference Zaidel1960; Erpenbeck Reference Erpenbeck1962; Fowles & Swan Reference Fowles and Swan1973; Fowles Reference Fowles1981; Fraley Reference Fraley1986; Bates & Montgomery Reference Bates and Montgomery2000; Wouchuk & López-Cabada Reference Wouchuk and López Cabada2004; Bates Reference Bates2004, Reference Bates2007, Reference Bates2012, Reference Bates2015; Konyukhov et al.
Reference Konyukhov, Likhachev, Fortov, Khishchenko, Anisimov, Lomonosov and Konyukhov2009; Wetta, Pain & Heuzé Reference Wetta, Pain and Heuzé2018; Touber & Alferez Reference Touber and Alferez2019) because adiabatic shocks in perfect gases are known to be inherently stable (Lapworth Reference Lapworth1959). When a shock moves steadily and is isolated from piston-like coupling effects, linear theory dictates that any perturbation to the shock surface decays on a long time scale in accordance with the power law
$t^{-3/2}$
when
where
${\mathcal{R}}_{d}=\unicode[STIX]{x1D70C}_{d}^{\prime }/\unicode[STIX]{x1D70C}_{o}^{\prime }$
is the mass-compression ratio across the shock and
${\mathcal{M}}_{d}=u_{d}^{\prime }/a_{d}^{\prime }$
is the Mach number of the shock front measured from the compressed gas reference frame. The quantity that distinguishes the type of asymptotic behaviour is
which measures the slope of the RH curve relative to the Rayleigh–Mickelson line. The subscripts
$d$
and
$o$
identify downstream and upstream properties, respectively, and the functions
$p^{\prime }$
and
$V^{\prime }$
define the dimensional pressure and the specific volume, respectively. Beyond the lower limit, namely,
the shock front enters the DK oscillatory mode, thereby generating constant-amplitude sonic waves downstream. It must be emphasized that unsteadiness in the base-flow variables, such as that found in blast waves (Vishniac & Ryu Reference Vishniac and Ryu1989), and piston-coupling effects that unavoidably appear to keep the adiabatic shock in the steady state (Bates Reference Bates2015) can modify these stability thresholds. Therefore, the application to the evolution of detonations induced by an instantaneous energy deposition is restricted to the long-time regime. In this type of configuration, the detonation strength diminishes whenever the front progresses up to find sonic conditions right behind the front. From that specific time, finite in cylindrical and spherical geometries, the detonation becomes effectively decoupled from downstream influences. The limits of validity of the planar and isolated shock assumptions in such conditions can be deduced in canonical configurations (see Levin & Chernyi Reference Levin and Chernyi1967; Liñan, Kurdyumov & Sánchez Reference Liñan, Kurdyumov and Sánchez2012; Clavin & Denet Reference Clavin and Denet2018). By properly including geometrical factors and coupling conditions, the model could also be adapted to address expanding accretion shocks (Velikovich et al. Reference Velikovich, Murakami, Taylor, Giuliani, Iwamoto and Wouchuk2016; Coughlin, Ro & Quataert Reference Coughlin, Ro and Quataert2019).
Following the ideas set forth by Bates & Montgomery (Reference Bates and Montgomery2000), Wouchuk & López-Cabada (Reference Wouchuk and López Cabada2004) and Bates (Reference Bates2012), who predicted the possibility of DK instability in gases governed by van der Waals forces, the present study shows that permanent DK oscillations and SAE can also occur in exothermic reactive shocks (Clavin & Searby Reference Clavin and Searby2016), such as detonations. The term ‘reactive shock’ is employed here to distinguish regular detonations, which are known to be very unstable (Zaidel Reference Zaidel1961; Erpenbeck Reference Erpenbeck1964; Fowles Reference Fowles1976; Majda & Rosales Reference Majda and Rosales1983; Bourlioux, Majda & Royburd Reference Bourlioux, Majda and Royburd1991, Reference Bourlioux, Majda and Royburd1992; Clavin & He Reference Clavin and He1996; Clavin, He & Williams Reference Clavin, He and Williams1997; Sharpe Reference Sharpe2001; Sharpe & Sharpe Reference Sharpe and Sharpe2001; Clavin & Denet Reference Clavin and Denet2002; Daou & Clavin Reference Daou and Clavin2003; Kasimov & Stewart Reference Kasimov and Stewart2004; Radulescu et al. Reference Radulescu, Sharpe, Law and Lee2007; Clavin & Williams Reference Clavin and Williams2009, Reference Clavin and Williams2012; Shepherd Reference Shepherd2009), from supersonic combustion fronts that may behave in a stable manner.
A linear normal-mode analysis, used to describe the multidimensional instability throughout the reaction zone in the detonation front (Short & Quirk Reference Short and Quirk1997; Short & Stewart Reference Short and Stewart1997, Reference Short and Stewart1999; Tumin Reference Tumin2007; Faria, Kasimov & Rosales Reference Faria, Kasimov and Rosales2015; Kabanov & Kasimov Reference Kabanov and Kasimov2018), provides a parametric domain in which planar detonation becomes unstable to initial perturbations. Three important bifurcation parameters were identified: the detonation overdrive, which measures how strong the piston-like driving effect produced by downstream conditions is with respect to the rocket-type effect produced by heat expansion; the amount of chemical heat released; and the chemical kinetics. The last of these can be conveniently characterized by the dimensionless activation energy for a one-step Arrhenius reaction model. Despite the wide variety of approaches used to describe the detonation behaviour, it can be concluded that the role of the acoustic field in a burnt gas is pivotal in triggering the instability. With regard to the perturbation wavelength, linear stability theory predicts unconditionally unstable detonations for short wavelengths within the high-activation-energy limit (Buckmaster & Ludford Reference Buckmaster and Ludford1988) and neutrally stable fronts for wavelengths much longer than the detonation thickness (Majda & Rosales Reference Majda and Rosales1983; Bourlioux et al. Reference Bourlioux, Majda and Royburd1992; Sharpe & Sharpe Reference Sharpe and Sharpe2001; Clavin & Williams Reference Clavin and Williams2012). Although most real detonations exhibit a multidimensional instability, theoretical analysis reveals that regular detonations are also unstable against planar disturbances (Lee & Stewart Reference Lee and Stewart1990; Sánchez et al. Reference Sánchez, Carretero, Clavin and Williams2001). In practice, these so-called galloping detonations are usually observed in small-diameter tubes as the wall constraints force the flow to be dominantly longitudinal (Fickett & Davis Reference Fickett and Davis2000; Gao, Ng & Lee Reference Gao, Ng and Lee2015).
For regular detonations, the amount of heat delivered to the fluid particles depends on the local thermodynamical properties, since the number of species involved can change with pressure and temperature conditions. Consequently, different fluid particles moving through the reaction zone may be subject to different energy changes. These variations, which can be incorporated in advanced numerical simulations capable to describe complex transport phenomena along with chemical kinetics and compressible effects (Sharpe & Sharpe Reference Sharpe and Sharpe2001; Martínez-Ferrer et al. Reference Martínez-Ferrer, Buttay, Lehnasch and Mura2014), are usually omitted in theoretical analyses that describe the reactive shock as a pure discontinuity with constant heat release. The aim of this work is to include the effect of variable heat release for an infinitely thin detonation, with the enthalpy changes being modelled as a shock-intensity-dependent function. The formulation, therefore, cannot be used to predict phenomena that are strongly bound to the shape of the internal profiles within the reaction zone.
With the strong simplification of a perfect-gas equation of state in mind, the results might also be applicable to thermonuclear detonations for which the number of nuclei burnt depends on the post-shock state (Gamezo, Khokhlov & Oran Reference Gamezo, Khokhlov and Oran2004). Nuclear detonations, as those occurring in Type Ia supernovae, can involve electron captures,
$\unicode[STIX]{x1D6FD}$
-decays, photo-disintegrations, two-body and three-body reactions (Ami Glasner et al.
Reference Ami Glasner, Livne, Steinberg, Yalinewich and Truran2018; Röpke & Sim Reference Röpke and Sim2018), thereby making the global heat release a shock-intensity-dependent function. For sufficiently strong shocks, further effects that can change the net energy gained by the fluid particles across the reacting shock might be dissociation, ionization or radiation, with the amount of energy absorbed depending on local opacity and post-shock density.
When the heat of reaction across the detonation front is positively correlated with the shock intensity, as when the heat released increases with pressure, density and/or temperature in a shocked gas, the detonation front may oscillate permanently. The possibility of strongly damped oscillations, which were previously associated with non-ideal gases (Bates Reference Bates2004), is also predicted when the exothermicity is negatively correlated with the shock strength. For sufficiently strong shocks, even when the sensitivity of the global heat released to the shock strength is very low, the shock dynamics has been found to switch to the permanent-oscillation regime, consequently causing SAE.
The methodology employs linear perturbation analysis in the thin-detonation limit or, equivalently, the long-wavelength regime. The theoretical model is then restricted to perturbation wavelengths
$\unicode[STIX]{x1D706}$
much longer than the detonation thickness
$\ell$
, which comprises the precursor adiabatic shock and the associated reaction layer. Following the analytical formulation employed in Huete, Sánchez & Williams (Reference Huete, Sánchez and Williams2013) and Huete et al. (Reference Huete, Jin, Martínez and Luo2017), the exact transient evolution of the detonation front and the associated disturbance patterns generated behind it, including explicit expressions for the distributions of density, pressure and vorticity, are provided. The acoustic energy flux generated behind the reactive shock is found to move nearly parallel to the shock surface. The isolated-shock assumption employed in this study, which places the shock acoustically far from piston-like supporting conditions, stems from the fact that the reacting gas is usually expanded to conditions that approach the sonic limit.
The study begins by formulating, in § 2, the equations that govern the detonation front and the burnt-gas flow. Next, § 3 shows the asymptotic long-time response of the detonation front, with the stability boundaries explicitly defined by the weakly and strongly overdriven limits. The initial value problem that determines the temporal evolution of the detonation front, along with the burnt-gas flow variables, is investigated by means of the Laplace transform in § 4, with the details of the mathematical treatment given in appendix A. Finally, concluding remarks are offered in § 5.
2 Problem formulation
2.1 Base-flow Rankine–Hugoniot equations
Consider a planar detonation front whose relative speed with respect to the upstream flow is
$u_{o}^{\prime }>a_{o}^{\prime }$
, where the speed of sound
$a_{o}^{\prime }$
is modelled with the perfect-gas assumption
$a_{o}^{\prime 2}=\unicode[STIX]{x1D6FE}p_{o}^{\prime }/\unicode[STIX]{x1D70C}_{o}^{\prime }$
. The changes in the flow variables across the front are obtained by integrating the conservation equations in the streamwise direction, yielding
provided that exothermic effects occur within a thin layer right behind the shock wave. Ideally, the total energy per unit mass delivered to the fluid particles is given by the difference between the formation enthalpies of the products and reactants
$\unicode[STIX]{x0394}h_{f}^{o}$
.
Although the heat release rate is strongly affected by local properties behind the precursor shock wave, mainly the temperature, the overall heat released is slightly affected by the thermodynamic conditions in the reaction layer. Nonetheless, sufficiently strong changes in pressure or temperature could promote the formation of different radicals, in which case the global thermodynamic equilibrium involving a combination of different species could affect the overall amount of heat released. One example is nitrogen in air, which cannot be considered inert when the temperature is sufficiently high. Direct computations of one-dimensional detonations for stoichiometric methane–air mixtures, performed independently with models (McBride, Zehe & Gordon Reference McBride, Zehe and Gordon2002; Browne, Ziegler & Shepherd Reference Browne, Ziegler and Shepherd2008), show that temperature tends to reduce global heat release while pressure contributes with an opposite effect. Other phenomena beyond chemical reactions may affect the overall heat released, such as the final state (liquid or gas) of the H
$_{2}$
O, or radiation losses that depend on the temperature. In particular, phase changes can affect the heat released when the fuel is a diluted spray (Williams Reference Williams1961, Reference Williams1985; Smolders & van Dongen Reference Smolders and van Dongen1992; Guha Reference Guha1994; Martínez-Ruiz et al.
Reference Martínez-Ruiz, Urzay, Sánchez, Liñan and Williams2013; Smirnov et al.
Reference Smirnov, Betelin, Kushnirenko, Nikitin, Dushin and Nerchenko2013; Vimercati, Gori & Guardone Reference Vimercati, Gori and Guardone2018). Along with these applications, the astrophysical context opens a wide range of possibilities for non-adiabatic waveforms, such as the thermonuclear detonations formed in Type Ia supernovae whose potential nuclear fuel depends on the density, which peaks just behind the precursor shock (Khokhlov Reference Khokhlov1989; Khokhlov, Oran & Wheeler Reference Khokhlov, Oran and Wheeler1997; Gamezo et al.
Reference Gamezo, Khokhlov and Oran2004; Wright et al.
Reference Wright, Nagaraj, Kneller, Scholberg and Seitenzahl2016; Ami Glasner et al.
Reference Ami Glasner, Livne, Steinberg, Yalinewich and Truran2018; Röpke & Sim Reference Röpke and Sim2018).
When the detonation is driven by reactions that depend on local thermodynamic functions, the heat released depends on the shock strength. Certainly, any particular case would call for a specific model to compute the variations in the global heat released with the shock intensity. Because this dependence is presumably small, a linear perturbation analysis could be employed to write
as the variable heat of reaction delivered to the fluid particles. Reactions that can depend on the bulk density, such as those involved in thermonuclear detonations, can be easily incorporated into (2.3) with the corresponding equation of state.
The heat released is conveniently scaled with the fresh-gas properties. Therefore,
${\mathcal{H}}_{0}=(\unicode[STIX]{x1D6FE}^{2}-1)\unicode[STIX]{x0394}h_{f}^{o}/(2a_{o}^{\prime 2})$
defines the dimensionless order-of-unity heat released. Because the pressure and temperature scale with the shock strength squared for sufficiently strong shocks, such as those expected in detonations, variations of the global heat released with the thermodynamic variables can be collected in the small term
$\unicode[STIX]{x1D716}$
, where
defines the shock-dependent dimensionless heat released. The ratio
${\mathcal{M}}_{o}=u_{o}^{\prime }/a_{o}^{\prime }$
defines the Mach number of the detonation wave, and
$q$
is the constant heat released contribution, i.e. when it does not depend on the shock strength. Note that the factor
$\unicode[STIX]{x1D716}\ll 1$
can be either positive, when the heat released is positively correlated with the shock strength, or negative, when it is negatively correlated.
Simple manipulation of (2.1a )–(2.1c ) provides the jump equations across the reaction front, namely,
for density and pressure, respectively. The Mach number of the burnt-gas flow is
with the function
accounting for the dimensionless non-adiabatic parameter
${\mathcal{H}}_{0}$
. For non-reacting shock waves,
$\unicode[STIX]{x1D705}=1-{\mathcal{M}}_{o}^{-2}$
, which simplifies (2.5)–(2.7) to the well-known adiabatic RH relationships. Note that the form of (2.5)–(2.8) is similar to those found in any classical reference on combustion (see Williams Reference Williams1985; Clavin & Searby Reference Clavin and Searby2016); the only difference here is the way the function
${\mathcal{H}}_{0}$
is modelled in (2.4).
The so-called Chapman–Jouguet (CJ) regime defines the minimum detonation strength associated with a given energy release, namely,
where
${\mathcal{M}}_{cj}^{\unicode[STIX]{x1D716}=0}=\sqrt{1+q}+\sqrt{q}$
is the CJ propagation Mach number when the heat released does not depend on the shock intensity,
$\unicode[STIX]{x1D716}=0$
. In this regime, the flow behind the reacting shock is sonic,
${\mathcal{M}}_{d}=1$
, and the detonation front decouples from the downstream condition; thereby, its motion is self-sustained by the expansion of the products resulting from the exothermic transformation. As expected, a positive correlation of the heat release with the shock intensity increases the self-propagating velocity
${\mathcal{M}}_{cj}$
, while the contrary occurs for
$\unicode[STIX]{x1D716}<0$
.
The RH curve relates the pressure jump to the mass compression ratio (or specific volume) across the detonation wave. This function,
is explicitly written for convenience in defining the auxiliary function (1.2) when it is written in terms of dimensionless variables; see (2.16).
The RH curves are plotted in figure 1. The intersection of the (straight and thin) Rayleigh–Mickelson lines with the RH curves provides the final state after the detonation front. For sufficiently strong waves under constant-heat conditions,
$\unicode[STIX]{x1D716}=0$
, the exothermic effect is negligible, and the corresponding curves approach the adiabatic RH curve. For
$\unicode[STIX]{x1D716}\neq 0$
, however, as the energy released to the fluid particles scales with the shock strength squared, the effect is always noticeable for strong shocks. It is readily seen that positively correlated heat releases shift the curves upwards, whereas the opposite applies for negative values of
$\unicode[STIX]{x1D716}$
. Although the shape of the curve is not significantly affected, this minor change induced by
$\unicode[STIX]{x1D716}$
may have a pivotal impact on the perturbation dynamics, as will be shown in § 3.
Rankine–Hugoniot curves for exothermic (coloured lines) and adiabatic (black line) shocks. Positive and negative correlations are represented by dashed and dotted lines, respectively. The straight narrow black lines are Rayleigh–Mickelson lines.

2.2 Linear perturbation analysis
Assuming that the background-flow variables are constant, which restricts the analysis to the stable detonation fronts commonly associated with overdriven detonations and/or low activation energies, regular perturbation analyses can be performed for the detonation front and the burnt-gas flow. Then, an initial ripple of the form
$\unicode[STIX]{x1D713}_{d}(t=0)=\unicode[STIX]{x1D713}_{d0}\cos (k^{\prime }y^{\prime })$
is considered, where
$\unicode[STIX]{x1D713}_{d0}$
is the initial amplitude and
$k^{\prime }=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}$
is the associated wavenumber. Together with the stability prerequisite in the base-flow variables, linear theory and the fast-reaction-limit assumptions set the following limits on validity:
$\ell \ll \unicode[STIX]{x1D713}_{d0}\ll \unicode[STIX]{x1D706}$
, with the corresponding lengths sketched in figure 2.
Sketch of the corrugated detonation front and the representative lengths (a), with
$\ell \ll \unicode[STIX]{x1D713}_{d0}\ll \unicode[STIX]{x1D706}$
. Qualitative description of the shock-strength sensitivity effect (b).

The wavenumber is employed to scale the spatial variables,
$x=k^{\prime }x^{\prime }$
and
$y=k^{\prime }y^{\prime }$
, and time is measured by the downstream speed of sound,
$\unicode[STIX]{x1D70F}=a_{d}^{\prime }k^{\prime }t$
. The dimensionless order-of-unity detonation ripple is defined as
$\unicode[STIX]{x1D709}_{d}=\unicode[STIX]{x1D713}_{d}/\unicode[STIX]{x1D713}_{d0}$
. The functions of interest are perturbed, and the small factor
$\unicode[STIX]{x1D713}_{d0}k^{\prime }$
is used to scale the amplitude of the perturbations, namely,
Integration domain and distinguished regions. Blue dotted lines refer to the particle paths, grey dashed lines indicate the positive and negative characteristic paths, respectively, and the red solid line defines the detonation trajectory.

The problem reduces to that of integrating the linearized Euler equations or, equivalently, the wave equation (2.13) for
$\unicode[STIX]{x1D70F}\geqslant 0$
over the domain delimited by the leading reflected sonic wave travelling backwards,
$x=-\unicode[STIX]{x1D70F}$
, and the reaction front moving forwards,
$x={\mathcal{M}}_{d}\unicode[STIX]{x1D70F}$
, as sketched in figure 3. One boundary condition is provided by the isolation assumption, which translates into neglect of the effect of the acoustic waves reaching the detonation from behind; this is detailed mathematically in the discussion after (A 7) in appendix A. Note that this condition does not restrict the existence of right-travelling waves in the burnt gas, i.e. waves propagating in the detonation propagation direction, but it impedes their ability to reach the detonation front from behind because they are all generated at the detonation front and there is no reflecting boundary downstream. The other boundary condition is given by the flow properties at the detonation front,
for the pressure, density and velocity perturbations, respectively, which are given by the linearized reactive RH equations, namely,
Note that, unlike the one-dimensional RH equations, the detonation ripple function, which is to be determined, calls for an extra independent equation: the momentum conservation equation in the transverse direction, which ultimately translates into the conservation of the tangential velocity.
The mass compression factor
${\mathcal{R}}_{d}$
and the post-detonation Mach number
${\mathcal{M}}_{d}$
are given in (2.5) and (2.7) as functions of the gas polytropic index
$\unicode[STIX]{x1D6FE}$
, the upstream Mach number
${\mathcal{M}}_{o}$
and the non-adiabatic contribution
${\mathcal{H}}_{0}$
. The factor that accounts for the RH slope downstream, (1.2), which is found to govern the stability of the reactive shock, is
and relates to the so-called DK parameter through
$h=-\unicode[STIX]{x1D6E4}_{d}$
(see D’Yakov Reference D’Yakov1954; Kontorovich Reference Kontorovich1957). It is readily seen that this factor tends to
${\mathcal{M}}_{o}^{-2}$
in the adiabatic limit,
$q\ll 1$
, to yield a perfect-gas equation of state.
Manipulation of the linearized RH equations leads to the following system of equations for the detonation boundary conditions:
and
to be employed when integrating the Euler equations for the burnt gas. The factors accompanying the pressure perturbation in (2.17) and (2.18) are conveniently expressed it terms of the functions
3 Long-term dynamics
3.1 General asymptotic theory
The qualitative response of the detonation front can be inferred with the aid of figure 2, where the shock-strength dependence of the heat released is depicted for positive and negative values of
$\unicode[STIX]{x1D716}$
. Positive correlations with the shock strength, which, in the linear regime, scales with the growth of the detonation shape as
$\unicode[STIX]{x1D6FF}{\mathcal{M}}_{o}=\dot{\unicode[STIX]{x1D709}}_{d}{\mathcal{M}}_{o}/({\mathcal{R}}_{d}{\mathcal{M}}_{d})$
, translate into amplified expansion of the hot products within the reaction zone when the detonation accelerates,
$\dot{\unicode[STIX]{x1D709}}_{d}>0$
, and a diminished expansion-driven push when it decelerates,
$\dot{\unicode[STIX]{x1D709}}_{d}<0$
. Generally, for very small values of
$\unicode[STIX]{x1D716}$
, this positive feedback decreases the damping of the oscillation amplitude towards a planar shape. It is shown later that, for sufficiently large values of
$\unicode[STIX]{x1D716}$
, yet still much smaller than unity,
$\unicode[STIX]{x1D716}\sim {\mathcal{M}}_{o}^{-2}\ll 1$
, the detonation oscillates with constant amplitude in the long term. This positive feedback that heat release may exert on the detonation dynamics resembles the Rayleigh criterion for the thermoacoustic instability. The opposite case,
$\unicode[STIX]{x1D716}<0$
, has a stabilizing counter-effect, with the endothermic perturbation decreasing the detonation ripple’s growth, thereby promoting faster decay towards the final planar shape.
The quantitative response that sets the limits of the distinguished asymptotic regimes is less straightforward. For the detonation to spontaneously and constantly radiate acoustic waves downstream, two conditions must be fulfilled simultaneously: the reactive shock must oscillate with a non-decaying amplitude, and the frequency of the oscillations must be sufficiently high. The latter condition is equivalent to saying that the acoustic wavenumber vector in the detonation reference frame points backwards. This condition is easily deduced by anticipating that sonic waves travelling in the burnt-gas region are functions of
$(\unicode[STIX]{x1D701}_{a}\unicode[STIX]{x1D70F}-k_{a}x)$
. The dimensionless frequency
$\unicode[STIX]{x1D701}_{a}$
and wavenumber
$k_{a}$
are obtained from the adiabatic dispersion relation
$\unicode[STIX]{x1D701}_{a}^{2}=k_{a}^{2}+1$
, as derived from (2.13), along with the compatibility condition at the reaction front,
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D701}_{a}-{\mathcal{M}}_{d}k_{a}$
, which provides
and
Acoustic perturbations become evanescent when the radicand in (3.1) and (3.2) becomes negative, i.e. when the detonation frequency satisfies
$\unicode[STIX]{x1D714}<\sqrt{1-{\mathcal{M}}_{d}^{2}}$
. In contrast, constant-amplitude waves occur when
$\unicode[STIX]{x1D714}\geqslant \sqrt{1-{\mathcal{M}}_{d}^{2}}$
, with the sign of
$k_{a}$
determining the propagation direction of the travelling waves. Positive values of
$k_{a}$
, or
$\sqrt{1-{\mathcal{M}}_{d}^{2}}<\unicode[STIX]{x1D714}<1$
, refer to waves travelling to the right, whereas the opposite describes waves travelling to the left,
$\unicode[STIX]{x1D714}>1$
. Then, the singular case
$\unicode[STIX]{x1D714}=1$
corresponds to standing acoustic waves in the longitudinal direction.
The Laplace transforms of the detonation ripple and pressure disturbances, given in (A 10) and (A 11), respectively, provide the following dispersion relationship:
with the nature of the poles describing the long-term response (Clavin & Williams Reference Clavin and Williams2012; Clavin & Searby Reference Clavin and Searby2016).
The value of
$\unicode[STIX]{x1D70E}_{b}$
relative to
$\unicode[STIX]{x1D70E}_{c}$
, which is shown in (A 12), defines the distinguished limits. The acoustic nature of the burnt-gas perturbation field is responsible for the term
$\sqrt{s^{2}+1}$
in (3.3), which introduces a branch cut in the complex plane of the Laplace variable
$s$
. As a direct consequence, not all purely imaginary poles that are roots of this equation are translated into undamped oscillations on the long time scale (Clavin & Searby Reference Clavin and Searby2016).
The condition that sets the limits for stable oscillations is
$s=\pm \text{i}$
, which occurs when
$\unicode[STIX]{x1D70E}_{b}=\unicode[STIX]{x1D70E}_{c}$
, as easily verified from (3.3). The dispersion relation is then recast as
$\sqrt{s^{2}+1}(s+\unicode[STIX]{x1D70E}_{b}\sqrt{s^{2}+1})=0$
, which, in the temporal domain, translates into damped oscillations whose asymptotic decay rate is proportional to
$\unicode[STIX]{x1D70F}^{-1/2}$
. For adiabatic shocks in perfect gases, this singular condition is only met in the strong-shock limit (Fraley Reference Fraley1986). Permanent oscillations at the reaction front occur when
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
, that is, when the imaginary poles in (3.3) lie outside the branch cut. The asymptotic detonation-oscillation frequency is
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D701}_{-}\sqrt{1-{\mathcal{M}}_{d}^{2}}$
, where
tends to unity in the limit
$\unicode[STIX]{x1D70E}_{b}-\unicode[STIX]{x1D70E}_{c}\ll 1$
.
The DK stability limit depends on the local slope of the RH curve in the burnt-gas state, as noted in § 1. This is commonly presented as a function of the parameter
which separates the acoustically radiating regime from the non-radiating regime. The equation
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
is equivalent to
$\unicode[STIX]{x1D70E}_{b}=\unicode[STIX]{x1D70E}_{c}$
, which is derived above.
When
$\unicode[STIX]{x1D70E}_{b}>\unicode[STIX]{x1D70E}_{c}$
, the detonation front oscillates towards the asymptotic planar solution with an amplitude that decays in time in accordance with the power law
$\unicode[STIX]{x1D70F}^{-3/2}$
. However, two distinguished scenarios exist: when
$\unicode[STIX]{x1D70E}_{b}>\unicode[STIX]{x1D70E}_{c}+1/(4\unicode[STIX]{x1D70E}_{c})$
(or
$\unicode[STIX]{x1D6EC}<0$
, as defined in Bates (Reference Bates2004)), the approach towards the asymptotic decay rate occurs faster than it does under regular conditions, namely,
$\unicode[STIX]{x1D70E}_{c}+1/(4\unicode[STIX]{x1D70E}_{c})>\unicode[STIX]{x1D70E}_{b}>\unicode[STIX]{x1D70E}_{c}$
. Regular conditions are associated with finite-strength shocks moving adiabatically in perfect gases. The limit that distinguishes fast or regular damping can also be expressed as
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
, where
A supersonic front becomes unstable when
$\unicode[STIX]{x1D6E4}_{d}<-(1+2{\mathcal{M}}_{d})$
or
$\unicode[STIX]{x1D70E}_{b}<-1$
. For this condition to be met, the RH slope must reverse, which does not occur in detonations moving through perfect gases. Within the parametric domain that defines strength-sensitive detonation,
${\mathcal{M}}_{o}\sim \sqrt{q}$
,
$q\sim 1{-}10$
, and
$|\unicode[STIX]{x1D716}|\ll 1$
, four regimes are identified and delimited in figure 4. The colour code used here is the same as that used in the computations of the asymptotic temporal evolution in figure 9.
General stability limits for detonation waves (Clavin & Williams Reference Clavin and Williams2012; Clavin & Searby Reference Clavin and Searby2016).

The above-mentioned parameters are computed in figure 5 for
$\unicode[STIX]{x1D716}>0$
(a) and
$\unicode[STIX]{x1D716}<0$
(b) as functions of
${\mathcal{M}}_{o}$
and
$q$
. Solid lines correspond to
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
, dashed lines to the limit
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
and dotted lines to
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
. The region above the solid lines is discarded because it corresponds to
${\mathcal{M}}_{o}<{\mathcal{M}}_{cj}$
. The region below the dashed/dotted lines corresponds to the regular oscillating regime, which is similar to that found in adiabatic shock waves. It is found that only positively correlated shock-strength sensitivities (
$\unicode[STIX]{x1D716}>0$
) lead to constant-amplitude oscillations, which correspond to the region enclosed by the lines
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
and
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
. In contrast, strongly damped oscillations are restricted to negative shock-strength sensitivities (
$\unicode[STIX]{x1D716}<0$
), with the region bounded by the curves
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
and
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
being the parameter space associated with this regime. For positively correlated shock-strength sensitivities (figure 5
a), it is found that the point of detachment occurs before order-of-unity energy variations, as observed for
$\unicode[STIX]{x1D716}=0.005$
and
${\mathcal{M}}_{o}\sim 4$
, which corresponds to energy changes of less than 10 % of the global heat released. For sufficiently weak detonations, the energy change induced by the shock strength is negligible; therefore, the curves
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
and
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
approach the boundary
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
. For finite shock strengths, however, the shock-strength sensitivity is no longer negligible, and the curves detach from their respective boundaries, enabling the possibility of permanent and strongly damped oscillation regimes for
$\unicode[STIX]{x1D716}>0$
and
$\unicode[STIX]{x1D716}<0$
, respectively. For negatively correlated cases, the shock strength needed to switch to strongly damped oscillations is higher, making this condition harder to meet. Results agree with previous analysis developed for stable detonations, i.e. characterized by a sufficiently low activation energy, in the long-wavelength limit (Bourlioux et al.
Reference Bourlioux, Majda and Royburd1992; Clavin & Williams Reference Clavin and Williams2012; Clavin & Searby Reference Clavin and Searby2016).
Dynamic detonation regimes for positively (a) and negatively (b) correlated shock-strength sensitivities as functions of
${\mathcal{M}}_{o}$
and
$q$
. Solid lines correspond to values of
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
. The curves have been computed with
$\unicode[STIX]{x1D6FE}=1.3$
.

The foregoing stability limits could have been obtained omitting the previous analysis on the poles in the complex plane. Direct inspection of (1.1) and (1.3) presented in the introduction, with the function
$\unicode[STIX]{x1D6E4}_{d}$
being already given in (2.16) and (2.10), provides the stability limits of the planar reacting front. Besides, previous works have been devoted to investigate the stability limits of planar reactive shocks and they offer a comprehensive and exhaustive analysis on the topic, including the extension to nonlinear perturbations (Majda & Rosales Reference Majda and Rosales1983; Clavin & Searby Reference Clavin and Searby2016). It is, however, convenient to repeat them here in the terms presented above, since part of the analysis is useful to describe the transient evolution of the reactive front and the properties of the burnt-gas variables. Moreover, the position of the poles in the complex plane is paramount to determine the perturbation decay rates and to obtain the corresponding strongly and weakly overdriven limits addressed below.
3.2 Strongly and weakly overdriven detonation limits for
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
The transition to acoustically radiating detonations can be investigated in the limits of strongly,
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}\gg 1$
, and weakly,
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1\ll 1$
, overdriven detonations. As inferred from figure 5(a), the iso-curve
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
shows that the value of
$q$
required to switch to the oscillatory mode decreases when the Mach number is sufficiently large and increases.
Performing a regular expansion for large Mach numbers,
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}\gg 1$
, and small sensitivity parameters,
$\unicode[STIX]{x1D716}\ll 1$
, yields the following expression for the iso-curve
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
:
which shows that the limiting Mach number increases as
$q$
and
$\unicode[STIX]{x1D716}$
decrease. It is also found that the curve
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
approaches the asymptote
$q=0$
for any
$\unicode[STIX]{x1D716}>0$
. For the singular case of
$\unicode[STIX]{x1D716}=0$
, the condition
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
is only met for
${\mathcal{M}}_{o}\gg 1$
, as found for adiabatic shocks in the strong-shock limit (Fraley Reference Fraley1986), which leaves no room for stable oscillations.
Weakly overdriven conditions are evaluated by the functions
$\unicode[STIX]{x1D70E}_{b}$
and
$\unicode[STIX]{x1D70E}_{c}$
. Regular expansions in the limit
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1\ll 1$
reveal that
$\unicode[STIX]{x1D70E}_{b}$
approaches unity in the form
$1+O(\sqrt{{\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1})$
, whereas the function
$\unicode[STIX]{x1D70E}_{c}$
converges to a finite value
$\unicode[STIX]{x1D70E}_{c}^{\ast }+O(\sqrt{{\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1})$
, where
The limiting condition that separates stable and decaying oscillations
$\unicode[STIX]{x1D70E}_{b}=\unicode[STIX]{x1D70E}_{c}$
provides
as the value above which the reacting shock undergoes stable oscillations in the CJ regime. As before, the small sensitivity parameter in (3.9) appears as a singularity in the leading-order term. Then, for reactive shocks whose sensitivity to the heat released is strictly null,
$\unicode[STIX]{x1D716}=0$
, the oscillatory mode is infeasible.
The strongly and weakly overdriven limits, which are found to be independent of the adiabatic index
$\unicode[STIX]{x1D6FE}$
, are computed in figure 6(a) as a function of
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1$
for the conditions in figure 5(a). Figure 6(a) shows that weakly overdriven detonations call for a minimum amount of heat to be released to achieve the stable oscillatory mode; this critical value is inversely proportional to
$\sqrt{\unicode[STIX]{x1D716}}$
. Similarly, strongly overdriven detonations have critical Mach numbers that are inversely proportional to
$\unicode[STIX]{x1D716}$
.
3.3 Strongly and weakly overdriven detonation limits for
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
The strongly overdriven limit
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}\gg 1$
corresponding to the condition
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
is
to leading order, which presents a pole at
$q=(1-\unicode[STIX]{x1D6FE}^{2/3})/(4\unicode[STIX]{x1D716})$
. Unlike the radiating condition shown in (3.7), the iso-curve
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
approaches a constant corresponding to the minimum heat of reaction needed to enter the strongly damped oscillation regime.
Weakly and strongly overdriven limits for
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
and
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
for
$\unicode[STIX]{x1D6FE}=1.3$
.

The opposite weakly overdriven limit for low sensitivities is
provided that
$1=\unicode[STIX]{x1D70E}_{c}^{\ast }+1/(4\unicode[STIX]{x1D70E}_{c}^{\ast })$
in the CJ regime, with
$\unicode[STIX]{x1D70E}_{c}^{\ast }$
being given in (3.8). Because
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{dam}$
is a monotonically decreasing function of
${\mathcal{M}}_{o}/{\mathcal{M}}_{cj}-1$
, any amount of dimensionless heat released greater than that in (3.11) makes the detonation oscillations attenuate quickly. Unlike the value obtained in (3.9), the factor
$q\sim 1/\unicode[STIX]{x1D716}$
, which demands that much more heat be released to reach the strongly damped regime under CJ conditions.
Both limits are computed in figure 6(b) for the conditions in figure 5(b). It is found that small deviations in the overdrive with respect to the CJ regime quickly decrease the critical heat release required to enter the strongly damped oscillatory mode. Convergence to the limiting value (3.11) is very slow. When
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
, linear perturbation theory in the fast-reaction limit predicts that
$\unicode[STIX]{x1D6E4}_{d}=\unicode[STIX]{x1D711}_{rad}$
is unconditionally satisfied for any value of
$q$
, which makes the front oscillations decay with
$\unicode[STIX]{x1D70F}^{-1/2}$
.
Note that the case
${\mathcal{M}}_{o}={\mathcal{M}}_{cj}$
is singular, i.e. any positive deviation of the velocity field with respect to the base flow may generate a locally supersonic zone, and the same applies to negative temperature perturbations. In principle, the detonation travelling is restricted to an asymmetric propagation condition: the detonation can move faster than the CJ velocity but not slower. When analysed in the strong-shock limit, both post-detonation gas velocity and hot-gas speed of sound scale linearly with the shock intensity,
$u_{d,cj}\sim a_{d,cj}\sim u_{o,cj}\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D6FE}-1)$
, thereby allowing both positive and negative small variations of the propagating velocity, with the governing constraint on the sonic condition downstream (
$M_{d}=1$
) unaltered. However, this adds a restriction to the perturbation amplitudes in the weakly overdriven limit, as the reacting flow cannot be expanded supersonically downwards. A proper bound on the fluctuation amplitudes in the weakly overdriven limit calls for a nonlinear treatment of the perturbation analysis, such as the one developed by Majda & Rosales (Reference Majda and Rosales1983). It must also be noted that the speed of sound is taken as an average property of the burnt gas; different values of the speed of sound may coexist in gas mixtures.
4 Exact temporal evolution
4.1 Transient evolution of the detonation front
There exist various analytical techniques to compute the temporal evolution of the reactive shock in the linear perturbation regime. Once the Laplace transform of the perturbations at the detonation front is explicitly determined, as shown in (A 10) and (A 11), with use made of Cauchy’s integral theorem and anticipating the presence of poles that lead to constant-amplitude detonation oscillations, the inverse Laplace transform can be employed to obtain the temporal evolution of the detonation front. Some examples of this type can be found in Wouchuk, Huete Ruiz de Lira & Velikovich (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009), Huete, Wouchuk & Velikovich (Reference Huete, Wouchuk and Velikovich2011, Reference Huete, Wouchuk and Velikovich2012) and Huete et al. (Reference Huete, Sánchez and Williams2013, Reference Huete, Jin, Martínez and Luo2017) for the evolution of adiabatic and detonation fronts moving through non-uniform flows, respectively. An alternative technique is the use of Bessel functions, which naturally appear as solutions of the hyperbolic sound wave equation (2.13).
By choosing the latter technique, the temporal evolution of the pressure disturbances at the reacting shock and the amplitude of the front corrugation are written as
and
respectively, where the values
$N_{\unicode[STIX]{x1D708}}$
and
$D_{\unicode[STIX]{x1D708}}$
are provided in (A 15) and (A 18).
Although the methods of the inverse Laplace transform and the series of Bessel functions are equivalent, from the point of view of their limits of validity, the latter presents an advantage: the recurrence expression that provides the coefficients accompanying the Bessel functions in (4.1) and (4.2) – see appendix A – does not depend on the type of asymptotic regime. We refer to Bates (Reference Bates2004) for a complete description of the family of solutions found in isolated shocks when the inverse Laplace transform method is employed.
Computations for corrugation amplitude
$\unicode[STIX]{x1D709}_{d}(\unicode[STIX]{x1D70F})$
are displayed in figure 7 for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D6FE}=1.3$
and a negative correlation parameter
$\unicode[STIX]{x1D716}=-0.01$
(figure 7
a), and for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D6FE}=1.3$
and positive correlation parameters
$\unicode[STIX]{x1D716}=0.06$
(figure 7
b),
$\unicode[STIX]{x1D716}=0.037$
(figure 7
c) and
$\unicode[STIX]{x1D716}=0$
(figure 7
d). It is found that the highly damped oscillating regime,
$\unicode[STIX]{x1D70E}_{b}>\unicode[STIX]{x1D70E}_{c}+1/(4\unicode[STIX]{x1D70E}_{c})$
, is obtained in figure 7(a). The case of regular detonations, which is found to lie within the domain
$\unicode[STIX]{x1D70E}_{c}>\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}+1/(4\unicode[STIX]{x1D70E}_{c})$
, is displayed in figure 7(b). In these two regimes, the amplitude of the detonation ripple oscillations decays with time with the power law
$\unicode[STIX]{x1D70F}^{-3/2}$
. The strong attenuation of the detonation-front vibrations in figure 7(a) is associated with the negative correlation with the shock intensity. The higher the shock perturbation towards the fresh mixture, the stronger is the mechanical energy subtracted from the bulk and the weaker is the overall heat release.
Temporal evolution of the detonation-front ripple amplitude
$\unicode[STIX]{x1D709}_{d}(\unicode[STIX]{x1D70F})$
computed for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D6FE}=1.3$
(
$\unicode[STIX]{x1D716}=-0.01$
), and for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D6FE}=1.3$
(
$\unicode[STIX]{x1D716}=0$
,
$0.037$
,
$0.06$
).

Note that the asymptotic decay rate of an isolated Bessel function is
$\text{J}_{\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D70F}\gg 1)\sim \unicode[STIX]{x1D70F}^{-1/2}$
, not the power law
$\unicode[STIX]{x1D70F}^{-3/2}$
that the complete series in (4.2) displays for adiabatic shocks. For
$\unicode[STIX]{x1D70E}_{c}>\unicode[STIX]{x1D70E}_{b}$
, a condition that applies for the two first cases (figures 7
a and 7
b), the leading terms cancel and the remaining leading order dominates the asymptotic expansion (Fraley Reference Fraley1986). When the condition
$\unicode[STIX]{x1D70E}_{b}=\unicode[STIX]{x1D70E}_{c}$
is satisfied, a case computed in figure 7(c), the leading terms do not vanish, which renders a decay rate with the power law
$\unicode[STIX]{x1D70F}^{-1/2}$
for
$\unicode[STIX]{x1D70F}\gg 1$
. For
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
, there exists a constant-amplitude contribution added to the time-decaying function. This is readily observed in figure 7(d), where the oscillations at the front reach a non-decaying regime for sufficiently long times, i.e. the DK regime.
A clear distinction of the family of solutions distinguished in figure 4 is also provided in Bates (Reference Bates2004), where the inverse Laplace transform technique renders different late-time asymptotic functions for
$\unicode[STIX]{x1D709}_{d}(\unicode[STIX]{x1D70F})$
depending on the equations of state. The highly damped oscillating regime is found to present an exponential contribution that decays faster, and a power-law contribution of the form
$\unicode[STIX]{x1D70F}^{-3/2}$
that dominates the long-time state. In particular, we refer to equations (57)–(60) in Bates (Reference Bates2004), which were contrasted against FAST2D hydrodynamical simulations and that have been employed here to check the results given by (4.2). To the best of the authors’ knowledge, similar computations of the type considered here have not been performed.
The expressions that describe the long-time evolution of the detonation front in the DK regime can be obtained by direct inspection of the residues in (A 10) and (A 11), which yields
and
for the detonation ripple and pressure perturbations, respectively, with the asymptotic amplitude of the detonation shape deviations with respect to the planar shape being
The natural oscillation frequency in the temporal variable
$\unicode[STIX]{x1D70F}$
is
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D701}_{-}\sqrt{1-{\mathcal{M}}_{d}^{2}}$
, where the value of
$\unicode[STIX]{x1D701}_{-}$
is provided in (3.4). The asymptotic long-time oscillation amplitude
$\unicode[STIX]{x1D709}_{d}^{\infty }$
and frequency
$\unicode[STIX]{x1D701}_{-}$
are computed in figures 8(a) and 8(b), respectively, for
$\unicode[STIX]{x1D6FE}=1.3$
and
$\unicode[STIX]{x1D716}=0.02$
. It is seen that
$\unicode[STIX]{x1D70E}_{b}=\unicode[STIX]{x1D70E}_{c}$
yields
$\unicode[STIX]{x1D701}_{-}=1$
or
$\unicode[STIX]{x1D714}=\sqrt{1-{\mathcal{M}}_{d}^{2}}$
corresponding to sound waves running parallel to the detonation surface. It is also observed that, although the amplitude grows with the heat release, the oscillation frequency barely changes since
$\unicode[STIX]{x1D701}_{-}$
is found to be roughly unity within the domain
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
explored in figure 8. Employing other values of the same order for
$\unicode[STIX]{x1D716}$
does not produce any qualitative change in the results.
Asymptotic detonation ripple amplitude (a) and oscillation frequency (b) as a function of
${\mathcal{M}}_{o}$
and
$q$
. The curves have been computed with
$\unicode[STIX]{x1D6FE}=1.3$
and
$\unicode[STIX]{x1D716}=0.02$
.

A simple form to inspect the power law that describes the long-time asymptotic decay is to display the peak values of the pressure function in a log–log scale for a sufficiently long temporal domain, and then to compute the corresponding regression lines. This analysis is performed in figure 9 with the same input parameters as those employed in figure 7. The results confirm the previous discussion on the asymptotic long-time limit. It is, however, instructive to check that some configurations, such as that found for
$\unicode[STIX]{x1D716}=0$
, demand very long times to fit on the asymptotic trend. In some conditions, the shock must cover distances of the order of 100 perturbation wavelengths. This fact limits the reproducibility of the asymptotic analysis in real configurations, where base-flow changes may occur in the formation and propagation of the reactive front. Besides, the long-time dynamics computed here is based on the isolated-front assumption. It is known, however, that piston coupling effects can effectively modify the stability limits obtained in figure 5, as shown in Freeman (Reference Freeman1955) and Bates (Reference Bates2012, Reference Bates2015). Although the isolated assumption is certainly non-physical for a non-reactive steady shock, detonations can propagate with a high degree of self-sustainability, or can even be totally self-sustained in the CJ regime. The closer to this regime the front is, the weaker are the influences from the downstream boundary conditions, as the coupling is determined by the acoustic field that becomes sonic in CJ conditions.
Long-time pressure perturbation peaks computed for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D6FE}=1.3$
(
$\unicode[STIX]{x1D716}=0.06$
,
$0.037$
,
$0$
), and for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D6FE}=1.3$
(
$\unicode[STIX]{x1D716}=-0.01$
).

Note that the scaling chosen to define the dimensionless spatial and temporal variables is not appropriate to reproduce CJ conditions, where the downstream flow is strictly sonic,
${\mathcal{M}}_{d}=1$
. Since the oscillation frequency factor in (4.4) is
$\sin (\unicode[STIX]{x1D701}_{-}r)=\sin (\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F})$
in the CJ regime, the value
$\unicode[STIX]{x1D714}$
approaches zero as
$\unicode[STIX]{x1D714}\sim 1-{\mathcal{M}}_{d}^{-2}$
. In other words, pressure perturbations cannot travel in the streamwise direction in sonic conditions, so the acoustic coupling of the pressure field at different positions along the moving detonation front is impossible. It is therefore convenient to derive an asymptotic theory for weakly overdriven detonations such as that provided in § 3.2.
4.2 Burnt-gas flow variables
The pressure field in the burnt-gas flow can be derived with the aid of (A 13) and undoing the hyperbolic variable change carried out in (A 1), namely
$\unicode[STIX]{x1D712}=\tanh ^{-1}(x/\unicode[STIX]{x1D70F})$
and
$r=\sqrt{\unicode[STIX]{x1D70F}^{2}-x^{2}}$
, to give
The spatial distribution of the pressure field is computed in figure 10 for the same conditions as in figure 7 for
$\unicode[STIX]{x1D70F}=40$
(orange),
$\unicode[STIX]{x1D70F}=80$
(brown) and
$\unicode[STIX]{x1D70F}=120$
(black). Two regimes are well distinguished: acoustically radiating and non-radiating conditions. The SAE condition is determined by the condition
$\unicode[STIX]{x1D701}_{-}>1$
, or equivalently
$\unicode[STIX]{x1D714}>(1-{\mathcal{M}}_{d})^{1/2}$
, a condition that is always satisfied in the stable oscillatory mode
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
. As previously mentioned, the case
$k_{a}<0$
represents acoustic waves travelling in the direction of the post-shock flow, while the waves moving in the opposite direction correspond to
$k_{a}>0$
. The solution corresponding to shock oscillation frequency
$\unicode[STIX]{x1D714}=1$
is associated with standing acoustic waves that separate the solution travelling to the left,
$\unicode[STIX]{x1D714}>1$
, from the one travelling to the right,
$(1-{\mathcal{M}}_{d}^{2})^{1/2}<\unicode[STIX]{x1D714}<1$
, in the reference frame comoving with the burnt-gas fluid.
Spatial vorticity field
$\hat{p}$
for
$\unicode[STIX]{x1D70F}=40$
(orange),
$\unicode[STIX]{x1D70F}=80$
(brown) and
$\unicode[STIX]{x1D70F}=120$
(black) for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D716}=-0.01$
(a), for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D716}=0$
(b), and for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D716}=0.06$
(c). The grey dashed line shows the pressure-detonation history.

Sufficiently far from the shock in the burnt gas, i.e. at a distance much longer than the perturbation wavelength, the asymptotic pressure perturbations when the reactive shock oscillates with constant amplitude are given by
where
$\hat{p}_{d}^{\infty }$
is the amplitude of the pressure perturbations at the detonation front. Besides the initial-value problem resolved in this work, normal-mode analysis can also be conducted to obtain the asymptotic flow-field variables in the burnt gas.
Associated with the SAE, there exists a sound energy radiation flux that reads
$\boldsymbol{j}=\unicode[STIX]{x1D70C}_{d}^{\prime }a_{d}^{\prime 2}|\hat{p}|^{2}\boldsymbol{v}_{a}^{\prime }$
in the burnt-gas reference frame. The acoustic wavenumber vector that points in the sound compression direction,
$\boldsymbol{v}_{a}^{\prime }=a_{d}^{\prime }\boldsymbol{k}_{a}$
, is
with the characteristic angle
$\unicode[STIX]{x1D703}_{a}$
being sketched in figure 11(a).
Computations of the acoustic energy flux can be employed to characterize the noise associated with the unsteadiness of the detonation, as done for regular shocks (Ribner Reference Ribner1969). In applications where the detonation wave remains anchored to the domain boundaries, as they are for supersonic combustors that employ detonations as a way of burning, the analysis is conveniently done in the detonation wave reference frame (see figure 11
b). For that case, the sound energy flux is given by
$\boldsymbol{v}_{ad}^{\prime }=a_{d}^{\prime }\boldsymbol{k}_{a}-u_{d}^{\prime }\hat{\boldsymbol{e}}_{x}$
.
Focused on the result of the most interest, the longitudinal projection is calculated scaled with the zero-order energy flux
$\unicode[STIX]{x1D70C}_{d}^{\prime }a_{d}^{\prime 3}$
, thereby yielding
as the dimensionless acoustic radiation, where the corresponding emission angle
is determined by the vectors
$\hat{\boldsymbol{e}}_{x}$
and
$\boldsymbol{v}_{ad}^{\prime }/|\boldsymbol{v}_{ad}^{\prime }|$
. The acoustic energy flux adds second-order corrections to the density, velocity and pressure behind the detonation wave, as well as mass, momentum and energy fluxes. Such deviations affect the overall velocity propagation, since the RH relations get modified by these second-order corrections (Velikovich, Huete & Wouchuk Reference Velikovich, Huete and Wouchuk2012).
Sketch of the sound-wave emission angle when measured in the burnt-gas reference frame (a) and detonation-front reference frame (b).

It is immediate to see that
$\unicode[STIX]{x1D701}_{-}\gtrsim 1$
implies
$\cos \unicode[STIX]{x1D703}_{a}\sim {\mathcal{M}}_{d}$
and
$\unicode[STIX]{x1D703}_{d}\sim \unicode[STIX]{x03C0}/2$
. In view of the computations displayed in figure 8, this is a representative condition within the oscillatory mode,
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
. Therefore, for scenarios where the rippled detonation oscillates with a non-decaying amplitude, i.e. when the overdrive is sufficiently low and
$\unicode[STIX]{x1D716}>0$
, the pressure field is dominated by transverse waves moving parallel to the detonation surface. Longitudinal acoustic radiation
$\unicode[STIX]{x1D703}_{a}\sim \unicode[STIX]{x1D703}_{d}\sim \unicode[STIX]{x03C0}$
would occur when the detonation front experiences high-frequency oscillations
$\unicode[STIX]{x1D701}_{-}\gg 1$
, such as those occurring in forced-oscillating conditions imposed by external excitation (Huete et al.
Reference Huete, Sánchez and Williams2013, Reference Huete, Jin, Martínez and Luo2017).
The isentropic density and temperature variations induced by the acoustic radiation are simply given by
$\hat{\unicode[STIX]{x1D70C}}_{a}(x,y,\unicode[STIX]{x1D70F})=\hat{p}(x,y,\unicode[STIX]{x1D70F})$
and
$\hat{T}_{a}(x,y,\unicode[STIX]{x1D70F})=(\unicode[STIX]{x1D6FE}-1)\hat{p}(x,y,\unicode[STIX]{x1D70F})$
, respectively, with the pressure function
$\hat{p}$
being previously given in (4.6) and displayed in figure 10.
The breakdown of the velocity field into acoustic and rotational contributions, namely
$\hat{u} (x,y,t)=\hat{u} _{a}(x,y,t)+\hat{u} _{r}(x,y)$
and
$\hat{v}(x,y,t)=\hat{v}_{a}(x,y,t)+\hat{v}_{r}(x,y)$
for the streamwise and lateral components, respectively, facilitates the problem resolution. As readily seen from the Euler equations, the acoustic contribution of the velocity fluctuations obeys
for the longitudinal and transverse components, respectively. It is, however, simpler to withdraw the acoustic velocity perturbations from the momentum conservation equations in (2.12d ). On the other hand, the steady-rotational perturbations, which are isobaric in the linear limit, are governed by
where the function
represents the dimensionless vorticity field. In the absence of diffusive effects, the vorticity–entropy perturbations in the burnt gas can only be generated at the detonation front.
In many applications, the upstream gas is not uniform and contains density and velocity disturbances (Huete et al. Reference Huete, Sánchez and Williams2013, Reference Huete, Jin, Martínez and Luo2017). For such cases, the vorticity deposited in the burnt gas is affected by other contributions, as the amplification of the vorticity ahead of the front might be as a direct result of the overall compression ratio and the baroclinic generation. The latter occurs when there are density and/or pressure gradients in the direction transverse to the detonation propagation direction.
Similar to the velocity perturbations, the total density field comprises acoustic and entropic contributions,
$\hat{\unicode[STIX]{x1D70C}}(x,y,t)=\hat{\unicode[STIX]{x1D70C}}_{a}(x,y,t)+\hat{\unicode[STIX]{x1D70C}}_{e}(x,y)$
. The latter can be readily obtained from the linearized RH equation (2.15c
) after subtracting the acoustic part, yielding
The vorticity and entropy perturbations generated in the far field can be determined by the asymptotic detonation evolution. While the planar detonation does not generate any vorticity–entropy perturbation in the long-time regime for
$\unicode[STIX]{x1D70E}_{b}>\unicode[STIX]{x1D70E}_{c}$
, there exist vortical and entropic disturbances in the burnt gas for
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
, which take the following form:
and
respectively.
Isocontours of the vorticity field
$\unicode[STIX]{x1D6FA}$
for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D716}=-0.01$
(a), for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D716}=0$
(b), and for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$\unicode[STIX]{x1D716}=0.06$
(c).

The spatial distribution of the vorticity left in the burnt-gas particles is computed in figure 12 for
${\mathcal{M}}_{o}=7$
,
$q=15$
,
$\unicode[STIX]{x1D716}=-0.01$
(a), and for
${\mathcal{M}}_{o}=3$
,
$q=1$
,
$0$
(b) and
$\unicode[STIX]{x1D716}=0.06$
(c). Figure 12(a,b) displays a vorticity field that decays with the distance
$x$
, as the detonation pressure decays with time as shown in figure 10. Figure 12(c), corresponding to DK conditions,
$\unicode[STIX]{x1D70E}_{b}<\unicode[STIX]{x1D70E}_{c}$
, shows a vorticity field that renders order-of-unity values of
$\unicode[STIX]{x1D6FA}(x\gg 1)$
, in consonance with figure 10. The maximum and minimum values are located at
$y=\pm \unicode[STIX]{x03C0}/2$
(
$\sin y=\pm 1$
) for
$-\unicode[STIX]{x03C0}<y<\unicode[STIX]{x03C0}$
, as they correspond to positions where the tangential velocity is the highest. The qualitative picture of the steady-density field is similar to that represented in figure 12 for the vorticity field, but shifted
$\unicode[STIX]{x03C0}/2$
in the transverse direction (see the dependence on the transverse coordinate) and with the iso-curve scaling being modified according to the factors accompanying the function
$\hat{p}_{d}$
in (4.13) and (4.14).
5 Conclusions
Planar reactive shocks, whose overall heat released is invariant to the shock strength, are stable to long-wavelength perturbations (Majda & Rosales Reference Majda and Rosales1983; Clavin & Williams Reference Clavin and Williams2012). Under these conditions, the initial corrugation of a sufficiently overdriven detonation decays with time following the power law
$t^{-3/2}$
, as occurs with regular shocks moving in perfect gases. In contrast, when the heat released is positively correlated with the shock strength, the so-called DK regime can be reached, which leads to constant-amplitude oscillations of the detonation front on a long time scale. Consequently, non-evanescent acoustic waves are spontaneously radiated downstream, and entropic rotational perturbations are deposited in the burnt gas. The positive feedback that heat release exerts on the detonation dynamics, when it is in phase with the pressure changes, presents a clear resemblance to the Rayleigh criterion for the thermoacoustic instability. The acoustic energy flux generated directly behind the detonation wave is dominated by transverse waves that move parallel to the detonation surface. The boundary between the two regimes corresponds to a singular oscillating regime in which perturbations at the detonation front decay with time following the power law
$t^{-1/2}$
, as found for adiabatic shocks in the strong-shock limit and under CJ conditions for planar reactive shocks (Fraley Reference Fraley1986; Clavin & Williams Reference Clavin and Williams2012). This boundary is presented as closed-form expressions for the limits of strongly and weakly overdriven detonations. When the heat released is negatively correlated with the shock strength, the reverse behaviour is found for the oscillating detonation dynamics. In this case, the heat released decreases with the detonation ripple growth rate, which acts as a stabilizing mechanism. Thus, the initially corrugated detonation decays to its ultimate planar shape more quickly. Analytical formulae are provided for strongly and weakly overdriven detonations.
Because the fast-reaction limit is employed in the problem formulation, the detonation-oscillation frequency is provided by the initial conditions through the front corrugation wavelength. Therefore, the model cannot be used to predict either the characteristic size of the cusp-like protuberances or the diamond patterns in the burnt-gas flow exhibited by regular detonations. The influence of the inner structure on the shock perturbations and vice versa should be contemplated in a more realistic scenario. In this fashion, the present model could be extended to conditions under which the factor
$\unicode[STIX]{x1D716}$
is purely two-dimensional, that is, when it depends on the shock curvature and the inner detonation structure. Similar to the results presented here, differences in the overall heat released along the corrugated front could be found, which may also lead to constant-amplitude oscillations or purely unstable behaviour. Such a study is left for a future work.
The results presented herein might be applicable to situations in which the type of reaction involved depends on global thermodynamic properties, such as temperature, density and pressure. Regardless of the scenario considered, the limits defining the distinguished regimes should be adapted to include the particular non-adiabatic model because others may not strictly scale with the propagating Mach number squared. Moreover, the effects neglected in this work, such as unsteadiness and acoustic coupling to the supported piston, should also be considered when required by the particular problem. As shown by Bates (Reference Bates2012, Reference Bates2015), the latter may alter the stability limits.
Acknowledgements
The work of C.H. was supported by grants ENE2015-65852-C2-1-R (MINECO/FEDER, UE) and BINV-hBbhOeJQ (Fundación Iberdrola España). C.H. is fully indebted to his PhD advisor and friend, Gustavo Wouchuk. The authors are also deeply grateful to Forman A. Williams for his encouragement and guidance in the preparation of this work and to A. Cuadra for his help in the review process.
Appendix A. Laplace transform
Following the same mathematical treatment as that employed originally by Zaidel (Reference Zaidel1960, Reference Zaidel1961) and in more recent publications by Wouchuk et al. (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009) and Huete et al. (Reference Huete, Jin, Martínez and Luo2017), the hyperbolic transformation,
is employed. The variable
$\unicode[STIX]{x1D712}=\text{const.}$
represents a planar front moving in the burnt gas along the
$x$
-axis, then sweeping the burnt-gas domain from the weak discontinuity
$x=0$
$(\unicode[STIX]{x1D712}=0)$
to the reacting-shock front
$x={\mathcal{M}}_{d}\unicode[STIX]{x1D70F}$
$(\tanh \unicode[STIX]{x1D712}_{d}={\mathcal{M}}_{d})$
. As sketched in figure 13, the initial-value problem formulated in Cartesian coordinates and depicted in figure 3 transforms into a boundary-value problem, with the moving boundary condition, the reacting shock, delimiting a rectangular-shaped domain. The hyperbolic transformation is used to open the vertex defining the origin, thereby modifying characteristic paths.
Integration domain in the hyperbolic variables
$r$
and
$\unicode[STIX]{x1D712}$
. Blue dotted lines refer to the particle paths, grey dashed lines indicate the positive and negative characteristic paths, respectively, and the red solid line defines the detonation trajectory.

The linearized RH equations (2.17) and (2.18) are rewritten as
must be particularized at the detonation front. The sound wave equation (2.13) reduces to
The analysis continues by applying the Laplace transform on the functions
$\hat{p}$
and
${\hat{h}}$
, which renders
respectively. The sound wave equation (A 4) and the function
${\hat{h}}(r,\unicode[STIX]{x1D712})$
can be written in terms of the variables
$\sinh \unicode[STIX]{x1D6FC}=s$
and
$\unicode[STIX]{x1D712}$
to give
After some algebra, the above system of equations can be integrated to provide
Considering a piston in the initial-value problem is found to significantly affect the stability threshold, as pointed out in Bates (Reference Bates2012, Reference Bates2015). Nonetheless, for the sake of simplicity, the isolated-shock boundary condition is assumed in this work. It reduces to
$F_{+}=\text{const.}$
, thereby being defined by the initial condition
$F_{+}=\hat{p}_{d0}/2$
.
With use made of the above relationships, and defining the Laplace function of the detonation shape, namely,
the boundary conditions at the reacting front are recast as
Straight manipulation in (A 9a ) and (A 9b ) gives
and
as the Laplace transforms of the front ripple and pressure perturbation functions, respectively. The factor
stands for the scaled
$\unicode[STIX]{x1D70E}_{a}$
factor.
Separation of variables yields Bessel functions as the family of solutions that satisfies (A 4). The result for the pressure field is then conveniently rewritten in the form
since Bessel functions of the second type must be excluded in order to avoid a divergent behaviour at the origin.
The Laplace transform is also applied to (A 13) at the detonation
$\unicode[STIX]{x1D712}=\unicode[STIX]{x1D712}_{d}$
to give
whose equivalence with (A 11) is employed to provide the factor
$N_{\unicode[STIX]{x1D708}}$
accompanying the Bessel functions, namely,
The even values for
$N_{2\unicode[STIX]{x1D708}}$
have been found to be zero and the series calls for the first values, the initial conditions
to be started. It is readily seen that
$\unicode[STIX]{x1D709}_{d}(s)$
can also be derived in a similar form, where
is related to (A 10) to obtain the factor
$D_{\unicode[STIX]{x1D708}}$
. The recurrence relation remains the same
but with the odd values for
$D_{2\unicode[STIX]{x1D708}+1}$
being all null, and with the initial factors being










