A gas-dynamic model for the dynamics of galloping detonations

Abstract

A model for galloping detonations is conceived as a sequence of very fast re-ignitions followed by long periods of evolution with quenched reactions. Numerical simulations of the one-dimensional Euler equations are conducted in this limit. While the mean speed and structure is found in reasonable agreement with Chapman–Jouguet theory, very strong pulsations of the lead shock appear, along with a train of rear-facing N-waves. These dynamics are analysed using characteristics. A closed-form solution for the lead shock dynamics is formulated, which is found in excellent agreement with numerics. The model relies on the presence of a single time scale of the process, the pulsation period, which controls the shock dynamics via the shock change equations and establishes a shock decay with a single time constant. These long periods of shock decay with known dynamics are punctuated by energy release events, with ‘kicks’ in the shocked speed controlled by the pressure increase and resulting lead shock amplification. Model predictions are found in excellent agreement with previous numerical results of pulsating detonations far from the stability limit.

1. Introduction

Detonation waves are supersonic combustion waves sustained by energy release. Their mean speed is approximately given by the Chapman–Jouguet (CJ) model, which assumes the wave is isolated from the rear expansion waves by a sonic surface. In gases, the structure of detonations is cellular. The global reactivity details in each cell are out of phase with those of its neighbouring cell. This gives rise to non-negligible transverse hydrodynamics. The adjacent cells are coupled gas-dynamically via transverse shocks and give rise to quasi-periodic longitudinal and transverse pulsations. The propagation mechanism of cellular detonations relies on the essential role played by transverse wave collisions in re-generating new cells by rapid energy release. Phenomenological models have been proposed relying on discrete explosion centres at the apex of each cell synchronising the dynamics (Vasiliev & Nikolaev 1978; Crane et al. 2021; Monnier et al. 2023), a view we pursue in the present study.

An analogous but simpler problem arises when detonations propagate close to their limits in tubes with a diameter significantly smaller than the characteristic dimension of the cell structure. A mostly one-dimensional galloping phenomenon is observed, whereby transverse dynamics is mostly suppressed during long decay phases, where the lead shock decays as an inert shock (Saint-Cloud et al. 1972; Edwards & Morgan 1977; Ul’yanitskii 1981; Haloua et al. 2000; Gao, Ng & Lee 2015; Jackson, Lee & Shepherd 2016; Huang et al. 2024). These long decays are punctuated by the sudden re-ignition of the gas behind the lead shock and a very rapid reaction of all gas accumulated by the lead shock, as clearly visualised by Edwards & Morgan (1977). The model we formulate in this study exploits this limit.

Sow, Lau-Chapdelaine & Radulescu (2023) observed a similar structure with a distinct long decay and fast re-amplification in numerical simulations of one-dimensional galloping detonations far from the stability limit of steady detonations. While their simulations indicated a pulsation length comparable with the cell length of cellular detonations, the gallop cycle observed in experiments is much longer, owing to the wall losses present in small-scale tubes.

In the present study, we focus on galloping detonations in the absence of wall losses. We wish to determine the gas-dynamic response of a gas to periodic forcing by fast reactions. Recognising that the galloping dynamics is characterised by two very distinct time scales of long shock decays with little exothermicity and very fast re-amplification transients, we study the dynamics in the limit of the re-ignition transient being infinitely fast, i.e. energy addition at constant volume. The model that we hence formulate and study is the inert hydrodynamics punctuated by constant-volume ignition of gases accumulated between the shock and the reactant–product interface since the previous energy release event.

The model we consider bears some resemblance to the phenomenological model of Ul’yanitskii (1981) for galloping detonations and the more recent numerical experiments of Mi, Timofeev & Higgins (2017). In those previous studies, the galloping detonations were treated by localised energy release sources followed by inert hydrodynamics. While the model of Ul’yanitskii (1981) was entirely phenomenological and relied on Taylor–Sedov blast waves, Mi et al. (2017) studied numerically the gas dynamics of shocks supported by concentrated energy sources with infinite density. The model we present is a spatially resolved version of Mi’s model, a feature which permits us to model the dynamics in closed form. Preliminary results reported in this paper were communicated at a conference (Radulescu & Shepherd 2015). Lau-Chapdelaine & Radulescu (2021) subsequently studied a similar model in the context of a toy model based on the reactive form of the Burgers equation. The model we investigate also falls into the broader class of models for cellular detonations with concentrated energy release (Vasiliev & Nikolaev 1978; Crane et al. 2021) for multiple dimensions. These models treat the decaying shocks from energy release in an ad hoc way as Taylor–Sedov blast waves. The results of the present study suggest a different avenue for closure of the cellular detonation problem, which we will address in a sequel.

The paper is organised as follows. Section 2 presents the model formulation and the numerical method we used for its solution. Section 3 presents numerical results. Section 4 formulates a simple closed-form gas-dynamic model for the dynamics; its predictions are compared with experiments and numerical results in § 5.

2. Model formulation

The model we wish to study is the one-dimensional propagation of a reaction wave in a compressible medium in a closed tube with pulsed energy release. Figures 1 and 2 illustrate the model’s formulation. The gas is initially at rest, at state 0. At time 0, a small fraction of gas is instantly reacted near the wall, which drives a lead shock. At subsequent increments of time $\tau$ , all the non-reacted gas between the leading shock and the contact surface separating the burned gases (from the previous cycle) and the newly shocked gas is reacted. This process continues until the wave front reaches a quasi-steady state of periodic pulsations, which we wish to analyse and model. The time delay $\tau$ is the model’s input.

Figure 1. The problem studied is the response of a gas to an instantaneous (very rapid) energy deposition followed by inert evolution; here the reaction in the layer of unburned gas is followed by a subsequent inert evolution, during which non-reacted gas accumulates between the lead shock and the contact surface separating reactants from products. The energy is released at time intervals $\tau$ .

Figure 2. Initial evolution of the pressure field following the release of energy at discrete times 0, 1, 2, 3, etc. between the lead shock and the contact surface separating the reacted gas from a previous cycle from non-reacted gas accumulated behind the lead shock, with $\gamma =1.2$ and $Q=50$ . See the supplementary material for animation.

We assume the gas to be a perfect gas with constant specific heats, characterised by the ratio of specific heats $\gamma$ . The energy released $Q$ per unit mass of reactant is modelled through the simple equation of state for the specific internal energy of the gas: $e(p,\rho ) = p \rho ^{-1} (\gamma -1)^{-1} +\lambda Q$ , where $p$ is the pressure, $\rho$ the gas density and $Q$ the heat release. The mass fraction of reactant $\lambda$ serves to describe the reaction progress, which takes a value of 1 in the reactants and 0 in the products.

The dynamic process we wish to study is the outcome of periodic release of energy over a very short time scale in the gas accumulated behind a leading shock. Between each ‘kick’, the gas evolves as an inert gas, reactants and product gases obeying the inert gas-dynamic equations.

Both the inert and reactive phases are described by the Euler equations, written in conservation form as

(2.1) \begin{align} \left ( \rho \right )_t + \left ( \rho u \right )_x=0; \quad \left ( \rho u \right )_t + \big ( \rho u^2 +p \big )_x = 0; \end{align}

(2.2) \begin{align} \left ( \rho \left (e+\frac {1}{2}u^2\right ) \right )_t + \left ( \rho u \left (e+\frac {1}{2}u^2\right )+p u\right )_x = 0 , \\[9pt] \nonumber \end{align}

where the evolution of the reactant mass fraction is modelled by a simple linear depletion rate:

(2.3) \begin{align} \left ( \rho \lambda \right )_t + \left ( \rho u \lambda \right )_x=- \rho K \lambda \end{align}

with the coefficient $K$ controlling the rate of the reaction. During the inert phase, the rate is turned off, while during the short reaction phase, the evolution of $\lambda$ is calculated explicitly in the gases behind the shock. The value of $K$ is chosen such that the local evolution of $\lambda$ is much faster than the inert cycle, which is controlled by the input forcing period $\tau$ .

To help in the interpretation of the numerical results, the evolution of the C+ characteristics was also calculated by evolving a label $\alpha$ for each C+ characteristic curve. The requirement that these labels remain constants along paths ${\rm d}x/{\rm d}t=u+c$ can be written as $\alpha _t + ( u+c )\alpha _x=0$ .

In all subsequent exposition, we normalise all variables of our problem by the characteristic scales $p_0$ , $\rho _0$ and $\tau$ . Non-dimensional speed becomes, for example, $\overline {u}=u/\sqrt {p_0/\rho _0}$ . Dropping the over-line from all non-dimensional terms, we recover the same system of equations and parameters, re-scaled as explained. The forcing period becomes $\tau =1$ and the initial state of the gas is $p_0=1$ and $\rho _0=1$ . For the calculations presented below, we take $\gamma =1.2$ and $Q=50$ . This choice of parameters represents approximately most gaseous detonations of interest and is the standard parameter set used in the literature.

Figure 3. The evolution of the cycle-averaged shock speed. The black line is the CJ speed; the red dashed line is the steady speed observed.

The governing equations were solved using the code MG developed by Falle at the University of Leeds (Falle 1991). The code uses Strang splitting of the reactive and hydrodynamics substeps. The calculations are readily performed on a personal laptop on a single thread. For reference, the detonation travelled approximately 300 grid points in one cycle of oscillation. A resolution study permitted the establishment of the convergence of our results. In the calculations reported below, we set $K=100$ , such that the energy deposition is 100 times faster than the rate at which the kick is repeated (formally integrating (2.3) at constant volume yields $\lambda =\exp (-K t)$ ). The reactive portion of the calculation was performed by reducing the time step set by the Courant–Friedrichs–Lewy condition based on the stable time step set by the hydrodynamics by a factor of 100, in order to ensure the proper resolution of the very fast transient associated with the energy release at nearly constant volume. The labels $\alpha$ were advected non-conservatively using an upwind approximation.

3. Numerical results

The evolution of the pressure field is shown in a space–time diagram in figure 2. Each energy release event generates a nearly instantaneous increase in pressure given by the constant-volume combustion of $\Delta p = (\gamma -1)\rho Q$ . This sudden pressure increase does not modify the density and speed of the gas inside the gas noticeably. This results in a rear-facing shock and sudden amplification of the lead shock, and two centred expansion waves propagating inward the high-pressure region, aimed to relieve the pressure. The lead shock subsequently decays, while a train of rear-facing waves is established. These reflect from the end wall, but do not catch up to the lead shock. An animation of this transient pressure evolution is provided as supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2025.10547.

The mean shock speed for every inert cycle was evaluated as the ratio of the length travelled by the lead shock over the time duration of one cycle, which is unity. The evolution of the cycle-averaged speed is shown in figure 3 for until approximately 40 kicks. A steady-state propagation appears to be established after approximately 7 kicks, reaching a speed in excess of the CJ value by 3 %. Mi et al. (2017) also reported speeds in excess of the CJ value for their randomised discrete system of the same magnitude. We attribute the excess speed to the tampering effect of the rear-facing shocks driven at each kick. The otherwise isentropic Taylor expansion in the idealised steady detonation (Taylor 1950) is now replaced by the decay of a rear-facing train of shocks, which form a train of N-waves. An extended discussion of the super-CJ propagation can be found in the Appendix.

Figure 4. The pressure distribution at $t=44$ . The green dashed line is the Taylor self-similar solution.

The pressure distribution at $t=44$ is shown in figure 4. The train of decaying rear-facing N-waves is clearly evident, in good accord with the model suggested by Wang, McDonald & Radulescu (2025). Also shown in the figure is the self-similar solution of Taylor (1950) for the pressure distribution behind a hypothetically infinitely thin detonation propagating away from a closed wall at the constant CJ speed. The pressure evolution is found in general good agreement, suggesting that the ideal model provides a correct representation of the phase-averaged dynamics.

Figure 5. The pressure distribution in an inertial frame moving with the mean lead shock speed. The red line denotes the interface between reacted and non-reacted gas while thin grey lines are forward-facing characteristics; the thick grey line is the limiting characteristic, while the thick dotted line is the lead signal coming from the sudden de-pressurisation of the high-pressure region. See the supplementary material for animation.

Figure 5 shows the pressure distribution in an inertial frame moving with the mean lead shock speed. An animation of the pressure evolution is provided as supplementary material. The rear-facing shock and expansion wave associated with the reactive event and its de-pressurisation are evident. Of interest is the head of the forward-facing expansion wave resulting from the de-pressurisation of high-pressure reacted gas towards the product side. This arrives at the lead shock only near the end of the cycle. We thus see that the lead shock decay takes its information from the entire high-pressure region only. Below, we develop a simple model exploiting this fact.

Figure 6. Lead shock speed decay in one cycle. Lines are exponential decays. The red line is closed-form model with a mean speed $\overline {D}=D_{CJ}$ ; the green line is for $\overline {D}=1.03 D_{CJ}$ .

Figure 6 shows the speed evolution in a single cycle. The decay is very well approximated by an exponential decay, i.e. by a single characteristic time. That a single time emerges from the dynamics is not unexpected, as the only time scale of the problem is $\tau$ in the limit of fast energy deposition at intervals $\tau$ . We have verified that the exponential time constant changes with $\tau$ in the same proportion, but this is inherent in our problem scaling.

4. Local decay model

The question that arises is what sets the decay rate and amplitude of the pulsations. The characteristic analysis above showed that the lead shock takes its information in the high-pressure region alone. The energy deposition is at constant volume, which also implies that the gas speed does not change during the reactive event. This signifies that the speed gradient remains constant across the jump. Through the shock change equations relating derivatives of the field variables with the dynamics of the shock, we have for planar strong shocks (Radulescu 2020):

(4.1) \begin{align} \frac {\partial u}{\partial x}=-\frac {6}{\gamma +1}\frac {\dot {D}}{D}. \end{align}

This means that the velocity gradient near the end of the cycle evaluated at the shock is equal to the gradient at the start of the cycle behind the shock. In the absence of another time scale, a likely situation is that the gradient remains constant during the entire cycle. This is in accord with our simulations. This provides the explanation of why $\dot {D}/D$ remains constant and the decay is exponential in one cycle. The time constant is found in terms of $\tau$ after the magnitude of each pulsation is determined, as developed below.

Note that the pulsed problem we have studied only has the forcing time scale as characteristic time as long as all shocks are strong. As can be shown using the shock change equations of Radulescu (2020), finite-strength shocks have a multitude of time scales involving the departure from the strong shock approximation. While the lead shock is strong in the present case, with a Mach number of 6.2, rear-facing shocks driven by each kick are not; this may explain the small departures from an exact exponential decay apparent in figure 6. An alternative explanation was provided by Lau-Chapdelaine & Radulescu (2021) for the Burgers equation involving the time histories of forward-facing characteristics arriving at the shock in the early or later phases of the dynamics.

The time constant of decay and the amplitude of the jump in speed are obtained from the Riemann problem of a new shock driven by a pressure increase from reactions, and the kinematic constraint between the speeds at the start and end of the cycle with the mean speed (controlled by energetics).

Figure 7. The Riemann problem for determining the new shock speed $D_0$ after the constant-volume energy deposition changing the post-shock state from 1 to 2.

Consider the diagram of figure 7. The shock speed prior to the energy release is $D_1$ . The post-shock state 1 is obtained from the Rankine–Hugoniot relations. We show the approach first for a strong shock for simplicity, but provide the results for a finite-strength shock at the end of this section. For a strong shock, the jump equations are

(4.2) \begin{align} p_1=\frac {2}{\gamma +1}\rho _0 D_1^2,\quad \rho _1= \frac {\gamma +1}{\gamma -1}\rho _0,\quad u_1 = \frac {2}{\gamma +1} D_1. \end{align}

Upon energy release at constant volume, state 2 becomes

(4.3) \begin{align} p_2 = p_1 + \rho _1 \left (\gamma -1 \right ) Q, \quad \rho _2= \rho _1, \quad u_2 = u_1, \quad c_2 = \sqrt {\gamma \frac {p_2}{\rho _2}}. \end{align}

State 4 is obtained from the invariance of the Riemann invariant $J_+=2c/(\gamma -1) + u$ on a C+ characteristic connecting states 2 and 4. Since the expansion wave is isentropic, we get

(4.4) \begin{align} p_4 = \left (\frac {c_4}{c_2}\right )^{2 \gamma/(\gamma -1)} p_2 = \left ( 1+ \frac {\gamma -1}{2} \left ( \frac {u_2-u_4}{c_2}\right ) \right )^{2 \gamma /(\gamma -1)} p_2 .\end{align}

While the pressure and flow speed across the new shock propagating at speed $D_0$ are

(4.5) \begin{align} p_5=\frac {2}{\gamma +1}\rho _0 D_0^2, \quad u_5 = \frac {2}{\gamma +1} D_0. \end{align}

The slip line separating states 4 and 5 is massless; hence, mechanical equilibrium requires $p_4=p_5$ and $u_4=u_5$ . Using the simple expressions developed, one can thus rewrite (4.4) as an implicit relation $D_0(D_1, \gamma , Q)$ , which takes the form

(4.6) \begin{align} 2 D_0^2-\beta \left (1+ (D_1 -D_0) \sqrt { \frac {\gamma -1}{\gamma \beta }} \right )^{2 \gamma /(\gamma -1)}=0 \end{align}

with $\beta =2D_1^2 +(\gamma +1)^2 Q$ . If we are dealing with a strong detonation, for which the CJ propagation speed is given by $D_{CJ}^2=2 (\gamma ^2-1 )Q$ , then $\beta$ can be rewritten in terms of $D_{CJ}$ :

(4.7) \begin{align} \beta =2D_1^2 + \frac {1}{2} \frac {\gamma +1}{\gamma -1}D_{CJ}^2 .\end{align}

The second constraint is the mean detonation speed being given by the CJ solution. For a constant shock decay rate leading to an exponential decay in time, the mean detonation speed is given by

(4.8) \begin{align} \overline {D}\equiv \ \frac {1}{\tau } \int _0^\tau D(t) \mathrm{d}t= \frac {D_0 - D_1}{\ln (D_0/D_1)}. \end{align}

Requiring that $\overline {D}\simeq D_{CJ}$ yields a closed system of equations for the shock speeds $D_0$ and $D_1$ compatible with the mean speed being given by the CJ solution and the jump being given by the energy release. For $\gamma =1.2$ and $Q=50$ , we obtain $D_0 = 8.12$ and $D_1= 5.34$ , while the CJ detonation speed assuming a strong detonation is $D_{CJ}=6.63$ . Good agreement is obtained, as indicated in figure 6.

A more accurate treatment is to relax the strong shock assumption to compute the solution. The treatment is identical but the expressions are more involved. For $\gamma =1.2$ and $Q=50$ , we obtain $D_0 = 8.09$ and $D_1= 5.67$ , while the exact CJ detonation speed is $D_{CJ}=6.81$ . Very good agreement can be seen for the speed distribution in a single cycle, as shown in figure 6. If we account for the observed speed overdrive of approximately 3 % in the pulsed problem, the agreement improves, as shown in figure 6.

Once the ratio $D_0/D_1$ is determined, the relation between the kick interval $\tau$ and shock decay rate $\zeta = - \dot {D}/D$ follows. If $\zeta$ is constant, as implied by the shock change equation argument, the shock decay is given by

(4.9) \begin{align} D=D_0 \exp \left (-\zeta t \right ). \end{align}

Since $D(\tau )=D_1$ at the end of cycle, we obtain

(4.10) \begin{align} \zeta =\frac {1}{\tau }\ln \left (\frac {D_0}{D_1} \right ). \end{align}

For the model parameters considered, we obtain $\zeta = 0.35/\tau$ . This result relies on the condition that the shock decay is only controlled by the time scale $\tau$ and other time scales are not relevant, i.e. chemical times are much shorter and resonant times are fixed by $\tau$ alone. This is obviously an approximation, which appears to capture the main dynamics correctly.

5. Model validation

Our numerical simulations suggest the mean speed of detonations with pulsed reactions to be slightly in excess of the CJ value, consistent with the finding of Mi et al. (2017). We suggest this effect is due to the tampering provided by the rear-facing shocks. Interestingly, this is also consistent with experiment, in that detonation speeds in the limit of large tubes and zero losses are usually in excess of the CJ value by the same proportion (Brochet et al. 1963). Further discussion is presented in the Appendix.

It is worthwhile to compare the pulsation amplitude obtained in the present study with available data in the literature for lossless conditions. We are only aware of the numerical study of Sow et al. (2023) far from the stability boundary. For the heat release parameter they used of $Q=9.1$ , our model predicts amplitudes of pulsations of $\pm$ 10 % of the CJ value, which is in excellent agreement with their results.

Jackson et al. (2016) studied the dynamics of galloping detonations experimentally and observed velocity excursions around the mean of approximately $\pm$ 30 % of the mean propagation speed. This is larger than predicted by the present model of $\pm$ 20 %, the difference attributable to the finite-time wave amplification occurring in the real process, which generates transient forward-facing shocks in the shock non-reacted gas. When only the low-speed excursions are considered, these are found to range between 20 % and 30 % below the mean propagation speed (Gao et al. 2015). The agreement with the model is satisfactory.

The closed-form solution we have obtained for the shock dynamics and jumps due to the sudden energy release is also compatible with the experiments of cellular detonations performed by Frederick et al. (2023). While we predict velocity excursions of the order of 20 % from the mean propagation speed, the experiments suggest departures of the order of approximately 30 %. The widening effect played by the irregularity of the front notwithstanding, the geometrical focusing of the shocks in the cellular problem is likely contributing to this small difference. Nevertheless, the general good agreement of our model with experiment suggests that the magnitude of velocity pulsations in cellular detonations is mainly due to the energy release from gas igniting from transverse wave collisions at the end of the cell, in good agreement with existing phenomenology (Vasiliev & Nikolaev 1978; Crane et al. 2021).

We have not attempted a quantitative comparison with the approximate galloping detonation models of Ul’yanitskii (1981) and Mi et al. (2017). These models are phenomenological and posit extreme assumptions on the energy release, which are incompatible with the physics of the model we develop, and with experiment. Their models assume point-source Taylor–Sedov blast waves and predict non-realistic infinitely strong shocks at the start of the cycle. There are several questionable modelling assumptions. Ul’yanitski tacitly assumes that the blast energy is equally partioned between forward and backward shock motion. Mi et al. (2017) recognised the importance of this partition and adopted Sakurai’s model to account for the density differences in reactants and product gases. They treated both blasts as strong, neglecting the kinetic and internal energy of the gas in front of the shocks. We believe this strongly violates the physics in our problem, where the rear-facing shock is found to be very weak (its Mach number is approximately 1.5 and overpressure of 0.5), since it results from a constant-volume energy addition in a finite kernel of shocked gas.

In conclusion, we have proposed a closed-form model for the pulsation dynamics in detonations with periodic release of energy in a spatially resolved zone. The dynamics of the lead shock is obtained with only the forcing period $\tau$ as parameter and the assumption that the time-averaged propagation is well approximated by the CJ solution if the wave were steady. Good to excellent agreement is obtained with experiment and numerics for the lead shock pulsations. We note that the results were obtained for lossless systems. To adapt them to realistic galloping phenomena, losses to the tube walls need to be incorporated. This is left for future study.