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 Reference Falle1991). 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. (Reference Mi, Timofeev and Higgins2017) 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 Reference Taylor1950) 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 (Reference Wang, McDonald and Radulescu2025). Also shown in the figure is the self-similar solution of Taylor (Reference Taylor1950) 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 Reference Radulescu2020):

(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 (Reference Radulescu2020), 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 (Reference Lau-Chapdelaine and Radulescu2021) 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. (Reference Mi, Timofeev and Higgins2017). 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. Reference Brochet, Manson, Rouze and Struck1963). 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. (Reference Sow, Lau-Chapdelaine and Radulescu2023) 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. (Reference Jackson, Lee and Shepherd2016) 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. Reference Gao, Ng and Lee2015). 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. (Reference Frederick, Gejji, Shepherd and Slabaugh2023). 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 Reference Vasiliev and Nikolaev1978; Crane et al. Reference Crane, Shi, Lipkowicz, Kempf and Wang2021).

We have not attempted a quantitative comparison with the approximate galloping detonation models of Ul’yanitskii (Reference Ul’yanitskii1981) and Mi et al. (Reference Mi, Timofeev and Higgins2017). 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. (Reference Mi, Timofeev and Higgins2017) 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.