1 Introduction
Detonation waves are combustion waves that move at supersonic speed with respect to the energetic medium through which they propagate. For a wave propagating at constant velocity, the steady conservation laws of mass, momentum and energy can be applied to a control volume that encloses and moves with the wave. Along with an equation of state, an additional criterion is required to close the set of governing equations and solve for the wave velocity and thermodynamic state at the exit of the wave. The condition of sonic flow (with respect to the wave-fixed frame) was proposed by Chapman (Reference Chapman1899) and Jouguet (Reference Jouguet1905), and this state can be shown to correspond to the minimum possible speed of a compressive wave with equilibrium products that satisfies the steady conservation laws. The Chapman–Jouguet (CJ) criterion can be justified on the grounds that subsonic flow at the exit of the control volume would allow rarefaction waves overtaking the detonation from downstream to weaken and decelerate the wave. Thus, a detonation overdriven by a piston from behind with subsonic flow at the end of the reaction zone (a so-called strong detonation) would eventually relax to the CJ solution when the driving piston is stopped. A supersonic state exiting a self-sustained wave (a weak detonation) might be possible under the scenario of complex chemical kinetics wherein the reaction zone is described by an exothermic phase, which brings the non-equilibrium flow to sonic, followed by an endothermic phase. In these so-called pathological detonations, rarely observed in experiments, the CJ criterion is replaced by an equivalent generalized CJ criterion that applies to the non-equilibrium flow in the reaction zone (Eyring et al. Reference Eyring, Powell, Duffy and Parlin1949; Wood & Kirkwood Reference Wood and Kirkwood1954; Higgins Reference Higgins and Zhang2012). The classical CJ criterion has been remarkably successful in predicting detonation velocity, with experimentally observed detonation waves in gases usually observed to propagate within 1 % of the equilibrium CJ speed in nearly all gaseous mixtures. The CJ criterion has had similar success in condensed phase energetic materials (i.e. solid and liquid high explosives), although uncertainty in the equation of state of condensed phase materials at extremely high pressure means that these calculations are not usually conducted from first principles and must rely on empirical correlations (often using the measured detonation velocity of the explosive itself as an input parameter).
The success of the CJ criterion is all the more remarkable given that all known detonations in gases are unstable, dominated by a transient and multi-dimensional cellular structure. The fact that detonations, governed by activated chemical reactions (i.e. Arrhenius kinetics), are unstable to perturbations was first demonstrated by Erpenbeck (Reference Erpenbeck1964). Short & Stewart (Reference Short and Stewart1998) proved that a detonation with a realistic value of heat release and activation energy would be unconditionally unstable in two dimensions, resulting in a cellular structure. The experimental evidence of detonation cells was first noted by Denisov & Troshin (Reference Denisov and Troshin1961), and these detonation cells are now recognized as being essential for the propagation of detonations in high activation energy mixtures (e.g. hydrocarbon fuels in air). In mixtures with high activation energy (e.g. methane/oxygen/diluent (Strehlow Reference Strehlow1968, Reference Strehlow1969)) or with gases featuring a large number of molecular degrees of freedom, resulting in low values of heat capacity ratio (e.g. halocarbons (Vandermeiren & van Tiggelen Reference Vandermeiren and van Tiggelen1988; Lefebvre, Nzeyimana & van Tiggelen Reference Lefebvre, Nzeyimana and van Tiggelen1993)), the detonation cellular structure becomes highly irregular, with explosive bursts of fine cellular structure and regions of apparently no cells at all. In such mixtures, deviations in local detonation front velocity can be as much as 50 % below and 100 % above the mean velocity. Despite these complexities, the CJ criterion is still remarkably accurate in predicting the overall average speed of the detonation front. In liquid explosives and some solid explosives such as trinitrotoluene (TNT), similar cellular instabilities are also noted (Fickett & Davis Reference Fickett and Davis2000; Lee Reference Lee2003). Other explosives, such as those typically used in commercial blasting, have a highly heterogeneous structure due to the coarse grains of energetic material used, with the scale of the heterogeneity being of the order of the detonation front thickness. Again, in these cases, the CJ criterion appears to be remarkably successful at predicting the detonation velocity of effectively infinite diameter charges (although the earlier caveat about equation of state still applies).
The seemingly paradoxical situation, in which the steady, one-dimensional CJ criterion appears successful in describing the propagation of what is now known to be a highly complex wave structure, was well articulated by Oppenheim shortly after the discovery of cellular structure in gases:
The detonation may form an essentially non-steady, non-uniform regime so that in order to explain its precise nature, multi-dimensional effects in space as well as its irregular behavior in time have to be taken into account. In fact, in view of this evidence one should express amazement that the one dimensional steady flow theory was so successful in rationalizing so many experimental observations (Oppenheim Reference Oppenheim1961).
A half-century later, this puzzle has not been completely resolved. The usual argument that is made is that, if a sufficiently large control volume is used to enclose the wave, then the steady, one-dimensional control volume analysis must still apply on average, even though the wave itself is multi-dimensional and unsteady. This reasoning is not entirely satisfactory, however. For one thing, it is not clear where to place the exit plane of the control volume a priori (too close to the initial shock front will locate the exit plane in a non-equilibrium flow while too far downstream places the exit plane in an unsteady rarefaction, the Taylor wave, where the steady conservation laws do not apply). In addition, for an unsteady flow, the concept of sonic flow is not particularly relevant (i.e. it is unclear which reference frame should be used in comparing the flow velocity to the local speed of sound). For unsteady detonation propagation in one dimension, the sonic surface is more accurately replaced by the limiting separatrix of characteristic waves propagating from behind the detonation toward the shock front, with characteristics in front of the separatrix able to influence the shock and those behind not able to reach the shock (Kasimov & Stewart Reference Kasimov and Stewart2004; Ng & Zhang Reference Ng and Zhang2012). The applicability of the steady, one-dimensional analysis to a two-dimensional, moderately unstable detonation wave was explored computationally by Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007) via spatial and temporal averaging of the unsteady structure in a steady reference frame following the wave. In their study, the effect of the interplay between the chemical exothermicity and hydrodynamic fluctuations on the location of the effective sonic surface in a two-dimensional, unsteady detonation structure was clearly illustrated by the resulting one-dimensional averaged wave profile (Radulescu et al. Reference Radulescu, Sharpe, Law and Lee2007).
Schematic of an explosive medium with (a) source energy (red dots) uniformly distributed throughout inert medium (blue dots) and with (b) source energy collected into planar sheets (layers of red dots of width
$W$
and spacing
$L$
) with inert media in between.
In this paper, we intend to put the CJ criterion to a rigorous test by examining detonation in a media in which the energy release is condensed into spatially discrete layers separated by regions of inert gas (as illustrated in figure 1). The energetic layers will release their energy when activated by a passing shock front after a fixed delay time. The new blast wave generated by the energy release becomes the mechanism to initiate additional sources, such that the wave propagates via a sequence of ‘sympathetic detonations’. The spatial scale of the sources and their spacing will be varied, examining the behaviour of the wave propagation in this system from continuous energy release (i.e. no discretization) to highly concentrated sources in the limit of
$\unicode[STIX]{x1D6FF}$
-function-like sources of energy separated by layers of inert gas. The overall energy release will be maintained as constant, so that the average wave velocity can always be compared to the equivalent CJ speed of the homogenized media. This work builds upon similar research examining flame propagation in discretized media, in which point-like sources release heat to diffuse outward and trigger subsequent sources (Goroshin, Lee & Shoshin Reference Goroshin, Lee and Shoshin1998; Beck & Volpert Reference Beck and Volpert2003; Tang, Higgins & Goroshin Reference Tang, Higgins and Goroshin2009; Goroshin, Tang & Higgins Reference Goroshin, Tang and Higgins2011; Tang, Higgins & Goroshin Reference Tang, Higgins and Goroshin2012; Mi et al.
Reference Mi, Higgins, Goroshin and Bergthorson2017). A wave propagation mechanism by which an energetic source generates a blast that, in turn, can initiate the next source was proposed by Stewart & Asay (Reference Stewart, Asay, Fife, Liñán and Williams1991) in modelling the response of propellant beds comprised of explosive grains to strong shock stimuli. They outlined a ‘theory of discrete interactions’ in which initiation of subsequent sources would be described by a nonlinear recursion relation. In spray detonations with millimetre-sized droplets of fuel, a similar propagation mechanism, i.e. blast waves generated by each individual droplet releasing its energy, merging with and reinforcing the leading shock wave, was experimentally identified by Dabora, Ragland & Nicholls (Reference Dabora, Ragland and Nicholls1969) and theoretically described by Pierce & Nicholls (Reference Pierce and Nicholls1973). Propagation of detonation in media with a sinusoidal variation in properties was explored computationally by Morano & Shepherd (Reference Morano and Shepherd2002) in one dimension and extended to two dimensions by Li, Mi & Higgins (Reference Li, Mi and Higgins2015). In these studies, no significant deviation greater than 2 % away from CJ was observed when the detonation propagated without losses present. Li et al. (Reference Li, Mi and Higgins2015) did observe enhanced wave speeds when the detonation propagated in a layer with inhomogeneities and yielding confinement (in comparison to a detonation propagating in a homogeneous layer with the same yielding confinement). The enhanced wave speeds were attributed to accelerated reaction rates resulting from the shock interactions caused by the inhomogeneities. Detonation propagation in a media described by the Euler equations with
$\unicode[STIX]{x1D6FF}$
-function-like sources has not been previously published.
The plan for this paper is as follows. In § 2, the problem statement and parameters will be defined. Section 3 describes the numerical methodology used to solve the governing Euler equations. Section 4 presents the results of a study where each of the model parameters was systematically varied and compared to the ideal CJ detonation velocity of the equivalent homogenous media. Section 5 analyses the results by temporally averaging select simulations to provide a quasi-steady one-dimensional description of the transient waves. The results are discussed in § 6 and summarized in the Conclusions (§ 7). An appendix presents an ad hoc, or heuristic, model that attempts to construct an analytic solution to this problem for the case of zero delay in triggering the sources.
2 Problem statement
This study considers a calorically perfect gas (i.e. fixed ratio of specific heats
$\unicode[STIX]{x1D6FE}$
) that has the potential to release energy with a heating value of
$\widetilde{Q}$
(
$\text{J}~\text{kg}^{-1}$
). The tilde ‘
${\sim}$
’ denotes a dimensional quantity. The flow variables, density, pressure and particle velocity are non-dimensionalized with reference to the initial state ahead of the leading shock, i.e.
$\unicode[STIX]{x1D70C}=\widetilde{\unicode[STIX]{x1D70C}}/\widetilde{\unicode[STIX]{x1D70C}}_{0}$
,
$p=\widetilde{p}/\widetilde{p}_{0}$
, and
$u=\widetilde{u}/\sqrt{\widetilde{p}_{0}/\widetilde{\unicode[STIX]{x1D70C}}_{0}}$
, respectively, and the space coordinate
$x$
by the spacing between two adjacent sources
$\widetilde{L}$
, i.e.
$x=\widetilde{x}/\widetilde{L}$
. The subscript ‘
$0$
’ indicates the initial state of the uniform reactive medium. The heat release
$\widetilde{Q}$
is non-dimensionalized as
$Q=\widetilde{Q}/\left(\widetilde{p}_{0}/\widetilde{\unicode[STIX]{x1D70C}}_{0}\right)$
. A detonation wave propagating through this uniform medium is expected to move at the non-dimensional CJ speed, given by,
$$\begin{eqnarray}V_{CJ}=M_{CJ}c_{0}=\left\{\sqrt{\frac{\unicode[STIX]{x1D6FE}^{2}-1}{\unicode[STIX]{x1D6FE}}Q+\sqrt{\left(\frac{\unicode[STIX]{x1D6FE}^{2}-1}{\unicode[STIX]{x1D6FE}}Q+1\right)^{2}-1}+1}\,\right\}\sqrt{\unicode[STIX]{x1D6FE}},\end{eqnarray}$$
where
$c_{0}$
denotes the non-dimensionalized initial speed of sound and equals
$\sqrt{\unicode[STIX]{x1D6FE}}$
. Equation (2.1) is the classical CJ detonation solution, and is used to compare the resulting propagation speeds from the discrete source simulations with that of the steady detonation in the equivalent homogeneous media. All the simulation results of wave speeds reported in this paper are normalized by the corresponding
$V_{CJ}$
.
To examine the effect of spatial discretization, the energy is initially collected in regions with width
$\widetilde{W}$
(source width) that are spaced a distance
$\widetilde{L}$
apart (i.e. the distance
$\widetilde{L}$
is, for example, from the right edge of a source to the right edge of its neighbour). The total energy release per unit mass of the medium remains fixed at
$Q$
, so the energy in each concentrated source is
$Q\widetilde{L}/\widetilde{W}$
. The initial discreteness of the system is described by a parameter
$\unicode[STIX]{x1D6E4}_{init}=\widetilde{W}/\widetilde{L}$
. As
$\unicode[STIX]{x1D6E4}_{init}\rightarrow 0$
, the sources of energy become highly spatially concentrated. In the limit of
$\unicode[STIX]{x1D6E4}_{init}\rightarrow 1$
, the energy release becomes continuous through the medium. In most of the simulations reported in this paper, the following mechanism of energy deposition is incorporated (unless otherwise specified). A discrete source of energy is instantaneously triggered by the passage of a rightward-propagating shock wave past the right edge of the source, independent of the strength of the shock, after a non-dimensional time delay
$t_{d}$
(time variables are non-dimensionalized as
$t=\widetilde{t}/(\widetilde{L}/\sqrt{\widetilde{p}_{0}/\widetilde{\unicode[STIX]{x1D70C}}_{0}})$
). This time delay is reported as another non-dimensional value
$\unicode[STIX]{x1D70F}=t_{d}V_{CJ}$
, where the time required for the CJ detonation to propagate over the distance
$\widetilde{L}$
in the equivalent homogenous system has been used to non-dimensionalize the delay time. Note that the non-dimensionalized spacing between two adjacent sources is unity.
Three different scenarios of sources are considered: fixed sources that deposit their energy instantaneously, convected sources that deposit their energy instantaneously and Arrhenius-rate-governed media in which the energy evolves from the compressible media itself. For fixed sources, the energy source is considered to be independent of the inert, gaseous medium and is not convected along with it. Conceptually, this would correspond to a fine mesh of explosive wire in a tube filled with inert gas, wherein the mesh spans the tube cross-sectional area and is attached to the tube wall. In this case, the shock-accelerated, inert gas is free to flow (without resistance) around the sources of energy, which remain fixed until they release their energy. In the second scenario (convected sources), the sources are assumed to be embedded in the inert, gaseous medium and convect along with it. This would correspond to a layer of explosive dust, suspended between sections of inert gas. As the successive blast waves propagate through the medium, the sources in this scenario are compressed and convected along with the post-shocked inert gas. The two above-described scenarios are illustrated schematically in figure 2. Finally, the case wherein the reactive Euler equations are used (with an Arrhenius-rate controlling the heat release) is considered. This would conceptually correspond to a scenario in which layers of explosive gas are separated by layers of inert gas and convected with the flow.
For the cases with instantaneous energy deposition, after the delay time has elapsed, the energy that is released by the source is deposited in the volume occupied by the source at that moment. Maintaining the overall amount of specific heat release of the medium at
$Q$
, the conservation laws require that the energy deposition results in an increase in the pressure in the region of an activated source by an amount given by,
$$\begin{eqnarray}\unicode[STIX]{x0394}p=\frac{Q\left(\unicode[STIX]{x1D6FE}-1\right)}{\unicode[STIX]{x1D6E4}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E4}$
is the effective spatial discreteness of a source when it instantaneously releases energy. It is important to note that, in the scenario of fixed sources, the source width remains constant as being processed by the leading shock and hence
$\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{init}$
; in the scenario of convected sources, the actual source width
$\widetilde{W}$
at the instant of energy release is smaller than its initial value, i.e.
$\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{init}$
, due to the compression of the medium in which the source energy is stored and the particle motion induced by the passage of the shock wave. The numerical implementation of delay time and convected source tracking and energy addition is described in § 3. This method of adding energy to a compressible flow is essentially the same as that originally used by Brode (Reference Brode1955) to numerically simulate the Taylor–Sedov point-blast problem.
Schematic of the numerical implementation of the energy release process and source tracking algorithm for (a) Scenario 1, laboratory-fixed and (b) Scenario 2, convected discrete sources at three different moments: (i) before the arrival of the leading shock wave
$x_{s}(t)$
, (ii) the source being shocked marking the onset of delay period, and (iii) the source instantaneously releasing energy.
The gas dynamics of the proposed system coupled with a mechanism of instantaneous energy deposition after a delay time can be described by the one-dimensional Euler equations with a source term
$\boldsymbol{S}$
in the laboratory-fixed reference frame,
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\boldsymbol{F}\left(\boldsymbol{U}\right)}{\unicode[STIX]{x2202}x}=\boldsymbol{S}\left(\boldsymbol{U}\right),\end{eqnarray}$$
where the conserved variable
$\boldsymbol{U}$
and the convective flux
$\boldsymbol{F}$
are, respectively,
$$\begin{eqnarray}\boldsymbol{U}=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}\\ \unicode[STIX]{x1D70C}u\\ \unicode[STIX]{x1D70C}e\end{array}\right),\quad \boldsymbol{F}\left(\boldsymbol{U}\right)=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}u\\ \unicode[STIX]{x1D70C}u^{2}+p\\ (\unicode[STIX]{x1D70C}e+p)u\end{array}\right).\end{eqnarray}$$
In the above equations,
$e=p/(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}+u^{2}/2$
is the non-dimensional specific total energy.
For the cases where energy is deposited instantaneously, the energy is incorporated in the energy balance equation via
$\boldsymbol{S}\left(\boldsymbol{U}\right)$
as follows,
$$\begin{eqnarray}\boldsymbol{S}\left(\boldsymbol{U}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\\ \displaystyle \mathop{\sum }_{n}\unicode[STIX]{x0394}p\unicode[STIX]{x1D6FF}\left(t_{i,n}-t_{d}\right)H\left(x-x_{L,n}\right)H\left(x_{R,n}-x\right)\end{array}\right),\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}$
and
$H$
denote the
$\unicode[STIX]{x1D6FF}$
- and Heaviside functions, respectively and
$n$
is the index of each discrete source of energy. The positions of the left and right edges of the
$n$
th source are denoted as
$x_{L,n}$
and
$x_{R,n}$
, respectively. The two Heaviside functions, i.e.
$H(x-x_{L,n})$
and
$H(x_{R,n}-x)$
, describe that the energy is deposited in the volume occupied by the
$n$
th source. Note that, in the case with extremely discrete sources, where
$(x_{R,n}-x_{L,n})=W$
tends to
$0$
, the spatially discrete sources approach the limit of a
$\unicode[STIX]{x1D6FF}$
-function in space, i.e.
$\lim _{W\rightarrow 0}H(x-x_{L,n})H(x_{R,n}-x)/W=\unicode[STIX]{x1D6FF}(x-x_{n})$
. For convected sources, the fluid elements marking the left and right edges of the
$n$
th source move at their corresponding particle velocities, hence, their motions are governed by the following equations,
$$\begin{eqnarray}\frac{\text{d}x_{L,n}}{\text{d}t}=u\left(x_{L,n},t\right)\quad \text{and}\quad \frac{\text{d}x_{R,n}}{\text{d}t}=u(x_{R,n},t).\end{eqnarray}$$
A variable,
$t_{i,n}$
, with its initial value 0, is introduced in (2.5) to measure the time elapsed since the source has been shocked. The time evolution of
$t_{i,n}$
is described by the following equation,
$$\begin{eqnarray}\frac{\text{d}t_{i,n}}{\text{d}t}=H\left(p\left(x_{R,n},t\right)-p_{0}\right).\end{eqnarray}$$
The Heaviside function indicates that the measurement of
$t_{i,n}$
starts once the local pressure at the right edge of the source becomes incrementally greater than the initial pressure of the medium.
The final case considered uses a more realistic, state-dependent mechanism of energy deposition, e.g. Arrhenius chemical kinetics, instead of an instantaneous energy addition after a prescribed delay time. In order to incorporate such a reaction mechanism, a reaction progress variable or normalized concentration of the reactant,
$Z$
, which evolves from 1 (unreacted) to 0 (fully reacted), associated with the overall amount of energy release
$Q$
, is added into the formulation of specific total energy as
$e=p/(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}+u^{2}/2+ZQ/\unicode[STIX]{x1D6E4}_{init}$
. An additional (a fourth) equation accounting for the balance of
$Z$
must be coupled with the system of governing equations ((2.3) and (2.4)). The fourth components of
$\boldsymbol{U}$
,
$\boldsymbol{F}$
and
$\boldsymbol{S}$
are respectively as follows,
$$\begin{eqnarray}U_{4}=\unicode[STIX]{x1D70C}Z,\quad F_{4}=\unicode[STIX]{x1D70C}uZ,\quad S_{4}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FA},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}$
is the reaction rate, i.e.
$\unicode[STIX]{x1D6FA}=\unicode[STIX]{x2202}Z/\unicode[STIX]{x2202}t$
. Note that, in this case, no source terms explicitly appear in the equations of mass, momentum and energy balance (
$S_{1}=0$
,
$S_{2}=0$
and
$S_{3}=0$
). The reaction rate
$\unicode[STIX]{x1D6FA}$
can be explicitly expressed in a single-step Arrhenius form as follows,
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}=-kZ\times \text{exp}\left(-\frac{E_{a}}{T}\right),\end{eqnarray}$$
where
$k$
and
$E_{a}$
are the dimensionless pre-exponential factor and activation energy, respectively, and dimensionless temperature is
$T=p/\unicode[STIX]{x1D70C}$
for an ideal gas. A summary of the cases and the range of parameters considered is given in table 1. For convenience, the numbers 1, 2 and 3 (in the first column of table 1) are used as indices of the scenarios of the discretized energy sources simulated in this study.
Summary of the scenarios and the range of parameters considered in this study.
3 Numerical method
Two independently written, one-dimensional, finite-volume Euler codes were used in this study. One of them was based on a fixed (non-adaptive) uniform grid. This solver, however, encountered computer memory limitations as the sources were progressively made more discrete due to the requirement for a minimum number of computational cells within each source to properly resolve all wave phenomena induced by energy release (see discussion below). For this reason, a second code using adaptive mesh refinement was used to both extend the results into the limit of highly spatially discrete sources and verify the results of the fixed grid code. Both codes solve the one-dimensional Euler equations using the monotonic upstream-centred scheme for conservation laws (MUSCL)-Hancock total variation diminishing (TVD) node-centred Godunov-type finite-volume scheme (Toro Reference Toro2009) with an exact Riemann solver and the van Leer non-smooth slope limiter. The scheme is of second-order accuracy in space and time on smooth solutions.
In the adaptive code, the background (initial) grid is also uniform. Grid adaptation is performed at each time step via hierarchical h-refinement. Each refinement level reduces the local grid step by factor two. The number of refinement levels is determined from the requirement of having the smallest grid step corresponding to the desired number of grid nodes within the energy source (see discussion below). The refinement/coarsening sensor is based on normalized second derivatives of density and pressure. Prior to energy release via pressure modification for a source, the grid is refined in the vicinity of the source, regardless of the sensor values in the region, so that the energy release and subsequent induced wave motion would be always resolved on the finest mesh. The data structure and adaptation procedure of the adaptive code were adopted, with suitable simplification, from the two-dimensional unstructured Euler code described by Saito et al. (Reference Saito, Voinovich, Timofeev, Takayama and Toro2001) to a one-dimensional code. Preliminary numerical trials demonstrated that simple linear solution interpolation used in the code (Saito et al. Reference Saito, Voinovich, Timofeev, Takayama and Toro2001) when inserting a new node led to significant (of the order of 5 %) error in wave speeds when long distance propagation typical for the problem under consideration was simulated. For this reason, node insertion and deletion procedures were modified to be fully conservative, i.e. the conservation of mass, momentum and energy were strictly enforced when assigning gas dynamic parameters at a newly inserted node or removing a node.
In both fixed and adaptive mesh refinement codes, the instantaneous energy deposition (in Scenarios 1 and 2) by each discrete source after a finite delay period is implemented as follows. The trajectory of the leading shock front,
$x_{s}(t)$
, is tracked in the simulations by finding the location where
$p$
first increases to 1.01 from its upstream initial state
$p_{0}=1$
at each time step. A discrete source is considered as being shocked once the leading shock passes the right edge of the space occupied by this source, i.e. the pressure at the right edge of the source
$p\left(x_{R},t\right)$
becomes greater than
$p_{0}$
(here and below the source index
$n$
is omitted). The time elapsed since the source has been shocked,
$t_{i}$
, is measured by numerically integrating equation (2.7). The algorithm of delay time update at each time step
$\unicode[STIX]{x0394}t$
is formulated as follows,
$$\begin{eqnarray}t_{i}^{t+\unicode[STIX]{x0394}t}=t_{i}^{t}+\unicode[STIX]{x0394}t_{i},\quad \text{where}\quad \unicode[STIX]{x0394}t_{i}=\left\{\begin{array}{@{}ll@{}}0, & p\left(x_{R},t\right)=p_{0}\\ \unicode[STIX]{x0394}t, & p\left(x_{R},t\right)>p_{0}.\end{array}\right.\end{eqnarray}$$
Once
$t_{i}$
reaches
$t_{d}$
, the simulation is stopped before advancing to the next time step, and energy is deposited within the space occupied by this triggered source, i.e. between the left and right edges of this source (
$x_{L}\leqslant x\leqslant x_{R}$
), via increasing the pressure by
$\unicode[STIX]{x0394}p$
,
$$\begin{eqnarray}p_{updated}\left(x,t\right)=p\left(x,t\right)+\unicode[STIX]{x0394}p\quad \text{at}\quad t_{i}=t_{d},\end{eqnarray}$$
with the value of
$\unicode[STIX]{x0394}p$
being given by (2.2). Unlike fixed sources in Scenario 1, the location and volume occupied by a convected source change with time in Scenario 2. The trajectories of the fluid elements marking the right and left edges of each convected source, i.e.
$x_{R}(t)$
and
$x_{L}(t)$
, respectively, are tracked via numerical integration of the corresponding particle paths described by (2.6a,b
) using a Lagrangian approach,
$$\begin{eqnarray}x_{R}^{t+\unicode[STIX]{x0394}t}=x_{R}^{t}+\unicode[STIX]{x0394}t\,u_{R}^{t}\text{and}\quad x_{L}^{t+\unicode[STIX]{x0394}t}=x_{L}^{t}+\unicode[STIX]{x0394}t\,u_{L}^{t}.\end{eqnarray}$$
With this source tracking algorithm, the code is able to add pressure into the space occupied by a convected source at the moment when it releases energy after the delay period. The numerical implementation of energy release and source tracking algorithm for laboratory-fixed and convected sources is illustrated in figure 2.
For the cases with Arrhenius kinetics (Scenario 3), the method of source term splitting is adopted for the numerical integration to treat separately the hydrodynamic process ((2.3) and (2.4)) and the reactive process ((2.8) and (2.9)). The Strang splitting method (Strang Reference Strang1968) is used in order to maintain second-order accuracy.
The length of the entire simulated problem is over 100 spatial units, i.e. containing more than 100 discrete sources. In order to make the simulation time efficient, the computation, at every time step, is only performed in a window enclosing the leading wave complex in a laboratory-fixed reference frame, instead of the entire simulation domain. A minimum eight unit wide window (i.e. containing eight sources) is necessary to capture all of the dynamics contributing to the propagation of the leading wave front for the longest delay times considered in this study. Once the leading front reaches the end of this computational window, the window frame (i.e. left and right boundaries) is advanced by one spatial unit while maintaining the same width of the domain. A transparent boundary condition is applied on both boundaries of the computational window (Toro Reference Toro2009). In order to rule out any possible influence of the boundary conditions, a 10 unit wide window was used in all the simulations reported below. As the numerical resolution is increased from 50 to 200 computational cells within a discrete source, the simulation results do not show any significant difference. The results reported in this paper were all performed with a numerical resolution of 100 computational cells per source.
Pressure profiles showing evolution of flow field with instantaneous release of source energy (
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=1.5$
,
$Q=50$
,
$\unicode[STIX]{x1D6FE}=5/3$
and laboratory-fixed sources): (a) triggering of the first three sources and (b) later time evolution after triggering 280 sources. The symbols plotted on the horizontal axis indicate the location of energy sources. An open diamond represents a source before being shocked, a solid diamond a source undergoing the delay period after being shocked and a solid circle a source after releasing its energy.
4 Results
The results of a sample calculation are shown in figure 3, showing the pressure profile of the computational domain for a simulation of Scenario 1 with the following parameters:
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=1.5$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=5/3$
. Both early time (the first three sources) and later time (after the triggering of 280 sources) are shown in parts (a) and (b), respectively, and symbols indicate the location and status of the sources (i.e. unshocked, shocked with delay time elapsing and energy released). The local shock dynamics following the initial release of energy resembles a classical blast wave profile. At later times (figure 3
b), a large number of sawtoothed waves, which are residual, decaying blast waves from earlier sources, can be seen superimposed upon the flow field. Figure 4 shows
$x$
–
$t$
diagrams constructed from simulation results (
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=0.1$
,
$Q=50$
,
$\unicode[STIX]{x1D6FE}=1.1$
and sources held fixed), with the shading indicating pressure in the flow field obtained over the region from the 280th to the 284th source. The dark V-shaped regions in the pressure fields are the blast waves generated by the energy release of the sources, propagating both forward and backward in the flow. The forward propagating blast generated by a newly triggered source catches up and accelerates the leading shock ahead of it. After the leading shock has travelled a sufficiently long distance, the process of new sources releasing their energy and the forward-propagating blast waves catching up to the leading shock becomes nearly periodic. In figure 4(b), the spatial coordinates are transformed into a reference frame moving at the average velocity of the leading shock, i.e.
$x^{\prime }=x-V_{avg}t$
. Right-running characteristics are constructed in figure 4(b) by integrating an ordinary differential equation (ODE) describing the path along which acoustic signals move through the domain at the particle velocity (with respect to the wave-attached frame) plus the local sound speed. The use of
$x$
–
$t$
diagrams and characteristics analysis was first explored by McVey & Toong (Reference Mcvey and Toong1971) and recently by Kasimov & Stewart (Reference Kasimov and Stewart2004) and Leung, Radulescu & Sharpe (Reference Leung, Radulescu and Sharpe2010) to illustrate the physical mechanisms governing the instabilities in hypersonic exothermic flow and pulsating detonations, respectively. This construction of the
$x$
–
$t$
diagram is a post-processing exercise performed upon the already computed solution. In figure 4(b), there is a characteristic (plotted as a thickened line) such that the characteristics upstream (on the right) of it can eventually reach the leading shock, while those downstream (on the left) can never reach the leading shock. This limiting characteristic is known as the separatrix, such that only the flow field upstream of it is able to influence the shock front. It can be clearly seen in figure 4(b) that the separatrix oscillates around a location which is significantly downstream from the loci of sources releasing energy.
$x$
–
$t$
diagram of flow field in (a) laboratory-fixed and (b) wave-fixed reference frame (
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=0.1$
,
$Q=50$
,
$\unicode[STIX]{x1D6FE}=1.1$
and laboratory-fixed sources). The curves in (b) are right-running characteristics generated by integrating the trajectory of an acoustic wave through the already computed flow field, and the thickened curve indicates the separatrix. The symbols plotted on the horizontal axis indicate the location of energy sources.
Figure 5(a) shows the evolution of the instantaneous velocity,
$V_{inst}$
, normalized by
$V_{CJ}$
, at which the leading wave front propagates over the first
$8$
laboratory-fixed sources for the case of
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
. As described in § 3, the trajectory of the leading shock front is tracked in the simulations.
$V_{inst}$
as a function of
$x_{s}$
can thus be obtained by numerically differentiating
$x_{s}(t)$
using a central difference scheme.
$V_{inst}$
exhibits large fluctuations as the successive sources are triggered and the blast waves they drive reach the leading shock front. Each time a blast wave reaches and merges with the shock front,
$V_{inst}$
is brought to a peak value, and then decays to a local minimum velocity just before the next blast wave reaches the front. In figure 5(a), it can be observed that both the local minimum and maximum
$V_{inst}$
, just before and after the instant of the next blast wave reaching the front respectively, monotonically increase to steady values as the sources are successively triggered. Also shown in figure 5(a) is the prediction of a simple heuristic model that is based on the classic similarity solution of Taylor and Sedov for planar blast waves, the concept of energy partitioning of new source energy released at the blast front, and displacement of the flow field resulting from the residual influence of prior sources; see the appendix for a complete development of this model.
The history of (a) the instantaneous wave velocity normalized by
$V_{CJ}$
over the first eight laboratory-fixed sources for the case of
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
as a function of the leading shock position
$x_{s}$
and (b) the source-to-source average wave velocity over 100 laboratory-fixed sources for the cases of
$\unicode[STIX]{x1D6E4}=0.001$
to 1,
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
as a function of the leading shock position
$x_{s}$
. The heuristic model prediction is plotted as a dashed curve. In (b), the symbols are plotted only for the first 10 sources.
With
$x_{s}(t)$
obtained from the simulation results, the average velocity at which the leading wave front travels from one source to the next can be calculated. Figure 5(b) shows the evolution of the source-to-source average wave velocity,
$V_{avg,source}$
, as a function of the distance travelled by the leading wave front. The results shown in this figure are for the cases of various values of discreteness (
$\unicode[STIX]{x1D6E4}=0.001$
to continuous),
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
over
$100$
laboratory-fixed sources. For all values of
$\unicode[STIX]{x1D6E4}$
,
$V_{avg,source}$
gradually increases and converges to a nearly constant value. A plateau value of
$V_{avg,source}$
indicates that the dynamics of the shock front triggering the energy release of the sources and the resulting blast waves reaching and supporting the shock front have reached a quasi-periodic state with a constant time interval between the triggering of subsequent sources. As the sources become more discrete (
$\unicode[STIX]{x1D6E4}\rightarrow 0.001$
), the plateau
$V_{avg,source}$
converges to a value that is 15 % in excess of the CJ velocity of the equivalent homogeneous system (note the horizontal line indicating CJ velocity) for
$\unicode[STIX]{x1D6FE}=1.1$
. The distance required for
$V_{avg,source}$
to reach a quasi-steady value decreases from of the order of 100 for a continuous energy source to that of the order of 10 for
$\unicode[STIX]{x1D6E4}=0.001$
. The heuristic model (see the appendix) prediction of how
$V_{avg,source}$
evolves as sources are triggered successively is also shown in figure 5(b). In order to show the results for a large range in the value of source discreteness, the simulation results shown in figure 5(b) are those obtained by using the adaptive mesh refinement code.
The following subsections systematically study the average velocity of the wave as the various parameters in the model are varied. The influence of the source scenarios (laboratory-fixed versus convected sources) and energy deposition mechanism (instantaneous deposition versus single-step Arrhenius kinetics) is also investigated. All average velocities, subsequently reported in this paper and denoted as
$V_{avg}$
, were measured after the leading shock front had reached a quasi-steady value. The length of the simulation domain, i.e. the total number of discrete sources, over which the detonation wave needs to propagate to reach a quasi-periodic state, varied approximately from 20 to 100 as
$\unicode[STIX]{x1D6E4}$
increases from 0.001 to 1. In the simulations using the fixed, uniform grid code, which are the cases with relatively less discrete sources (
$\unicode[STIX]{x1D6E4}\geqslant 0.05$
), the leading front propagated through 300 sources, and the average velocities calculated over the last 100, 50 or 10 sources differed by less than 0.1 %. Hence, for the simulations using the fixed, uniform grid code, average velocity over the last 50 sources, i.e. from the 250th to 300th source, is reported in the following subsections of § 4 as
$V_{avg}$
. For the cases with extremely discrete sources (
$\unicode[STIX]{x1D6E4}\leqslant 0.01$
) which were only simulated using the adaptive mesh code, the length of the simulation domain cannot be extended beyond 30 sources due to computational time and memory limits. This length is, however, sufficient for the wave propagation to reach a quasi-periodic state as shown in figure 5(b). Hence, for those cases, it is sufficient to report the average velocity over the last 10 sources as
$V_{avg}$
. All values are normalized by the CJ velocity of the equivalent homogenized media.
4.1 Scenario 1: fixed sources with instantaneous energy deposition
4.1.1 Discreteness
The effect of the spatial concentration of the sources, i.e. the effective discreteness parameter
$\unicode[STIX]{x1D6E4}$
, which is the same as the initial discreteness in this scenario with fixed sources (
$\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{init}=\widetilde{W}_{t=0}/\widetilde{L}$
), on the average shock front velocity (normalized by CJ velocity) is reported in figure 6(a). For these calculations, the value of
$Q$
was held constant at
$50$
, the delay time was a constant
$\unicode[STIX]{x1D70F}=0.01$
and the sources were assumed to remain fixed after the passage of the shock. Note that in all simulations, as the width of the source was expanded to fill the entire space between sources (
$\unicode[STIX]{x1D6E4}\rightarrow 1$
), the average wave velocity converged to the CJ velocity (to within 0.5 %). Also note that in this limit, the energy of each source is still released instantaneously, only now there is no gap between subsequent sources. The fact that the CJ velocity is recovered in the limit of homogeneous energy released in a piecewise continuous manner is good confirmation in the overall implementation of this model.
Simulation results of the average wave velocity normalized by the CJ velocity as a function of (a) effective discreteness, (b) specific heat ratio, (c) heat release and (d) delay time with laboratory-fixed sources. Simulation data obtained using the adaptive refinement code are plotted as open circles. In (a), the simulation results with
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
, 1.4 and
$5/3$
are labelled (i), (ii) and (iii), respectively. In (b), the simulation results with
$\unicode[STIX]{x1D6E4}=0.05$
,
$\unicode[STIX]{x1D70F}=0.01$
and
$Q=50$
are plotted as diamonds and compared to the heuristic model prediction plotted as a solid curve (see the appendix). In (c), the simulation results with
$\unicode[STIX]{x1D6E4}=0.05$
and
$\unicode[STIX]{x1D70F}=0.01$
are plotted as circles and diamonds for
$\unicode[STIX]{x1D6FE}=1.1$
and
$5/3$
, respectively. In (d), the simulation results with
$\unicode[STIX]{x1D6E4}=0.05$
and
$Q=50$
are plotted as circles and diamonds for
$\unicode[STIX]{x1D6FE}=1.1$
and
$5/3$
, respectively.
As the sources are made progressively more discrete by concentrating them in space, an increasing deviation in the average velocity away from the CJ solution is observed. For
$\unicode[STIX]{x1D6FE}=5/3$
, the deviation from CJ reaches a plateau value of
$V_{avg}/V_{CJ}=1.025$
as
$\unicode[STIX]{x1D6E4}$
decreases below
$0.2$
. For lesser values of
$\unicode[STIX]{x1D6FE}$
, greater deviations from CJ are observed and the deviation continues to increase as the sources are made increasingly discrete (
$\unicode[STIX]{x1D6E4}\rightarrow 0$
). For the cases with highly discrete sources, i.e.
$\unicode[STIX]{x1D6E4}<0.05$
, the numerical code with a fixed grid cannot adequately resolve the sources with sufficient resolution (100 computational cells per source) due to computer memory limitations, motivating the use of the adaptive refinement code. The uniform, fixed grid and adaptive mesh refinement codes exhibit good agreement, as seen by comparing the solid and open symbols in figure 6(a), respectively, providing confidence that the results of this study, and the observed super-CJ velocities in particular, are not an artefact of the particular solver used.
For the case of
$\unicode[STIX]{x1D6FE}=1.1$
, the average wave speed achieves a value of 15 % greater than CJ (i.e.
$V_{avg}/V_{CJ}=1.15$
) as the sources become concentrated into 5 % of the available volume (i.e. discreteness
$\unicode[STIX]{x1D6E4}=0.05$
). This is the greatest deviation away from the CJ solution observed in this study. The effect of the ratio of specific heats is explored further in the next section. As the sources are made progressively more discrete (
$\unicode[STIX]{x1D6E4}<0.05$
) in the
$\unicode[STIX]{x1D6FE}=1.1$
case, the average wave speed appears to reach a plateau value. We interpret this plateau as an asymptotic limit in spatial discreteness, and we explore the concept of an asymptotic limit of
$\unicode[STIX]{x1D6FF}$
-function-like sources further in the appendix.
4.1.2 Specific heat ratio
The average wave speed measured for quasi-periodic propagation is reported as a function of the specific heat ratio
$\unicode[STIX]{x1D6FE}$
in figure 6(b). For these simulations, the value of heat release was fixed at
$Q=50$
, the delay time at
$\unicode[STIX]{x1D70F}=0.01$
and the discreteness was fixed at
$\unicode[STIX]{x1D6E4}=0.05$
. From figure 6(a), this value of the discreteness parameter is seen to be sufficiently small that the wave is propagating in the asymptotic limit of highly discretized sources, such that decreasing the value of
$\unicode[STIX]{x1D6E4}$
further will no longer influence the result. The average velocity (as normalized by the CJ velocity) is seen to increase as the value of
$\unicode[STIX]{x1D6FE}$
is decreased. Also shown in this figure is the result of the heuristic model developed in the appendix. The qualitative agreement between the model and the simulations will provide a basis to explain the dependence on
$\unicode[STIX]{x1D6FE}$
, as discussed in § 6.
4.1.3 Heat release
The value of the heat release
$Q$
was varied in a way such that the Mach number of the CJ detonation in the homogenous media remained between values of
$M_{CJ}=4$
and
$M_{CJ}=10$
in order to be representative of real detonable mixtures. This resulted in the value of
$Q$
varying between 50 and 250 for
$\unicode[STIX]{x1D6FE}=1.1$
and between 10 and 60 for
$\unicode[STIX]{x1D6FE}=5/3$
. The average wave speed measurements of these simulations are shown in figure 6(c). The observed wave speed normalized by the CJ speed is not particularly sensitive to the heat release of the mixture.
4.1.4 Delay time
The assumption of a fixed delay time, which is not influenced by the strength of the shock wave that triggers the source or the post-shock thermodynamic states, may be viewed as an unrealistic assumption in the model. Thus, it is important to examine the influence of varying the delay time
$\unicode[STIX]{x1D70F}$
in separate calculations. The effect of varying the delay time, which remained fixed for each simulation, is reported in figure 6(d) for discreteness
$\unicode[STIX]{x1D6E4}=0.05$
and
$Q=50$
. Note that, as the delay time is varied from
$\unicode[STIX]{x1D70F}=0.004$
to
$\unicode[STIX]{x1D70F}=3$
, the observed average wave speed does not exhibit any overall trend. A small deviation is observed around
$\unicode[STIX]{x1D70F}=0.03$
and
$\unicode[STIX]{x1D70F}=0.1$
for
$\unicode[STIX]{x1D6FE}=1.1$
and
$5/3$
, respectively, however, we are unable to attribute any significance to this particular result.
4.2 Scenario 2: convected sources with instantaneous energy deposition
The simulations presented so far are for the case of sources that remain spatially fixed after the passage of the shock. In order to examine the influence of this assumption, a series of simulations were performed in which the energy source is embedded in the medium through which the shock wave propagates, and following the passage of the shock wave, is convected along with the flow. Since a convected source is compressed by the shock wave, the width of a source when it releases energy after a finite delay time
$t_{d}$
is smaller than its initial width before the passage of the shock. Hence, the effective discreteness (
$\unicode[STIX]{x1D6E4}=\widetilde{W}_{t_{i}=t_{d}}/\widetilde{L}$
) is smaller than that in the initial distribution of sources (
$\unicode[STIX]{x1D6E4}_{init}$
), and this effective discreteness is used here to compare these results with those from § 4.1.1 in figure 7(a). The simulations were performed with
$Q=50$
,
$\unicode[STIX]{x1D70F}=0.01$
and
$\unicode[STIX]{x1D6FE}=1.1$
and
$5/3$
. The result of an increasing average wave speed for the cases of convected sources, as the value of
$\unicode[STIX]{x1D6E4}\rightarrow 0$
, coincides with that for the cases of fixed sources. For
$\unicode[STIX]{x1D6E4}$
approximately below 0.1, the average wave speed reaches a plateau value of 15 % and 3 % greater than the CJ value for
$\unicode[STIX]{x1D6FE}=1.1$
and
$5/3$
, respectively.
Simulation results of the average wave velocity normalized by the CJ velocity as (a) a function of effective discreteness
$\unicode[STIX]{x1D6E4}$
for the cases of instantaneously triggered laboratory-fixed (Scenario 1, circles) and convected (Scenario 2, diamonds) sources with
$\unicode[STIX]{x1D70F}=0.01$
,
$Q=50$
and (i)
$\unicode[STIX]{x1D6FE}=1.1$
(solid symbols and trend line) and (ii)
$\unicode[STIX]{x1D6FE}=5/3$
(open symbols and dashed trend line) and (b) a function of initial discreteness
$\unicode[STIX]{x1D6E4}_{init}$
for the cases of
$Q=50$
and
$\unicode[STIX]{x1D6FE}=1.1$
with single-step Arrhenius kinetics (Scenario 3,
$E_{a}=20$
,
$L=10$
, open squares with dashed trend line) and instantaneously triggered, convected sources (Scenario
$2$
,
$\unicode[STIX]{x1D70F}=0.01$
, solid diamonds and trend line).
4.3 Scenario 3: convected sources with single-step Arrhenius kinetics
In order to examine whether the effect of spatial discretization of energy on the wave propagation speed depends on the mechanism of energy deposition, a series of simulations were performed with single-step Arrhenius kinetics instead of the instantaneous energy deposition after a prescribed delay time. In these simulations, the specific heat ratio is selected to be
$\unicode[STIX]{x1D6FE}=1.1$
, which results in the largest deviation from CJ propagation as presented in § 4.1.3. The dimensionless activation energy and pre-exponential factor are chosen to be
$E_{a}=20$
and
$k=637.4$
, respectively, so that the half-reaction-zone length of the homogeneous system, i.e. the distance over which the reaction progress variable
$Z$
evolves from 1 to 0.5 for a steady detonation wave, is unity. Each two consecutive discrete sources are spaced by a distance that is ten times the homogeneous half-reaction-zone length, i.e.
$L=10$
, in these simulations with Arrhenius kinetics. The averaged propagation velocity for the cases of single-step Arrhenius kinetics with the above-mentioned parameters is plotted as a function of the initial spatial discreteness
$\unicode[STIX]{x1D6E4}_{init}$
(open squares) in figure 7(b). Since the reactive layers (in Scenario 3) are continuously reacting once shocked, the measurement of effective discreteness
$\unicode[STIX]{x1D6E4}$
is somewhat arbitrary. Thus, for these simulations, the value of
$\unicode[STIX]{x1D6E4}_{init}$
is used and compared to the results of § 4.2 using
$\unicode[STIX]{x1D6E4}_{init}$
as well. As
$\unicode[STIX]{x1D6E4}_{init}$
decreases, the resulting
$V_{avg}$
gradually deviates from
$V_{CJ}$
and asymptotically approaches a plateau value that is nearly 15 % greater than
$V_{CJ}$
. As shown in figure 7(b), at the plateau super-CJ velocity where
$\unicode[STIX]{x1D6E4}_{init}\leqslant 0.3$
, the result of
$V_{avg}$
for the cases with Arrhenius kinetics coincides with that for the cases with instantaneously triggered, convected sources (plotted as solid diamonds).
5 Analysis
The fact that a spectrum of quasi-periodic solutions that propagate at average speeds significantly greater than the CJ speed (i.e. super-CJ solutions) when the source energy is sufficiently concentrated in space is a result that could be met with a large degree of scepticism given the considerable success of the CJ solution in comparison with experiments. In order to provide an explanation of this finding, the results of select simulations reported in the prior section were analysed using the conventional pressure/specific volume (
$p$
–
$v$
) representation of thermodynamic states. Following the approach developed by White (Reference White1961), Lee & Radulescu (Reference Lee and Radulescu2005), Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007) and Sow, Chinnayya & Hadjadj (Reference Sow, Chinnayya and Hadjadj2014), the flow field of the wave propagating in a discrete source detonation will be analysed in a moving reference frame and the flow at each point in the domain will be averaged over time. This will enable the average structure of the wave to be compared to the classical structure of the Zel’dovich–von Neumann–Döring (ZND) model of a detonation, as visualized in the
$p$
–
$v$
diagram. The analysis shown in this section was only performed for the simulation results for the cases with instantaneously triggered, laboratory-fixed sources.
Similar to the analysis performed by Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007), the averaging is done in a reference frame moving at the average velocity of the wave calculated from the 250th to the 300th source. Thus, in the moving reference frame, the spatial coordinate and particle velocity are transformed as
$x^{\prime }=x-V_{avg}t$
and
$u^{\prime }=u-V_{avg}$
, respectively. For convenience,
$u$
denotes the particle velocity with respect to the moving frame in this section.
A simple time averaging, or Reynolds averaging, procedure is then applied to density and pressure (and other quantities in (5.3)–(5.7)) as follows,
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D70C}}(x^{\prime })=\frac{1}{t_{2}-t_{1}}\int _{t_{1}}^{t_{2}}\unicode[STIX]{x1D70C}(x^{\prime },t)\,\text{d}t, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{p}(x^{\prime })=\frac{1}{t_{2}-t_{1}}\int _{t_{1}}^{t_{2}}p(x^{\prime },t)\,\text{d}t, & \displaystyle\end{eqnarray}$$
where
$t_{1}$
and
$t_{2}$
indicate the starting and ending time of the period, respectively, over which
$\unicode[STIX]{x1D70C}$
and
$p$
are averaged. The overbar ‘
$\overline{~~}$
’ indicates time-averaged variables. Favre averaging (i.e. density-weighted averaging) is applied to particle velocity and specific energy as follows,
$$\begin{eqnarray}u^{\ast }=\frac{\overline{\unicode[STIX]{x1D70C}u}}{\bar{\unicode[STIX]{x1D70C}}}\quad \text{and}\quad u=u^{\ast }+u^{\prime \prime },\end{eqnarray}$$
$$\begin{eqnarray}e^{\ast }=\frac{\overline{\unicode[STIX]{x1D70C}e}}{\bar{\unicode[STIX]{x1D70C}}}\quad \text{and}\quad e=e^{\ast }+e^{\prime \prime },\end{eqnarray}$$
where superscripts
$\ast$
and
$^{\prime \prime }$
indicate Favre-averaged variables and fluctuating quantities, respectively. The average structure of the wave is therefore governed by the one-dimensional, stationary Favre-averaged Euler equations as follows,
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x^{\prime }}\left(\bar{\unicode[STIX]{x1D70C}}u^{\ast }\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x^{\prime }}\left(\bar{\unicode[STIX]{x1D70C}}{u^{\ast }}^{2}+\bar{p}+\overline{\unicode[STIX]{x1D70C}{u^{\prime \prime }}^{2}}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x^{\prime }}\left(\bar{\unicode[STIX]{x1D70C}}e^{\ast }u^{\ast }+\bar{\unicode[STIX]{x1D70C}}\left(e^{\prime \prime }u^{\prime \prime }\right)^{\ast }+\overline{pu}\right)=0. & \displaystyle\end{eqnarray}$$
For an averaged fluid element traversing the wave structure, its averaged mass and momentum are conserved according to,
$$\begin{eqnarray}\overline{M}^{2}=\frac{\bar{p}-1+f}{\unicode[STIX]{x1D6FE}\left(\bar{v}-1\right)},\end{eqnarray}$$
where
$\bar{v}=1/\bar{\unicode[STIX]{x1D70C}}$
is the averaged specific volume,
$\overline{M}=V_{avg}/c_{0}$
is the Mach number of the average wave velocity and
$f=\overline{\unicode[STIX]{x1D70C}{u^{\prime \prime }}^{2}}$
is the Reynolds stress, which measures the intensity of fluctuations in momentum. Note that when
$f=0$
, i.e. there is no fluctuation in fluid momentum over time, equation (5.8) reverts to the ideal Rayleigh line that a fluid element traversing a one-dimensional steady detonation structure is expected to follow.
The average sound speed, which is assumed to be independent of the intensity of fluctuation, can be calculated as,
$$\begin{eqnarray}c^{\ast }=\sqrt{\frac{\unicode[STIX]{x1D6FE}\bar{p}}{\bar{\unicode[STIX]{x1D70C}}}}.\end{eqnarray}$$
The sonic point in the one-dimensional averaged wave structure is located at the position where
$u^{\ast }+c^{\ast }=0$
. The importance of using a mean steady detonation profile to determine the location of the effective sonic plane was highlighted by Lee & Radulescu (Reference Lee and Radulescu2005).
Figure 8 shows the simulation results for the case of continuous (i.e.
$\unicode[STIX]{x1D6E4}=1$
) laboratory-fixed energy sources,
$Q=50$
,
$\unicode[STIX]{x1D70F}=0.1$
and
$\unicode[STIX]{x1D6FE}=1.1$
analysed using the Favre-averaging approach and
$p$
–
$v$
representation of thermodynamic states. Note that while the source energy is uniform throughout the medium, it is still released in piecewise-continuous segments with width
$\widetilde{L}$
, which is non-dimensionalized to a value of unity when reported here. In figure 8(a), the evolution of the flow field is plotted as a grey scale contour of pressure in an
$x^{\prime }$
–
$t$
diagram. Figure 8(b) shows the spatial profile of
$f$
in the wave-attached reference frame, which measures the intensity of mechanical fluctuations around the Favre-averaged value in the flow field, normalized by
$Q$
. The one-dimensional steady structure of the time-averaged pressure
$\bar{p}$
in the wave-attached reference frame is plotted in figure 8(c) below the
$x^{\prime }$
–
$t$
diagram. The thermodynamic path taken by a fluid element traversing this one-dimensional steady, averaged wave structure is plotted in a
$p$
–
$v$
diagram as shown in figure 8(d). The location of the sharp increase in
$\bar{p}$
from the upstream initial state to its peak value occurs in the region between the leading shock and the loci of the source releasing energy in the
$x^{\prime }$
–
$t$
diagram. The peak in the averaged pressure matches the von Neumann pressure, i.e. the downstream pressure behind a normal shock of strength
$M_{CJ}$
, shown in figure 8(d). As shown in figure 8(a), the
$u+c$
characteristics near the separatrix are quasi-straight lines with slight wiggles as they pass through backward-running shock waves and contact surfaces.
For the case of continuous (i.e.
$\unicode[STIX]{x1D6E4}=1$
) laboratory-fixed energy sources,
$Q=50$
,
$\unicode[STIX]{x1D70F}=0.1$
and
$\unicode[STIX]{x1D6FE}=1.1$
, (a) pressure contour plotted as an
$x^{\prime }$
–
$t$
diagram showing the evolution of the flow field and right-running characteristics, (b) the spatial profile of the intensity of mechanical fluctuation
$f$
normalized by
$Q$
, (c) the one-dimensional steady structure of the time-averaged pressure
$\bar{p}$
in the wave-attached reference frame and (d) the
$p$
–
$v$
representation of the one-dimensional averaged wave structure. In (d), curve (i) is the equilibrium Hugoniot curve of complete energy release
$Q=50$
, curve (ii) is the Hugoniot curve of zero energy release
$Q=0$
, curve (iii) is the ideal Rayleigh line corresponding to a leading shock strength of
$M_{CJ}$
and curve (iv) is the
$p$
–
$v$
thermodynamic path taken by a fluid element traversing the one-dimensional averaged wave structure.
Comparing figures 8(a) and 8(c), the sonic point in the averaged profile is found to match the trajectory of the limiting characteristic of the separatrix, which appears downstream from the loci of sources releasing energy. The thermodynamic state at this averaged sonic point nearly coincides with the CJ state of a one-dimensional steady detonation wave propagating at
$M_{CJ}$
, as seen in figure 8(d). This is consistent with the result that the
$\overline{M}$
obtained from the simulation with continuous energy sources is only 0.3 % greater than the
$M_{CJ}$
for the same total energy release
$Q$
. Since the CJ point is the intersection of the ideal Rayleigh line and the equilibrium Hugoniot curve of complete energy release, the averaged flow field relaxes to equilibrium at the sonic point. (It is important to note that, in this problem, the term ‘equilibrium’ refers to mechanical (or hydrodynamic) equilibrium that the flow reaches at the sonic point in the classical CJ solution, rather than chemical equilibrium of a system with reversible reactions.) In figure 8(b), it can be seen that the mechanical fluctuation
$f/Q$
is significantly positive only in a narrow region which corresponds to the sharp increase in averaged pressure shown in figure 8(c). In figure 8(d), the averaged
$p$
–
$v$
path from the upstream initial state to the peak averaged pressure (or the von Neumann) state significantly departs from the ideal Rayleigh line corresponding to a wave speed of
$M_{CJ}$
, which can be interpreted as the fact that the mechanical fluctuation is concentrated between the shock and energy release, as shown in figure 8(b). As
$f/Q$
rapidly decays to and remains at zero downstream of the location of peak averaged pressure to the sonic point, the averaged
$p$
–
$v$
path connecting the peak averaged pressure state to the sonic point coincides with the ideal Rayleigh line. In the case of continuous energy sources, despite the mechanical fluctuation caused by the interaction between the forward-running blast generated by each individual source releasing its energy and the preceding shock front, the highly unsteady wave propagation and energy deposition process has a temporal average structure that is nearly equivalent to a steady one-dimensional detonation structure with the same amount of heat release. The excellent agreement obtained between the averaged structure of the wave extracted from the simulations (with continuous energy release) and the steady ZND structure can be taken as validation of the averaging method used in this study.
Figure 9 shows the simulation results for the case of highly discrete (
$\unicode[STIX]{x1D6E4}=0.05$
) laboratory-fixed energy sources,
$Q=50$
,
$\unicode[STIX]{x1D70F}=0.1$
and
$\unicode[STIX]{x1D6FE}=1.1$
analysed using the Favre-averaging approach and the
$p$
–
$v$
representation of thermodynamic states. As seen from figure 6(b),
$\unicode[STIX]{x1D6E4}=0.05$
approximately corresponds to the limit of discreteness where the wave velocity becomes independent of the value of
$\unicode[STIX]{x1D6E4}$
(i.e. making the sources even more discrete does not result in greater deviation from the CJ solution). The
$x^{\prime }$
–
$t$
diagram of flow field evolution, the spatial profile of
$f/Q$
, the one-dimensional steady profile of
$\bar{p}$
and the
$p$
–
$v$
diagram of the averaged wave structure are plotted in figures 9(a), 9(b), 9(c) and 9(d), respectively. Comparing figures 9(a) and 9(c), it can be found that the increase in averaged pressure from the initial state to its peak value spans over the region where the forward-running blast waves interact with the leading shock front in the
$x^{\prime }$
–
$t$
diagram. As shown in figure 9(b), the maximum intensity of mechanical fluctuation is associated with wave interactions in this region.
For the case of highly discrete (
$\unicode[STIX]{x1D6E4}=0.05$
) laboratory-fixed energy sources,
$Q=50$
,
$\unicode[STIX]{x1D70F}=0.1$
and
$\unicode[STIX]{x1D6FE}=1.1$
, (a) pressure contour plotted as an
$x^{\prime }$
–
$t$
diagram showing the evolution of the flow field and right-running characteristics, (b) the spatial profile of the intensity of mechanical fluctuation
$f$
normalized by
$Q$
, (c) the one-dimensional steady structure of the time-averaged pressure
$\bar{p}$
in the wave-attached reference frame and (d) the
$p$
–
$v$
representation of the one-dimensional averaged wave structure. In (d), curve (i) is the equilibrium Hugoniot curve of complete energy release
$Q=50$
, curve (ii) is the Hugoniot curve of zero energy release
$Q=0$
, curves (iii) and (iv) are the ideal Rayleigh lines corresponding to a leading shock strength of
$M_{CJ}$
and
$\overline{M}$
, respectively, and curve (v) is the
$p$
–
$v$
thermodynamic path taken by a fluid element traversing the one-dimensional averaged wave structure.
In contrast to the continuous case (figure 8), as shown in figure 9(a), the
$u+c$
characteristics exhibit a ‘zig-zag’ pattern as they are processed by much stronger backward-running blast waves. The trajectory of the limiting characteristic (i.e. separatrix) oscillates back and forth while its overall zig-zag pattern remains in a fixed region relative to the leading wave front over time. For this particular calculation, the sonic point found in the one-dimensional Favre-averaged wave structure coincides with the right (or upstream) end of the limiting characteristic. Similar to the continuous case, the location of the sonic point is significantly downstream from the loci of discrete sources releasing energy. For instantaneous energy deposition, this result, wherein the forward-most point of characteristic propagation corresponds to the sonic point found in the averaging procedure, was verified for all cases considered.
Plotting the
$p$
–
$v$
path taken by a fluid element traversing the Favre-averaged wave structure in figure 9(d), it can be seen that the peak averaged pressure is greater than the von Neumann pressure behind a shock of strength
$M_{CJ}$
. The averaged sonic point departs from the CJ point, and is not found on the equilibrium Hugoniot curve for
$Q=50$
. Hence, equilibrium is not reached as the separatrix and the effective sonic surface in the flow field are encountered. Equilibrium is eventually attained in the far field, downstream from the sonic point, as this downstream state coincides with the intersection between the equilibrium Hugoniot curve of
$Q=50$
and the Rayleigh line of a wave speed
$\overline{M}$
. The averaged pressure at this equilibrium state is significantly less than that of the CJ state. The
$\overline{M}$
for
$\unicode[STIX]{x1D6E4}=0.05$
discrete case, which is 14 % greater than
$M_{CJ}$
, is associated with a weak detonation solution. As shown in figure 9(b), the mechanical fluctuation retains a significantly positive value from the leading wave front to a location beyond the sonic point, such that the entire averaged
$p$
–
$v$
path departs from the ideal Rayleigh line based on a one-dimensional steady detonation structure. The quantity
$f/Q$
eventually decays to zero in the far-field downstream, consistent with the result that the far-field downstream thermodynamic state reaches the intersection between the equilibrium product Hugoniot curve and the Rayleigh line associated with the wave speed
$\overline{M}$
.
6 Discussion
The results presented in this paper demonstrate that the classical CJ solution for detonation velocity breaks down and is no longer valid as the source of energy is spatially concentrated into highly discretized sources. In order to explain this result, the analysis in § 5 shows that flow crossing the separatrix (i.e. the limiting characteristic that can influence the shock front) has not reached equilibrium due to on-going shock wave interactions that are occurring in the flow for the case of a highly discretized medium. This separatrix is the unsteady flow analogue to a sonic surface, and the fact that the flow is not in equilibrium at the effective sonic surface is a hypothesis as to why non-CJ speeds are observed. As seen in figure 9(d), the flow does eventually reach the equilibrium Hugoniot downstream of the effective sonic surface. This result suggests that the super-CJ speeds in this study are properly designated as weak detonations. Prior examples of weak detonations include so-called pathological detonations, in which competing exothermic and endothermic reactions permit the flow in the reaction zone to pass through a sonic point while in a non-equilibrium state (von Neumann Reference von Neumann1942; Fickett & Davis Reference Fickett and Davis2000; Higgins Reference Higgins and Zhang2012), and ultrafast detonations (Anderson & Long Reference Anderson and Long2003), in which chemical reactions are sufficiently fast that significant reaction occurs within the shock front itself, enabling the wave to bypass the von Neumann state and proceed directly to the weak branch of the product Hugoniot. The mechanism suggested in this paper, namely, a high degree of flow non-equilibrium resulting from shock interactions, is believed to be a new mechanism able to realize weak (non-CJ) solutions.
The deviations from CJ speed depend significantly upon the ratio of specific heats,
$\unicode[STIX]{x1D6FE}$
, as seen in figure 6(b). The heuristic model developed in the appendix, which examines the limiting case of the time delay going to zero (
$\unicode[STIX]{x1D70F}\rightarrow 0$
) and the spatial discreteness also approaching zero (
$\unicode[STIX]{x1D6E4}\rightarrow 0$
), such that the energy release of the first source will result in a classical point blast solution, may provide an explanation for this trend. In this model, subsequent sources will release their energy on the shock front of the prior source, and as discussed in the appendix, the partitioning of energy release of that source into forward and backward directed blast waves is related to the density ratio in the limit of a strong shock front (
$(\unicode[STIX]{x1D6FE}+1)/(\unicode[STIX]{x1D6FE}-1)$
). Since this function increases rapidly as the value of
$\unicode[STIX]{x1D6FE}$
approaches unity (
$\unicode[STIX]{x1D6FE}\rightarrow 1$
), we would expect the energy release of subsequent sources to be increasingly directed in the direction of propagation. In simpler language, in the limit of
$\unicode[STIX]{x1D6FE}$
going to unity, the strong shock front becomes a solid wall, so that energy released on that wall results in a blast wave directed entirely forward. This picture may provide a qualitative explanation for the influence of
$\unicode[STIX]{x1D6FE}$
, and indeed the simple heuristic model developed in the appendix performs better in predicting the average wave speed than might be expected, as seen in figure 6(b).
As the value of dimensionless heat release
$Q$
varies in a way such that the CJ detonation velocity for an equivalent homogeneous reactive medium remains in the range between
$M_{CJ}=4$
and
$M_{CJ}=10$
, the ratio of the resulting
$V_{avg}$
in the discrete reactive system to the homogenous equivalent
$V_{CJ}$
exhibits no significant change (see figure 6(c)). An insight obtained from the heuristic model is that the key dynamics of the wave propagation in a medium of discrete sources can be fairly well pictured as a complex superposition of blast waves generated by a series of point sources. According to the Taylor–Sedov blast wave solution, the average speed at which the blast wave front propagates over a fixed distance has a square root dependence on the energy released by a planar source, namely,
$\sqrt{Q}$
. In the limit of large heat release, the CJ velocity also linearly depends on
$\sqrt{Q}$
. Hence, this dependence of
$\sqrt{Q}$
might be expected to cancel out when
$V_{avg}$
is normalized by
$V_{CJ}$
, and this may explain the results presented here.
The result that the deviations of the average wave speed away from the CJ speed do not depend on delay time is another result from this study, as seen in figure 6(d). The delay time was varied between
$\unicode[STIX]{x1D70F}=0.004$
, meaning that the sources release their energy effectively instantaneously upon contact with the leading shock, and
$\unicode[STIX]{x1D70F}=3$
, meaning that the blast wave released by a source must approximately propagate a distance of three source spacings before it reaches the location where the shock was when the source released its energy, by which time the leading shock had moved further downstream as well, resulting in it taking even longer for the blast to reach the leading shock. In this limit of large
$\unicode[STIX]{x1D70F}$
, given the large number of blast wave interactions that occur before the blast from a source actually reaches the leading shock and contributes to its sustenance, it is perhaps remarkable that the average speed does not depend upon the delay time. Note, however, that in the classical CJ detonation model, the details of the reaction zone, such as kinetic rates, do not affect the detonation wave velocity. In this connection, it is perhaps plausible that the results do not depend on the delay time of the discrete sources.
As shown in figure 7(b), in the limit of highly discretized sources (
$\unicode[STIX]{x1D6E4}\rightarrow 0$
), the averaged wave speed resulting from the cases with single-step Arrhenius kinetics, where chemical energy of the sources is released at a state-dependent finite rate and reaches the same plateau super-CJ value, i.e.
$V_{avg}/V_{CJ}=1.15$
for
$\unicode[STIX]{x1D6FE}=1.1$
, as that resulting from the cases with an instantaneous addition of energy after a state-independent delay time. This result suggests that the effect of energy discretization on the wave propagation speed is very likely independent of the detailed mechanism of how energy is deposited into the flow. The mechanism of instantaneous energy deposition used in most of this study is likely to play a role similar to point-blast theory (Taylor–Sedov blast), which defines the limiting case when the blast source is sufficiently concentrated in space and time.
The results of this study may have application to resolving the conundrum that was discussed in the Introduction (§ 1), namely, resolving why the CJ criterion is so successful in describing the average propagation velocity of transient and multi-dimensional detonation waves. Unstable detonations feature pockets of gas that may become compressed by multiple shock waves and explode, while other portions of the detonation cell might be described by a non-reacting shock wave. The study of Kiyanda & Higgins (Reference Kiyanda and Higgins2013) estimated that nearly the entire second half of a detonation cell in low pressure methane/oxygen does not react coupled to the shock front, but rather burns as shock-compressed pockets detached from the leading front. The effective discreteness of such unstable detonations, as defined in this paper, is unlikely to be less than
$\unicode[STIX]{x1D6E4}=0.3$
. As shown in figure 7, it is unlikely that deviations from the ideal CJ speed would be more than a few per cent for this degree of discreteness, a difference that is hardly experimentally resolvable. This result may also have relevance to the phenomenon of galloping detonation, wherein steady detonation propagation is not possible due to the tube diameter being smaller than a characteristic cell size required for propagation, such that the detonation fails to approximately half CJ speed and periodically reinitiates to an initially overdriven wave via a process similar to deflagration to detonation transition. Despite the fact that the cycle of the galloping detonation occurs over hundreds of tube diameters and the wave should experience significant heat and momentum losses due to the relatively small size of the tube used, galloping detonations continue to propagate remarkably close to the CJ speed on average (Lee et al.
Reference Lee, Dupré, Knystautas and Lee1995; Gao, Ng & Lee Reference Gao, Ng and Lee2015; Jackson, Lee & Shepherd Reference Jackson, Lee and Shepherd2015). This aspect of the discrete source detonation model has been further explored by Radulescu & Shepherd (Reference Radulescu and Shepherd2015). The possibility that all real detonations might be weak detonations (if only slightly away from the CJ solution) is intriguing, as this idea was suggested by Davis (Reference Davis1987) but has not attracted significant attention in detonation research.
The results of the present study are fundamentally different from those of Mi & Higgins (Reference Mi and Higgins2015), which examined the discrete source detonation problem in a Burgers analogue system. The one-dimensional, inviscid scalar Burgers equation with a source term was proposed by Fickett & Davis (Reference Fickett and Davis2000) and Majda (Reference Majda1981) as an analogue to the reactive Euler equations in order to explore detonation dynamics. In recent years, study of this scalar analogue system has generated a number of interesting results, including the existence of pulsating, chaotic solutions (Radulescu & Tang Reference Radulescu and Tang2011; Kasimov, Faria & Rosales Reference Kasimov, Faria and Rosales2013; Tang & Radulescu Reference Tang and Radulescu2013; Faria & Kasimov Reference Faria and Kasimov2015). The study of Mi & Higgins (Reference Mi and Higgins2015) was a precursor to the present study using the inviscid Burgers equation with periodically spaced
$\unicode[STIX]{x1D6FF}$
-function sources that were triggered by the passage of the leading shock front. Both regularly spaced sources with fixed delay and randomly spaced sources with randomly generated delays were considered. The resulting wave dynamics involved the interaction of a number of sawtooth-profiled blast waves generated by the sources. For the Burgers equation the trajectories of the blast waves could be solved analytically. In all cases considered, the average propagation velocity was found to be within 1 % of the CJ detonation velocity of the equivalent CJ system. That result is in seeming contradiction to the results of the present study, which exhibited super-CJ solutions, so comparing these two studies warrants further discussion. The Euler equations have three families of propagating characteristics, namely, right-running, left-running and particle characteristics, while the Burgers equation has only the family of right-running characteristics. This feature of the Burgers equation system might explain the unusual dynamics noted in Mi & Higgins (Reference Mi and Higgins2015) (see § 4.2 and figure 7(b)) wherein the limiting characteristic became coincident with the locus of new sources in the case of regularly spaced sources with a fixed delay, while in the present study, the limiting characteristic is spatially separated from the locus of new sources (see figures 4 and 9 in the present paper). This aspect of the Euler equations means that the dynamics of the shock front influences the post-shock flow, which in turn influences the shock again, resulting in the entire flow being a non-simple wave region. In the present study, the super-CJ velocity in a discrete reactive medium is a direct result of the significantly intense mechanical fluctuation near the limiting characteristic (as explained in § 5 and shown in figure 9), which is caused by the downstream-propagating blast waves generated by the sources releasing energy. In contrast, due to the lack of left-running characteristics, a source releasing energy cannot generate a blast wave influencing the flow field downstream from the source in the Burgers equation system (Mi & Higgins Reference Mi and Higgins2015). Thus, we can hypothesize that it is the existence of both right- and left-running characteristics that is the necessary feature to observe the super-CJ solutions found in this study.
Extension of the current, Euler equation-based study to higher dimensions is likely to reveal richer behaviour. Examination of the propagation of flames in three-dimensional systems of random, discrete sources has revealed a percolation-like regime that has demonstrated the ability of flames to propagate beyond the thermodynamic limit of the corresponding homogeneous medium (Tang et al. Reference Tang, Higgins and Goroshin2012), and similar behaviour may occur with detonations in discrete systems. Examination of detonation propagation in the limit of spatially randomized, point-like sources in three-dimensional clouds may have some relevance to the anomalous scaling between experiments in axisymmetric geometries (tubes) and two-dimensional slab geometries (rectangular geometries) (Gois, Campos & Mendes Reference Gois, Campos and Mendes1996; Petel et al. Reference Petel, Mack, Higgins, Turcotte and Chan2006, Reference Petel, Mack, Higgins, Turcotte and Chan2007; Higgins Reference Higgins2009, Reference Higgins2014). A preliminary step in this direction was made in the recent computational study of Li et al. (Reference Li, Mi and Higgins2015), in which a two-dimensional sinusoidal ripple in properties was introduced in a medium prior to detonation propagation through it. Extension of this study using perturbations with greater discreteness than sinusoidal would be of great interest.
Detonation theory has profitably explored the asymptotic limits of high activation energy (Buckmaster & Neves Reference Buckmaster and Neves1988; Short Reference Short1997), zero activation energy (Erpenbeck Reference Erpenbeck1965), high overdrive (Clavin & He Reference Clavin and He1996; Clavin, He & Williams Reference Clavin, He and Williams1997), low energy release (Short & Stewart Reference Short and Stewart1999) and the Newtonian limit (i.e. ratio of specific heats approaching unity) (Short Reference Short1996). The line of investigation suggested in this study examines a different type of asymptotic limit, the limit of spatially discrete energy sources approaching
$\unicode[STIX]{x1D6FF}$
-functions in space and time. It is unknown if a rational asymptotic solution to the problem, rigorously derived from the governing Euler equations can be found, in contrast to the ad hoc solution constructed in the appendix. This problem is left as an open question to the detonation theory community.
7 Conclusions
Detonation propagation in systems with the energy release of the medium spatially concentrated into discrete pockets was simulated computationally. The energy release of one such discrete source drives a blast wave which is capable of initiating the next source after a prescribed delay. The resulting ensemble of blast wave interactions propagates in a quasi-periodic manner, and the average wave speed was measured. Systematic variation of the model parameters identified that the average wave deviated significantly from the CJ detonation speed of the equivalent homogenous media, with speeds as great as 15 % in excess of the CJ speed being measured. This discrepancy is significant, given that experimental measurements of detonation speeds in gases usually agree to within 1 % of the equilibrium CJ speed. A systematic variation of model parameters found that the deviation away from CJ depended on the degree of spatial concentration of the sources (more discrete sources resulted in greater deviations above CJ), with the wave speed reaching an asymptotic limit as the sources were concentrated into a space occupying less than 1 % of the entire domain (limit of
$\unicode[STIX]{x1D6E4}\rightarrow 0$
). The deviation from CJ also depended on the ratio of specific heats (with greater deviations observed as
$\unicode[STIX]{x1D6FE}\rightarrow 1$
), but the deviation from CJ did not depend upon the delay time of the sources between being shocked and releasing their energy or the average value of the energy release. The results were interpreted via temporal averaging of the simulations onto a steady, one-dimensional projection that could be compared with the classical ZND structure of detonations. This analysis suggests that the existence of non-CJ solutions can be interpreted as weak detonations due to the non-equilibrium of the flow resulting from the ensemble of shock interactions as the flow crosses the effective sonic surface.
Acknowledgements
The authors would like to acknowledge extended discussions with S. Goroshin, J. Lee, M. Radulescu, V. Tanguay, H. D. Ng and C. Kiyanda that encouraged this study. Earlier simulations of the discrete detonation problem (the results of which, while not directly reported in this study, are consistent with the findings) were performed by A.J.H. using Clawpack (‘Conservation Laws Package’) by R. Leveque. Additional simulations by M. Javaid performed with E.V.T. provided further impetus and validation for this study.
Appendix A. Heuristic model
The collection of a large number of interacting blast waves and other unsteady flow features that comprise the dynamics of the discrete source detonation problem are unlikely to be amenable to an analytic solution. In the limit of point-like energy sources (discreteness
$\unicode[STIX]{x1D6E4}\rightarrow 0$
) and the delay time going to zero (
$\unicode[STIX]{x1D70F}\rightarrow 0$
), however, the problem becomes simpler. In this case, as shown in figure 10(a), the blast wave from a new source occurs on the shock front of the prior blast, initially with no other flow interactions involved. This problem, and the subsequent propagation of the blast, are considered in this appendix in an attempt to construct an analytic solution to this problem. We emphasize that this is a constructed solution, in which prior known solutions are patched together in an ad hoc, or heuristic, manner, rather than a solution to the discrete detonation problem derived rigorously from the governing conservation laws.
An approximate, analytical solution to this problem for the case where the delay time and discreteness are both taken to the limit of zero begins with the similarity solution for the planar version of the well-known point-blast problem of Taylor and Sedov (Jones Reference Jones1961). Treating a detonation wave as an interaction between blast waves draws some inspiration from the heuristic model of Vasil’ev & Nikolaev (Reference Vasil’ev and Nikolaev1978), which utilized the point-source blast wave solution to treat detonation cell structure. In this current model, the motion of the shock front
$x_{s}$
is given by the following non-dimensionalized equation,
$$\begin{eqnarray}\frac{\text{d}x_{s}}{\text{d}t}=\sqrt{\frac{Q}{Bx_{s}}},\end{eqnarray}$$
where
$B$
is a dimensionless energy parameter depending on the specific heat capacity ratio
$\unicode[STIX]{x1D6FE}$
.
Schematic illustration of the one-dimensional, discrete-source detonation propagation problem, as an
$x$
–
$t$
(space–time) diagram (a). In (b), the mechanism of source energy deposition is assumed to occur via the impulsive motion of two outward-facing, massless pistons. The approximation of the carried-over influence from the previous sources onto the
$(i-1)$
th source is shown in (c). In (d), the carried-over particle velocity imposed on the
$i$
th source,
$V_{p,0}^{i}$
, is approximated as the piston velocity required to sustain the instantaneous shock velocity,
$V_{s,0}^{i}$
.
Since the source releases its energy at the shock front originating from the previous source, the partitioning of blast energy into forward- and backward-propagating blast waves can be estimated as follows. As shown in figure 10(b), if the mechanism of source energy deposition is assumed to be two hypothetical massless pistons, one that pushes outward into the undisturbed gas ahead of the blast and the other that pushes into the gas behind the blast from the prior source, the pressure on both piston faces must be equal (since, being massless, they can exert no net force on the flow). This condition permits the ratio of piston velocities to be solved for and, since the pistons must act for the same duration of time to conserve momentum, this ratio also determines the ratio of work done. Via the Rankine–Hugoniot normal shock relations, the ratio of particle velocity in front of each piston (i.e. piston velocity) and the work done by each piston can be obtained, and thus the energy partitioned to contribute to the forward- and backward-running blast waves can be determined. In the limit of strong shock relations, the equal pressure condition permits the partitioning of energy into forward- and backward-propagating blasts to be expressed analytically as a function of
$\unicode[STIX]{x1D6FE}$
,
$$\begin{eqnarray}\frac{Q_{forward}}{Q}=\unicode[STIX]{x1D702}=\frac{1}{\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}+1}}}+1}.\end{eqnarray}$$
Since only
$\unicode[STIX]{x1D702}Q$
contributes to the forward-propagating blast, the equation describing the motion of the leading shock front needs to be modified,
$$\begin{eqnarray}\frac{\text{d}x_{s}}{\text{d}t}=\sqrt{\frac{\unicode[STIX]{x1D702}Q}{Bx_{s}}}.\end{eqnarray}$$
Interestingly, equation (A 2) is the same partition of blast energy found by Sakurai (Reference Sakurai1974) in examining planar blast waves generated by energy release at a stationary density interface (as opposed to energy release at a moving shock front in the current problem).
Following the release of source energy, the subsequent blast wave motion is additionally influenced by the particle velocity that was imposed by the blast wave from the previous source. This problem (in position–time space) is analogous via hypersonic similarity to the two-dimensional hypersonic flow past a blunted wedge, a problem previously treated by Chernyi (Reference Chernyi1961, pp. 218-221). The blunt leading edge of the wedge represents the instantaneous energy release of the source under the hypersonic blast wave analogy, and the surface of the wedge represents the particle motion imposed by the previous blast wave (as determined by the blast strength from the previous source as it reached the new source). Considering the energy partitioned into the forward-propagating blast and the particle motion imposed by the previous source, the motion of the blast wave propagation from one source to the next (
$x_{s}$
from
$0$
to
$1$
) is given by,
$$\begin{eqnarray}\frac{\text{d}x_{s}}{\text{d}t}=\sqrt{\frac{\unicode[STIX]{x1D702}Q}{Bx_{s}}}+V_{p}\left(x_{s}\right),\end{eqnarray}$$
where
$V_{p}$
is the particle velocity carried over from the previous source.
The carried-over particle velocity at the location of a newly triggered (the
$i$
th) source, i.e.
$V_{p}^{i}(x_{s}=0)$
or
$V_{p,0}^{i}$
, can be approximated as follows. The blast wave, which imposes
$V_{p,0}^{i}$
onto the
$i$
th source, is a result of the energy release of not only the
$(i-1)$
th source, but also all the previous sources (from the first one) with diminishing effect. The influence of the previous sources on the leading blast propagation after the trigger of the
$(i-1)$
th source is, however, too complicated to be solved exactly. As shown in figure 10(c), this influence is approximated as the kinetic energy possessed by the particle motion imposed on the
$(i-1)$
th source, i.e.
$(V_{p,0}^{i-1})^{2}/2$
, added to the forward-partitioned source energy
$\unicode[STIX]{x1D702}Q$
to propel the blast propagating towards the
$i$
th source. With this carried-over amount of energy from the previous sources, the velocity, at which the blast reaches the
$i$
th source, can be obtained from the Taylor–Sedov solution,
$$\begin{eqnarray}V_{s,0}^{i}=\sqrt{\frac{\unicode[STIX]{x1D702}Q+(V_{p,0}^{i-1})^{2}/2}{B}}.\end{eqnarray}$$
The particle velocity imposed on the
$i$
th source,
$V_{p,0}^{i}$
, as shown in figure 10(d), can be then approximated as the piston velocity required to sustain the blast front moving at its approximated instantaneous velocity
$V_{s,0}^{i}$
,
$$\begin{eqnarray}V_{p,0}^{i}=\frac{2}{\unicode[STIX]{x1D6FE}+1}V_{s,0}^{i}=\frac{2}{\unicode[STIX]{x1D6FE}+1}\sqrt{\frac{\unicode[STIX]{x1D702}Q+(V_{p,0}^{i-1})^{2}/2}{B}}.\end{eqnarray}$$
Note that (A 6) is a recursion relation, similar to the concept suggested by Stewart & Asay (Reference Stewart, Asay, Fife, Liñán and Williams1991). As the carried-over particle velocity is zero when the first source releases energy, the subsequent values of
$V_{p,0}$
can be calculated recursively. As the blast wave propagates to the next source, the influence of the previous source on the particle motion should diminish. In order to consider this effect,
$V_{p}$
is modelled to be inversely proportional to
$x_{s}$
,
$$\begin{eqnarray}V_{p}(x_{s})=\frac{1}{x_{s}+1}V_{p,0}.\end{eqnarray}$$
Note that the ‘1’ in (A 7) is the spacing between two adjacent sources in the non-dimensionalized coordinates. This simple approximation can be justified by the fact that the particle velocity profile of a blast wave is roughly linear. Thus, the approximate, analytic solution of the leading shock wave propagation from one discrete source to the next can be obtained by integrating (A 4) with (A 6) and (A 7).
The predictions of this model for the instantaneous velocity of the shock front are shown in figure 5(a). Note that, since the Taylor–Sedov similarity solution for a point blast (in planar geometry) is used, the shock velocity at each source is initially infinite. The velocity decays as the blast propagates forward, then jumps up again as a new source is triggered. The minimum velocity before each new source is triggered monotonically increases as more sources are encountered, in qualitative agreement with the computational simulations. In figure 5(b), the average velocity (from one source to the next) is plotted as a dashed line, and again exhibits good qualitative agreement with simulations in the limit as discreteness
$\unicode[STIX]{x1D6E4}\rightarrow 0$
, although it tends to over-predict the final value of the plateau velocity. This model can be used to examine the effect of
$\unicode[STIX]{x1D6FE}$
, which enters the model significantly via the partitioning of blast energy released at the shock front (A 2), as shown in figure 6(b). Again, the model captures the qualitative trend of an increasing deviation away from the CJ solution as
$\unicode[STIX]{x1D6FE}$
approaches unity, although the predicted deviation is more than twice that observed in the simulations.
