3.1. Hydrodynamics

In the operator splitting, or Strang splitting, approach that we use, the contribution of hydrodynamics while keeping a frozen composition is to solve the governing equations (2.2)–(2.5) without the chemical species evolution term. This yields

(3.2) \begin{align} \left (\frac {\partial v}{\partial t}\right )_\phi -\left (\frac {\partial u}{\partial \phi }\right )_t=0 ,\quad\end{align}

(3.3) \begin{align} \left ( \frac {\partial u}{\partial t} \right )_\phi + \left (\frac {\partial p}{\partial \phi }\right )_t=0,\quad \end{align}

(3.4) \begin{align} \left ( \frac {\partial e_{tot}}{\partial t} \right )_\phi + \left ( \frac {\partial pu }{\partial \phi } \right )_t =0, \end{align}

(3.5) \begin{align} \left ( \frac {\partial Y_i}{\partial t} \right )_\phi = 0.\end{align}

The hydrodynamic step uses a second-order HLLE scheme devised for a structured uniform grid in the $\phi$ dimension, as shown schematically in figure 3. The solution vector $U=(v, u, e_{tot}, Y_1, \ldots , Y_N)$ is stored at the cell centres and the ‘fluxes’ $F=(-u, p,pu, 0, \ldots , 0)$ are evaluated at the cell interfaces. The solution inside one of these cells is represented by piecewise limited linear functions. Central differences are used to reconstruct the slopes and the Van Albada limiter (Van Albada et al. Reference Van Albada, Van Leer and Roberts1997) has been applied to limit variations in conserved variables. The cell-averaged values of the conserved quantities $U_i$ are updated by computing the flux at the two numerical cell interfaces $\mathcal {F}_{i-1/2}$ and $\mathcal {F}_{i+1/2}$ , where $i$ is the index of cells. The inter-cell flux $\mathcal {F}$ is evaluated through the approximation solver of the Riemann problem using the HLLE flux functions with the reconstructed left and right state at the cell boundary described by Einfeldt et al. (Reference Einfeldt, Munz, Roe and Sjögreen1991). The left and right states are interpolated values at the cell on the left- and right-hand side of the cell interface. The second-order spatial–temporal step is performed as follows. First, a first-order approximation is obtained:

(3.6) \begin{equation} \hat {\bar {U}}_{i}^{n+1}=\bar {U}_{i}^{n}+\dfrac {\Delta t}{\Delta x}\left [\mathcal {F}_{i+1/2}^{n}-\mathcal {F}_{i-1/2}^{n}\right ]. \end{equation}

This solution $\hat {\bar {U}}_{i}^{n+1}$ together with the initial vector $\bar {U}_{i}^n$ are then used to update the second-order cell state $\widetilde {U}_{i}^{n+1}$ according to

(3.7) \begin{equation} \widetilde {U}_{i}^{n+1}=\bar {U}_{i}^{n}+\dfrac {\Delta t}{2\Delta x}\left [\mathcal {F}_{i+1/2}^{n}+\hat {\mathcal {F}}_{i+1/2}^{n}-\mathcal {F}_{i-1/2}^{n}-\hat {\mathcal {F}}_{i-1/2}^{n}\right ]. \end{equation}

When the reactions are coupled, $\widetilde {\bar {U}}_{i}^{n+1}$ is further updated after the reaction step. During the hydrodynamic substep, the ideal gas law is used with an isentropic coefficient, $\gamma$ , provided by Cantera at the end of the last chemical substep. The local value of $\gamma$ is assumed constant during a hydrodynamic step.

The boundary $\phi =0$ corresponds to the piston face, where the gas speed is that of the piston. In order to impose a flux into the first cell, the pressure is determined by solving the exact Riemann problem with the velocity of the piston at the cell boundary. Once the velocity and the pressure are evaluated in the cell boundary, the flux can be prescribed. The boundary condition at the right boundary of the computational domain uses an extrapolation of cubic interpolation; nevertheless, the calculation ends before the lead shock reaches the right boundary such that this local boundary condition does not have any role in the calculation.

Figure 3.

At a given time, the solution in neighbouring cell mass elements is represented by piecewise linear functions used to achieve second-order accuracy in the solution of the Riemann problem at cell interfaces for calculating inter-cell fluxes.

3.2. Chemical species evolution

In the operator splitting, or Strang splitting, approach that we use, the contribution of the change of chemical composition of the gas is to solve the governing equations (2.2)–(2.5) without the flux terms (derivatives with $\phi$ ). This yields

(3.8) \begin{align} \left (\frac {\partial v}{\partial t}\right )_\phi =0 ,\quad\end{align}

(3.9) \begin{align} \left ( \frac {\partial u}{\partial t} \right )_\phi =0,\quad \end{align}

(3.10) \begin{align} \left ( \frac {\partial e_{tot}}{\partial t} \right )_\phi =0,\,\, \end{align}

(3.11) \begin{align} \left ( \frac {\partial Y_i}{\partial t} \right )_\phi = v \omega _i . \\[6pt] \nonumber \end{align}

During the reactive step, the specific volume, speed and total energy density (for a fixed mass) are to be kept constant. This substep thus naturally lends itself to using Cantera’s computational framework already available for this canonical calculation at constant volume and energy. The composition and temperature of each cell are thus evolved as a constant-volume and energy reactor in Cantera during the chemical reaction step. The ordinary differential equations are solved using the SUNDIALS stiff ODE solver used by Cantera over the global time step $\Delta t$ dictated by the hydrodynamics.

In the Strang splitting approach used, the sequentially taken hydrodynamic and reactive substeps thus provide the entire solution vector $(v, u, e_{tot}, Y_1, \ldots , Y_N)$ at the end of a global time step. We use Cantera’s built-in solvers to determine the pressure and temperature corresponding to this uniquely defined thermodynamic state by its specific volume, energy and composition.

3.3. Numerical verification

3.3.1. Verification of the hydrodynamic solver

Two frozen hydrodynamic problems have been chosen to verify the correct implementation of the hydrodynamic solver and confirm the adequacy of the numerical method proposed. The first problem is a Riemann problem involving a weak shock and expansion fan in a perfect gas. An initial discontinuity at $x=1$ separates the left state $ (u_l=0, p_l=2, \rho _l=1 )$ from the right state $ (u_r=0, p_r=1, \rho _r=0.5 )$ with the isentropic coefficient $\gamma =1.4$ constant for both sides. Figure 4 shows the comparison of the calculated density with the exact solution at time $t=0.3$ . It is found that the numerical scheme treats satisfactorily the shock, contact surface and expansion wave typical of the HLLE scheme.

The second test problem is a strong shock case in a real gas, with frozen chemical composition. A piston with a constant speed of 1500 ${\rm m\,s}^{-1}$ moves into a stoichiometric mixture of hydrogen and oxygen at pressure $p_r=101\,325$ Pa and temperature $T_r=293$ K. Figure 4 shows the comparison of the calculated density with the exact solution at time $t=0.03$ s, where excellent agreement is observed between the computed and the exact solution. This confirms the reliability of the numerical solver.

Figure 4.

Density profile for shock-expansion problem (left); strong shock problem (right).

3.3.2. Verification of the reactive solver coupling to hydrodynamics

The coupling between the reactive and hydrodynamic solver was tested in stable detonation wave propagation. The test is whether the solver can propagate a stable Zel’dovich–von Neumann–Döring (ZND) detonation wave. We tested an over-driven ZND detonation in a stoichiometric mixture of hydrogen and oxygen initially at $1.01\times 10^5$ Pa and 300 K. A piston speed of 2583.6 m s⁻¹ drives an over-driven detonation wave propagating at 3500 m s⁻¹. For reference, the CJ self-sustained detonation wave speed in this gas is 2839 m s⁻¹ and the material speed at the CJ state is 1293 m s⁻¹. The driving piston being larger than the CJ material speed gives as expected an over-driven detonation propagating at a speed exceeding the CJ value.

The calculation was initialised with the over-driven ZND detonation profile calculated separately using the Shock and Detonation Toolbox (EDL 2023) in Cantera. Figure 5 shows the evolved wave structure at $t=4.8\times 10^{-7}$ s. The calculation used a Courant–Friedrichs–Lewy number of 0.7. Excellent agreement is found between the calculated detonation structure and the expected travelling-wave solution provided by the ZND solution, verifying the numerical coupling between the reactive and hydrodynamic solvers.

Figure 5.

Temperature profile for the ZND test.

4.1. Overview

The ignition of the gas induced by the lead shock in the presence of gas-dynamic fluctuations is controlled by the temperature and pressure variation in the induction zone, prior to the onset of exothermicity. It is thus useful to establish the fluctuations induced by the piston first, without accounting for exothermicity.

A numerical example of the temperature evolution in the induction zone is illustrated in the $\phi {-}t$ diagram of figure 6 by freezing the exothermicity. The mixture is $2{\rm H}_2+{\rm O}_2$ initially at $T_0=300$ K and $p_0=5.90 \times 10^{3}$ Pa, with a mean piston speed $u_{p0}=1.63\times 10^{3}$ m s⁻¹, fluctuation amplitude $A=0.2u_0$ and frequency $f=45.4$ kHz. This leads to a non-fluctuated post-shock temperature $T_1=1100\,{\rm K}$ and post-shock pressure $p_1=1$ atm. The reactive solution to this problem is discussed later.

Figure 6.

Temperature evolution in $\phi {-}t$ space induced by the impulsive piston motion ( $u_0=1626.35$ m s⁻¹, $A=0.2u_{p0}$ and $f=45.4$ kHz) into chemically frozen $2{\rm H}_2+{\rm O}_2$ initially at $T_0=300$ K and $p_0=5900$ Pa; dark dashed blue lines are C− characteristics while cyan dashed blue lines are C+ characteristics.

Figure 7.

Temperature profiles at several times; same conditions as figure 6.

While the impulsive piston motion generates instantly a lead shock, the subsequent compressive parts of the piston motion generate a waveform that steepens to form internal shocks. These internal shocks can be readily identified as right-facing fronts. This forward-facing waveform then reflects on the lead shock, generating reflected disturbances along the C− characteristics and along (vertical) particle paths. One of these evident entropy waves can be identified as originating at the lead shock when an internal shock wave overtook the lead shock at approximately $t=3.8\times 10^{-5}$ s. The particle path evolution can thus be seen to be modulated by the initial temperature obtained at the shock and the subsequent isentropic expansion or compression in regions sufficiently close to the piston before inner shock formation, or further wave by inner-shock heating.

The magnitude and wave shape of the fluctuations can also be observed in figure 7, which shows the temperature distributions at four different times. The last profile corresponds to an instant just prior of the lead shock arrival at the right boundary of the computational domain in this case. Inner shocks and contact surfaces discussed above can be clearly identified. Note, however, that it would have been very difficult to rationalise these features without the space–time diagram illustrating the wave dynamics.

Given the propensity to form inner shocks and their role in shock heating the gas in the induction zone, it is desirable to construct an approximate analytical model for the process to predict the timing of the various events and the temperature amplitude. The problem can be conceptualised as follows. The state generated by an impulsively started steadily moving piston (state 1) behind a constant-speed shock can be taken as the leading-order constant solution. Perturbations to this state are brought about by the harmonic piston motion. Denoting the over-pressure of the lead shock as $z_1=(p_1-p_0)/p_0$ , the Riemann variable $J^-$ across the shock varies as $O(z_1^3)$ while pressure, density and speed vary as $O(z_1)$ (the acoustic solution) (Chandrasekhar Reference Chandrasekhar1943; Whitham Reference Whitham1974). A weakly nonlinear description of the state evolution behind the lead shock can thus be sought assuming $J^-=\text{const.}$ to $O(z_1^2)$ . This constitutes a ‘simple wave’ solution first suggested by Chandrasekhar (Reference Chandrasekhar1943) and further exploited by Whitham (Reference Whitham1974) in different shock-formation problems. This is the first inert problem to which we seek a simple solution, labelled problem A. Its weak nonlinearity permits one to account for the dynamics of inner shock waves.

The problem B we wish to solve is the interaction of right-facing waves described in problem A with the lead shock, which leads to reflected waves along particle paths and C− characteristics. This reflection problem is solved in the linear regime and serves to correct the lead shock strength and the interior states provided by problem A.

We thus seek solutions of the form

(4.1) \begin{eqnarray} u=u_1+u_A+u_B\cdots , \\[-24pt] \nonumber \end{eqnarray}

(4.2) \begin{eqnarray} p=p_1+p_A+p_B\cdots ,\\[-24pt] \nonumber \end{eqnarray}

(4.3) \begin{eqnarray} T=T_1+T_A+T_B\cdots .\,\,\,\\[6pt] \nonumber \end{eqnarray}

4.2. Problem 0: the constant-speed shock driven by the steady piston

The leading-order solution is the shock driven by the impulsively started piston at steady speed $u_{p0}$ into quiescent gas at state 0. This generates a constant state between the piston and the lead shock labelled with subscript 1, which obeys the usual Rankine–Hugoniot shock jump conditions.

4.3. Problem A: right-facing wave train and inner-shock formation

The piston speed departure from the steady state $u_A(\phi =0)=u_p-u_{p0}=A\sin (2\pi f t)$ generates a train of waves propagating into the gas at state 1. The weakly nonlinear waveform generated by an oscillating piston behind the lead shock can be obtained by neglecting $O(z^3)$ variations in entropy and $J^-$ Riemann variable across the lead shock and all internal shocks. This simple wave problem is solved in closed form in Lagrangian coordinates for a perfect gas. We are not aware of its solution elsewhere, although similar arguments were presented by Chandrasekhar (Reference Chandrasekhar1943) and Whitham (Reference Whitham1974) in other problems in Eulerian coordinates, which makes the analysis more complicated.

As detailed in Appendix A, the characteristic equations for isentropic flow of a perfect gas are

(4.4) \begin{equation} \left (\frac {\partial J^\pm }{\partial t}\right )_\phi \pm \rho c \left ( \frac {\partial J^\pm }{\partial \phi }\right )_t=0 ,\end{equation}

where the Riemann variables are $J^\pm =2c/(\gamma -1)\pm u$ . Since all C− characteristics originate from the upstream state 0, $J^-=J^{-}_1=J^{-}_0$ is constant everywhere to the level of the current approximation and we have immediately the simple wave relation between flow and sound speed perturbations:

(4.5) \begin{equation} c=c_1+(u-u_1)\frac {\gamma -1}{2} .\end{equation}

Given the flow is isentropic in the current approximation, all variables can be related to the sound speed variation:

(4.6) \begin{equation} \frac {\rho }{\rho _1}=\left (\frac {p}{p_1}\right )^{1/\gamma }=\left (\frac {T}{T_1}\right )^{1/(\gamma -1)}=\left (\frac {c}{c_1}\right )^{2/(\gamma -1)} \end{equation}

and using (4.5), the propagation speed of forward-facing characteristics in $\phi {-}t$ space can be expressed in terms of $u$ only:

(4.7) \begin{equation} \frac {{\rm d} \phi }{{\rm d}t}=\rho c=\rho _1 c_1 \left (1+\frac {u-u_1}{c_1}\frac {\gamma -1}{2}\right )^{(\gamma +1)/(\gamma -1)}. \end{equation}

Similarly, using (4.5), the Riemann variable $J^+$ can also be expressed in terms of $u$ only. Given $J^+$ remains constant along C+ characteristics, it implies that $u$ remains constant along C+ characteristics (and all other variables by virtue of (4.5) and (4.6)) and C+ characteristics are straight lines. These characteristics hence communicate the values of $u$ , $c$ , $\rho$ , $\rho c$ , $T$ , $p$ , etc., from the face-forward piston. In short, all variables hence satisfy the simple advection equation for $\alpha =u, c, \rho , \rho c,T,p\ldots$ :

(4.8) \begin{equation} \left ( \frac {\partial \alpha }{\partial t} \right )_\phi + \rho _1 c_1 \left (1+\frac {u-u_1}{c_1}\frac {\gamma -1}{2}\right )^{(\gamma +1)/(\gamma -1)} \left ( \frac {\partial \alpha }{\partial \phi } \right )_t=0 .\end{equation}

The trajectory of any C+ characteristic can be obtained by integrating (4.7) from the piston $\phi =0$ and reference time $t^*$ , yielding

(4.9) \begin{equation} \phi =\rho _1 c_1 \left (1+\frac {u(t^*)-u_1}{c_1}\frac {\gamma -1}{2}\right )^{\frac {\gamma +1}{\gamma -1}} (t-t^*) \end{equation}

Since $u(t^*)$ is the piston speed, this expression provides implicitly the dependence $t^*(\phi , t)$ . Given $t^*$ , $u(t^*)$ is known and remains constant along that characteristic. The other variables are obtained from (4.5) and isentropic relations (4.6). Using the harmonic piston speed perturbation, (4.9) can be written as

(4.10) \begin{equation} \phi =\rho _1 c_1 \left (1+\frac {A\sin (2\pi f t^*)}{c_1}\frac {\gamma -1}{2}\right )^{(\gamma +1)/(\gamma -1)} (t-t^*), \end{equation}

which can be further simplified for small-Mach-number harmonic motion or in the Newtonian limit $\gamma \rightarrow 1$ to

(4.11) \begin{equation} \phi =\rho _1 c_1 \left (1+ \frac {\gamma +1}{2}\frac {A\sin (2\pi f t^*)}{c_1}\right ) (t-t^*). \end{equation}

The solution is complete but breaks down when characteristics intersect and the solution becomes multi-valued. When this happens, shocks need to be fitted. The shock formation time and location, as well as the fitted inner-shock trajectories can be obtained in closed form.

Figure 8.

Sketch illustrating the shock formation process along the characteristic passing through the point where $\rho c$ had the steepest rate of increase on the piston face (point X); shock forms at point ^* from characteristics merging.

Figure 8 illustrates the shock formation and the fitted shock trajectory, which approximately bisects the angle formed by the characteristics originating from either side of the shock. Shocks will form along the characteristics where the front was initially steepest at the piston face, i.e. ${\rm d}(\rho c)/{\rm d}t$ was initially maximum, since these steepest points on the waveform will be maintained on the same characteristic. At the piston face, we have

(4.12) \begin{equation} \rho c = \rho _1 c_1 \left (1+ \frac {\gamma +1}{2}\frac {A\sin (2\pi f t^*)}{c_1}\right ) \end{equation}

and the steepest points are located at $t=t_{X}=n/f$ with $n=0,1,2, 3\ldots .$ . Shocks will thus form along these characteristics, which are given by

(4.13) \begin{equation} (\rho c)_X = \rho _1 c_1 = \frac {\phi }{t-t_X}. \end{equation}

These characteristics will intersect with a neighbouring characteristic originating at point $\phi =0$ and $t=t_X+\delta t$ and having a propagation speed $(\rho c)_X+ ({{\rm d} \rho c}/{{\rm d}t}) \delta t$ at the shock formation time

(4.14) \begin{equation} t^*=\frac {n}{f} + \frac {1}{(\gamma +1)\pi f A/c_1} \end{equation}

and shock formation particle label

(4.15) \begin{equation} \phi ^*=\frac {\rho _1 c_1}{(\gamma +1)\pi f A/c_1} \end{equation}

in the limit of $\delta _t\rightarrow 0$ . For the example considered, the shock is predicted to form at particle label $\phi ^*=1.26\times 10^{-3}$ kg m $^{-2}$ , in very good agreement with the simulations shown in figure 9.

Figure 9.

Temperature evolution in chemically frozen $2{\rm H}_2$ + ${\rm O}_2$ initially at $T_1=1100$ K generated by a piston speed with fluctuation amplitude $A=325.27$ m s⁻¹; numerical result (top) and analytical prediction (bottom).

The fitted shock trajectory can be approximated to be the continuation of the characteristic along which the shock has first formed. This approximation is the leading-order solution for weak shocks (Whitham Reference Whitham1974), since the speed of weak shocks is the average of the speed of the characteristics on either side. It can be shown that the error in shock position location in this case is given by the square of the characteristic speed perturbation, i.e. $ ( ({\rho c - \rho _1 c_1})/{\rho _1 c_1} )^2$ . This is a higher-order effect than the simplification leading to (4.11), hence can be neglected. The analytical solution is now complete.

Figure 9 shows the prediction of the wave dynamics obtained numerically by the full hydrodynamic numerical solution and its comparison with the analytical prediction developed above. Select profiles are shown in figure 10. The analytical solution is found to be in very good agreement with the numerics, in spite of the various simplifications made. The analytical prediction nevertheless somewhat under-predicts the exact numerical solution once the inner shocks are formed. The reason is that the inner shock waves keep dissipating the heat to the local gas during the propagation. Although such heat is small comparing with the fluctuation energy, this mechanism is irreversible and will continuously heat the gas. Before the inner shocks form, this dissipation is negligible and the agreement between numerics and the analytical result is very good.

Figure 10.

Temperature profiles at select times corresponding to the profiles of figure 9.

4.4. Decay of the N-wave train

Once the shocks form, they rapidly asymptote into a train of sequential N-waves that progressively weaken as they travel away from the piston. While their speed is approximately constant, their amplitude is continuously decaying. It is of interest to establish the decay rate of their amplitude and dependence on forcing frequency, as this serves to understand the reactive solution described below and evaluate their penetration capability. Consider an N-wave at a sufficiently late time (after travelling a few wavelengths), such that the wave shape can be well approximated by a saw-tooth shape with linear profiles separating the shocks. The wavelength $\varDelta$ is controlled by the forcing frequency, i.e. $\varDelta =\rho _1 c_1/f$ . The wave evolution is given by (4.8); with $\alpha =\rho c$ , this is the inviscid Burgers equation:

(4.16) \begin{equation} \left ( \frac {\partial \alpha }{\partial t} \right )_\phi +\alpha \left (\frac {\partial \alpha }{\partial \phi }\right )_t=0. \end{equation}

In a wave fixed frame moving at speed $\alpha _1=\rho _1 c_1$ , the shock is stationary. After a time $\delta t$ , using the method of characteristics, the multi-valued solution requiring shock fitting has displaced the wave shape according to its amplitude $\alpha$ ; this is shown schematically in figure 11. Instead of the multi-valued solution denoted by the wave shape marked by the dotted lines in figure. 11, the position of the fitted shock coincides with the original one and its new peak amplitude is $\alpha _s'$ . The evolution of the peak amplitude can be easily determined. In figure 11, the slope of the linear profile at time $t+\delta t$ is $\alpha _s'/ (\varDelta /2 )$ . But since the peak has been displaced by $\alpha _s \delta t$ , the slope is also $\alpha _s/ (\varDelta /2+\alpha _s \delta t )$ . We thus have the relation

(4.17) \begin{equation} \alpha _s'/\left (\varDelta /2\right )=\alpha _s/\left (\varDelta /2+\alpha _s \delta t\right ). \end{equation}

Writing $\alpha _s'=\alpha _s+ ({\rm d}\alpha /{\rm d}t ) \delta t$ , in the limit of $\delta t \rightarrow 0$ , this becomes a differential equation for the wave amplitude:

(4.18) \begin{equation} \frac {1}{\alpha _s^2}\frac {{\rm d}\alpha _s}{{\rm d} t}=-\frac {2}{\varDelta }. \end{equation}

Its solution shows that the N-wave amplitude $\alpha _s$ decays as $\varDelta /(2t) =\rho _1c_1/ (2tf)$ .

Figure 11.

The decay of an N-wave of fixed wavelength $\varDelta$ ; red denotes the wave shape at time $t$ and blue its evolved form after time $\delta t$ . The shaded regions are the equal-area lobes justifying the location of the fitted shock at the same position as the initial one.

4.5. Corrections for shock dissipation

The analytical model described above neglects the energy dissipation of internal shocks. Our ignition simulations described below show that the cumulative effect of many internal shocks can have a non-negligible effect on ignition, particularly for cases of very high-frequency forcing, where an igniting particle undergoes repeated shocking before it ignites. It is thus of interest to correct the model above and incorporate the energy dissipation of internal shocks. To this effect, a simple correction scheme is to use the shock strength obtained in the model above, and estimate the irreversible temperature increase via the exact shock-jump equations. The shock Mach number is evaluated from the jump in particle speed across each shock obtained by the method of characteristics. The resulting temperature increment is the correction that is added to all the gas along that constant- $\phi$ line. Each successive shock along that line has the same strength and hence the same correction. For example, gas having been shocked by four internal shocks will have four times that temperature correction.

Figure 12 (left) shows the Mach number of the inner shocks evaluated as a function of particle label $\phi$ . Figure 12 (right) gives the corresponding temperature correction. Unity solution in figure 12 (left) denotes the region where a shock has not yet formed, in which region there is no shock dissipation and the correction vanishes.

Figure 12.

Evolution of the Mach number of the internal shock wave (left) and the corresponding irreversible increase across the shock (right).

Figure 13.

The irreversible temperature gain map due to inner shock waves in $2{\rm H}_2$ + ${\rm O}_2$ mixture at state 1 with $A=325$ m s⁻¹.

Figure 14.

Temperature profiles obtained numerically (solid lines) and analytically (broken lines) in $2{\rm H}_2$ + ${\rm O}_2$ mixture at state 1 with $A=325$ m s⁻¹.

The additional temperature correction is shown in figure 13. The dissipation becomes finite once the inner shocks form at approximately $\phi ^*$ = 0.00126 kg m $^{-2}$ in this case. The increasing of the fluctuation amplitude and the frequency would shorten this distance. To the right of this boundary is the highest region of shock dissipation, which decays as the inner N-wave structure decays.

The corrected temperature evolution is shown in figure 14. The temperature amplitude is now in excellent agreement with numerics.

4.6. Problem B: reflected waves on the lead shock

The disturbances generated along the C+ characteristics described above in problem A finally reflect on the lead shock, generating a wave of opposite nature on the C− characteristics (expansion for an incident compression, and vice versa) and a temperature (entropy) disturbance on the particle path. At the level of approximation considered, these reflected disturbances propagate along straight lines ${\rm d}\phi /{\rm d}t=-\rho _1 c_1$ for the C− characteristics and ${\rm d}\phi /{\rm d}t=0$ for the particle paths by definition and the lead shock path is a straight line. Appendix B provides the details of this problem for a single disturbance propagating along a C+ characteristic arriving at the lead shock. Note that the reflected acoustic disturbances are much weaker than the entropy waves, hence re-reflections between these waves with other incident waves are not treated. The treatment provided in Appendix B for a single wave can be generalised for the continuum of waves arriving at the lead shock. Say the lead shock strength at location $\phi$ is known. The C+ characteristic arriving at location $\phi +\delta \phi$ will change the shock strength and provide temperature perturbations along the C− and C0 waves. These perturbations are then added along these respective lines for the entire domain. Marching along the shock, the arrival of the next C+ characteristic provides the next perturbation and reflected perturbations along C− and C0 lines. We marched along the shock at regular intervals in order to obtain the entire interior solution. The interior solution is then interpolated on the grid with the desired resolution. The entire solution for reflected waves can thus be obtained by the method of characteristics. The temperature perturbations due to the reflected acoustic disturbances are shown in figure 15 while the temperature perturbations are shown in figure 16 for the chemically frozen 2H $_2$ + O $_2$ mixture with post-shock temperature of $T_1=1100$ K and fluctuation frequency of $f=45.35$ kHz.

Figure 15.

Temperature gain contribution due to C− reflected waves due to the reflection of C+ waves interacting with the lead shock in 2H $_2$ + O $_2$ .

Figure 16.

Temperature gain contribution due to entropy waves generated at the lead shock from C+ waves reflected at the lead shock in 2H $_2$ + O $_2$ .

These two solutions can be linearly combined along with the solution of problem A to obtain the desired solution for temperature evolution in the induction zone. The approximate solution illustrated in figure 17 is found in general good agreement with the full numerical solution of figure 6. The timing and amplitude of the various events are well reproduced. Note, however, that the central part of the hotspots are predicted to be approximately 5 $\%$ hotter than the full simulation results; likewise, cooler portions are under-predicted by the approximate model by a similar amount.

Figure 17.

Analytically reconstructed temperature field for 2H $_2$ + O $_2$ at conditions of table 1 and $A=0.2 u_{p0}$ and $f=45.4$ kHz without chemical reaction.

Figure 18.

Temperature field obtained numerically (top) and analytically (bottom) for 2H $_2$ + O $_2$ at conditions of table 1, $A=0.2 u_{p0}$ and $f=454$ kHz without chemical reaction.

The gas-dynamic model formulated revealed the most important effects to consider in the interpretation of the reactive problem, namely the propensity to form inner shocks at a finite distance from the piston. These inner shocks generate entropy and the cumulative effect may generate the strongest hotspots. Second in order of importance are the reflections of these shocks on the lead shock, generating temperature perturbations by strengthening the lead shock. Note that the inner shocks decay as N-waves as $1/t$ . In problems with very high frequency, the resulting temperature increase near the piston dominates the temperature increase away from the piston. An example is shown in figure 18 obtained numerically and using the approximate model developed. In this case, the fluctuation frequency is increased to 10 times the value as the case shown in figure 6. The strongest hotspot is dominated by the inner-shock dissipation. Overall, the analytical results predict very well the location and temperature of the hotspot.

Table 1.

Relevant thermochemical properties of the two reacting mixtures.

5.1. Fuels and their characteristic ignition and reaction times

Shock-induced ignition and transition to detonation was studied in two reactive mixtures spanning the different behaviours of ignition observed in practice, namely 2H $_2$ + O $_2$ and C $_2$ H $_4$ + 3O $_2$ at a temperature of approximately 1100 K, which is the lower-end temperature permitting sufficiently rapid auto-ignition and transition to detonation in practice. For the 2H $_2$ + O $_2$ mixture, a post-shock temperature of $T_1=1100$ K and post-shock pressure of $p_1=1$ atm are obtained by a piston with mean speed of 1626 m s⁻¹ moving into a gas at an initial temperature of $T_0=300$ K and initial pressure of $p_0=6$ kPa. This choice was also motivated by the abundance of experiments on reference ignition data, e.g. induction delay, for hydrogen–oxygen at 1 atm, as well as addressing the final stages of DDT where sufficiently strong shocks are generated in this range of temperatures. Moreover, approximately 1100 K is the lower limit of high-temperature ignition for hydrogen–oxygen ignition at 1 atm (Meyer & Oppenheim Reference Meyer and Oppenheim1971). The test conditions for C $_2$ H $_4$ + 3O $_2$ were $T_0=300$ K and $p_0=6.2$ kPa, with a mean piston speed of $u_{p0}=1259$ m s⁻¹, which brings the post-shock conditions to those studied by Saif (Reference Saif2016) in DDT experiments. Table 1 lists the relevant mixture properties under these conditions. We first characterised the relevant time scales in the two mixtures in constant-volume calculations; the piston forcing characteristics are referenced to the relevant chemical times. We define the ignition delay time $t_i$ at constant volume as the time elapsed until maximum thermicity. For a mixture of ideal gases, the thermicity reduces to (Fickett & Davis Reference Fickett and Davis1979; Williams Reference Williams1985; Kao & Shepherd Reference Kao and Shepherd2008)

(5.1) \begin{equation} \dot {\sigma }= \sum \limits _{i=1}^{N}\left (\frac {\bar {W}}{W_i}-\frac {h_i}{c_pT}\right ) \frac {{\rm D} Y_i}{{\rm D} t}, \end{equation}

where $W_i$ is the molecular weight of the $i\text {th}$ component, $\bar {W}$ is the mean molecular weight of the mixture, $h_i$ is the specific enthalpy of the $i\text {th}$ species, $c_p$ is the mixture frozen specific heat and $Y_i$ is the mass fraction of the $i\text {th}$ species in the mixture of $N$ total number of species. The other relevant chemical time is the exothermic reaction time $t_r$ , during which energy release affects the dynamics of the gas. It is defined as the inverse of the maximum thermicity. These two time scales are reported in table 1. The ethylene mixture is characterised by a much larger ratio of these two time scales $\varLambda =t_i/t_r$ . This value lies at the lower range of most hydrocarbons, which partly motivates our selection of this mixture. Very large values are computationally prohibitive, since the reaction time and length scales require sufficient resolution. Figure 19 shows the evolution of temperature and thermicity for the two fuels. For ethylene, we have also calculated the ignition process using the full mechanism. The reduced mechanism predicts well the ignition delay, with a minor under-prediction. Both models predict very short reaction times $t_r$ , the reduced model over-predicting it by a factor of 2. In the reactive calculations presented below, we ensured that the reaction zone thickness was sufficiently well resolved with at least 10 points. The grid spacing used $\varDelta_\phi$ was $2\times 10^{-6}$ kg m $^{-2}$ for the hydrogen mixture and $2.8\times 10^{-6}$ kg m^- $^{-2}$ for the ethylene mixture. A detailed grid sensitivity and convergence study is reported in Appendix C.

For reference, we also provide the calculations of the hypothetical CJ detonation and the corresponding von Neumann state. These values are tabulated in table 1. Of interest are the values of material velocity at the CJ state. Since in both cases the piston speed is larger than the material speed at the CJ state, this means that the detonation waves established by the piston speeds selected are over-driven.

Figure 19.

Temperature and thermicity evolution at constant volume for 2 ${\rm H}_2$ + O $_2$ (top) and C $_2$ H $_4$ + 3O $_2$ (bottom) taking state 1 (see table 1) as initial condition; ‘full’ and ‘reduced’ profiles denote solutions obtained with the full San Diego mechanism and the reduced mechanism used in this study.

Figure 20.

Shock-induced ignition in 2 ${\rm H}_2$ + O $_2$ (top) and C $_2$ H $_4$ + 3O $_2$ (bottom) without perturbations. The initial conditions in the pre-shock region are zero velocity and state 0, with piston velocities given in table 1. The post-shock state is state 1.

5.2. Ignition and transition to detonation without fluctuations

The dynamics of shock-induced ignition was first studied in the absence of fluctuations. Figure 20 shows the evolution of the temperature field. For the hydrogen mixture, the first ignition along the piston path occurs approximately at the constant-volume ignition time, since prior to this event, the lack of thermal evolution keeps the gas shocked by the lead shock at the constant post-shock state. The subsequent rapid acceleration of the reaction zone and formation of an internal shock on the time scale of energy deposition $t_r$ is characteristic of such dynamics; this internal dynamics is compatible with the model of Sharpe (Reference Sharpe2002). The rapid internal acceleration of the fast flame leads to the development of an internal CJ detonation wave followed by an interior Taylor wave (Fickett & Davis Reference Fickett and Davis1979). The arrival of the internal detonation wave at the lead shock transmits an over-driven detonation wave in the quiescent gas, which eventually decays towards the self-sustained over-driven solution supported by the piston speed.

The transition sequence for the ethylene mixture follows a similar sequence, but the rapid initial transition occurs faster, on time scales of $t_r$ , which are much shorter than the ignition delay time scale in this case.

5.3. Slow forcing

When the forcing period is longer than the induction delay, the sequence of events leading to auto-ignition resembles qualitatively the ignition sequence without fluctuations. An example is shown in figure 21 for the hydrogen system. The slow compression of the gas in the induction zone shortens the ignition delay to $2.5\times 10^{-5}$ s. One could modify the analysis of Sharpe (Reference Sharpe2002) to account for the non-uniform conditions in a straightforward, albeit algebraically complex manner. Note that the phase of the fluctuation is likely to be very important here. Initial cooling of the induction-zone gas would have the opposite effect to delay the ignition. The effect of the phase of the fluctuation, while interesting in its own right, is left for future study.

When the period of the fluctuation becomes comparable with the induction delay, it modifies the induction time gradient more significantly and can give rise to first ignition away from the piston face. An example is shown in figure 22 for the ethylene mixture. The gradient of ignition delay set up by the lead shock and modulated by the piston non-steadiness can be flatter than that provided by a steady shock. Its phase velocity becomes closer to the acoustic signals $\pm \rho _1 c_1$ , leading to more prompt acceleration by the Zel’dovich gradient mechanism (Sharpe & Short Reference Sharpe and Short2003) in both directions. While not shown, the slope of the C+ and C− characteristics in the induction zone is $\pm \rho _1 c_1$ , which, by virtue of the values of table 1 for the shock speed, are 2.5 times steeper than the lead shock. It can be speculated that the rapid formation of inner detonations is due to these more favourable gradients.

Figure 21.

Shock-induced ignition with 2 ${\rm H}_2$ + O $_2$ , $f=4.535\,{\rm kHz}$ , $A=0.2u_p$ .

Figure 22.

Shock-induced ignition in C $_2$ H $_4$ + 3O $_2$ with $f=0.2\,{\rm kHz}$ , $A=0.2u_{p0}$ .

5.4. Forcing on induction time scales

Figure 23.

Shock-induced ignition in $2{\rm H}_2+{\rm O}_2$ with $f=45.4$ kHz and $A=0.2 u_{p0}$ .

When the forcing frequency is further increased, the nonlinearity results in stronger lead shocks by the inert mechanism discussed above, but the phase velocity of the fast flames is now larger than the acoustic speeds, resulting in decoherence of the Zel’dovich mechanism. Figure 23 shows a striking example of this situation. The first ignition begins at $t=1.3\times 10^{-5}$ s. The spontaneous wave from the first ignition generates a forward and backward shock, which are out of phase with the fast flame. The forward-facing shock catches up to the lead shock and generates the second hotspot at $\phi =1.2\times 10^{-3}\,\rm kg\, m^-{^2}$ . The fast flame evolving from the second hotspot is now in phase with the acoustics, resulting in detonation formation. After formation of this one-dimensional detonation structure, the reaction zone length oscillates almost in phase with the fluctuation, but the reaction never decouples from the lead shock. Some weak shocks can be observed in the reacted gas with some complex forward and backward patterns. They are formed by interaction between the previous inner shock and the boundary (piston or lead shock wave). Since these shock waves are not in the reaction zone, they do not influence the detonation formation process.

Figure 24.

Shock-induced ignition in C $_2$ H $_4$ + 3O $_2$ with $f=20$ kHz and $A=0.2 u_{p0}$ .

A similar hot-spot ignition mechanism was also observed in the ethylene–oxygen system. The first ignition starts on the first hotspot at $t=2.80\times 10^{-5}$ s. Different from the hydrogen–oxygen ignition, the spontaneous wave generates two stronger shock waves and the forward-propagating shock wave keeps in phase with the reaction front until it reaches the lead shock. The detonation forms prior to the reaction zone reaching the second hotspot. The backward-propagating shock wave decouples from the reaction front and reflects on the piston prior to the reaction front.

Figure 25.

Shock-induced ignition in 2 ${\rm H}_2$ + O $_2$ with $f=454$ kHz and $A=0.2 u_{p0}$ . The purple line denotes the ignition locus prediction without coupling to the gas exothermicity predicted by coupling of (5.2) and the analytical temperature field of § 4.

5.5. Forcing periods between $t_i$ and $t_r$

When the forcing frequency is increased such that the nominal induction zone counts many forcing periods, a sequence of hotspots are observed. Figure 25 shows this interesting regime for the hydrogen–oxygen case for a frequency of $453.5\, \rm kHz$ and amplitude $A=0.2 U_{p0}$ . The higher frequency no longer significantly shortens the delay to first ignition at hotspots, since the inner-shock amplitude remains controlled by the fluctuation amplitude, which is kept constant. The slight reduction in ignition delay in this case is due to the residual repeated shock energy dissipation modelled in the previous section. Instead, a series of subcritical hotspots are generated along particle paths where internal shocks locally amplified the lead shock. These hotspots can be clearly identified as vertical particle paths originating from the confluence of internal shocks with the lead shock. Note that the higher frequency modulates the trajectory of the fast flame. The saw-tooth fast-flame trajectories are now significantly out of phase with the acoustics. This phenomenon is particularly striking in the more sensitive ethylene system illustrated in figure 26. Owing to the much larger effective activation energy of the induction kinetics, the hotspots are igniting much earlier than cold spots and the saw-tooth fast flames are even more out of phase with the acoustics.

To help us in the interpretation of the mechanism of ignition, we made use of the temperature field prior to ignition determined in § 4, in order to determine the locus of the fast flame without accounting for exothermicity. For this purpose, the induction layer is modelled in a conventional way by

(5.2) \begin{equation} \left (\frac {\partial \zeta }{\partial t}\right )_\phi =k \exp \left (-\frac {T_a}{T} \right ), \end{equation}

where $\zeta$ is the induction progress variable ranging from 0 (fresh gases) to 1 (end of induction time) and the constants $k$ and $T_a$ are calibrated to recover the correct ignition delay and sensitivity to temperature. These are calibrated for constant-volume calculations in Cantera. They yield $T_a =9562.2$ K and $k=1.6118\times 10^8\, \rm s^{-1}$ for the hydrogen–oxygen mixture and $T_a =23\,825.0$ K and $k=3.452\times 10^{12}\,\rm s^{-1}$ for the ethylene mixture under these operating conditions. Given the temperature field, integration of (5.2) yields the ignition delay $t_i(\phi )$ along a particle path implicitly:

(5.3) \begin{equation} 1=\int _{t_S}^{t_i}k \exp \left (-\frac {T_a}{T(t)} \right )\mathrm {d}t, \end{equation}

where $t_S(\phi )$ is the lead shock time.

Figure 25 shows how the sequence of hotspots are formed and compares their location with the zero-order prediction. While each sequential hotspot is formed from the coalescence of the internal shock with the lead shock, we note that the ignition delay of each hotspot is successively shorter than the previous one. This is very well predicted by the ignition model with no exothermicity coupling. The mechanism of this reduction is the cumulative energy dissipation heating along a particle path. Nevertheless, we observe that by the fourth hotspot, the energy release from previous hotspots has now played a sensible role in shortening the induction time. The inert gas-dynamic model begins to over-predict the ignition delay, since it neglects the effect of energy release on further gas compression in the induction zone. The mechanism suggested from figure 25 is through the strengthening of forward-facing shock waves passing through the main reaction zone. After the sixth hotspot, the sequence of hotspots become in phase with the motion of one of the inner shocks. We thus see a discrete-like gradient amplification mechanism where the sequence of hot-spot onset becomes in phase with acoustics. We label these hot-spot cascades.

Figure 26.

Shock-induced ignition in C $_2$ H $_4$ + 3O $_2$ with $f=200$ kHz and $A=0.2 u_{p0}$ . The purple line denotes the ignition locus prediction without coupling to the gas exothermicity predicted by coupling of (5.2) and the analytical temperature field of § 4.

More complex hot-spot cascades can be observed in a more sensitive mixture. Figure 26 shows the ethylene–oxygen ignition with fluctuation frequency of $200\, \rm kHz$ and amplitude of 20 % of the mean speed. The first ignition starts at $t=2.2\times 10^{-5}\, \rm s$ from the first hotspot – this is very well predicted by the uncoupled ignition model. Nine hot-spot ignitions can be observed prior to the reaction front merging to the lead shock and detonation formation. The second hotspot is also very well predicted by the model. The third, however, is accelerated by the energy release of the previous two, transmitted along the forward shocks originating from the first two. The three first hotspots thus follow the hot-spot cascade mechanism described above. Note that the inert model over-predicts the ignition delay more severely than for hydrogen, due to neglect of the gas-dynamic heating from previous reactions. By the fourth hotspot, the heating provided by the previous three hotspots gives rise to a more substantial decrease in the induction delay, as expected for this system with a higher activation energy.

Nevertheless, the arrival of the first shock amplified by the hot-spot cascade at the lead shock creates a strong enough disturbance to lead to prompt ignition at the hotspot along $\phi =0.006$ (hotspot 7). This localised energy release triggers the backwards sequence of promoting the ignition of hotspots 6 and 5 and 4. It also triggers the sequence of hotspots 8 and 9, aided by the shock from hotspot 3. Note that the internal transitions follow again a discrete version of the gradient or SWACER mechanism. We refer to this more complex situation as ‘bifurcated hotspot cascade’, since forward and backward hotspot avalanches are present due to hot-spot–lead-shock feedback.

5.6. Fast forcing

Once the ignition delay is much longer than the fluctuation period, the hot-spot cascade mechanism disappears. Figure 27 shows the pattern for hydrogen–oxygen ignition with frequency in $4535\,{\rm kHz}$ and amplitude of 20 % of the mean speed. The first ignition starts near the piston at $t=1.0\times 10^{-5}$ s and a slow reaction front develops in the first layer of gas of width of approximately $\delta \phi \simeq 0.0002$ kg m $^{-2}$ . Subsequently, the evolution follows a process very similar to the non-fluctuated case of figure 20. The zone marked by the slow acceleration of the reaction front is likely due to the finite region of dissipation created by the N-wave train. The N-wave train decays as $1/x$ and its influence is only felt a few wavelengths away from the piston. For reference, the inert solution can be inferred from the results shown in figure 18 by rescaling $t$ and $\phi$ by a factor of 10, controlled by the unique time scale in the inert problem, the fluctuation period. The inert N-wave train inferred by the rescaled results of figure 18 has decayed by nearly an order of magnitude at this penetration distance of $\delta \phi \simeq 0.0002$ , which corresponds to eight wavelengths of the train. The relatively slow reaction front observed in figure 27 is thus due to the energy dissipation gradient of the N-wave train decay. Further away, the fluctuations no longer play a substantial role directly. Although the first ignition delay was very short, the detonation front forms relatively late. Compared with the non-fluctuation case (figure 20), the inner supersonic reaction front forms almost at the same time. The place and time of merging between the inner supersonic reaction front and the lead shock does not have significant changes. This very high-frequency fluctuation does not make a significant contribution to reducing the time of the detonation formation in this one-dimensional ignition problem.

Figure 27.

Shock-induced ignition in 2 ${\rm H}_2$ + O $_2$ with $f=4540$ kHz and $A=0.2 u_{p0}$ .

Figure 28.

Shock-induced ignition in C $_2$ H $_4$ + 3O $_2$ with $f=2000$ kHz and $A=0.2u_{p0}$ .

A similar pattern was observed for the more sensitive ethylene–oxygen system. Figure 28 shows the ignition diagram with fluctuation frequency of $2000\, \rm kHz$ and amplitude of 20 % of the mean speed. In this case, the ignition delay has been shortened at $t=1.3\times 10^{-5}\, \rm s$ but the supersonic reaction wave takes longer to form. Two hotspots are formed in this case. The first hotspot is generated from the dissipation. The second hotspot is formed by the complex hotspot cascade mechanism which is the interaction between the lead shock and intensified inner shock. In this mixture, the ignition delay is much more sensitive than for the hydrogen–oxygen mixture. The intensified shock wave is strong enough to generate an ignition hotspot. This hotspot ignition locus quickly accelerates to an inner detonation, which catches up to the main shock. These inner dynamics shortens the formation distance of the final detonation by approximately 90 % as compared with the nominal non-fluctuated case and cannot be neglected.