Numerical modelling of steady detonations with a variational streamline approach

Previous studies (Watt et al. in J Eng Math 75(1):1, 2012; Cartwright and Falle in J Eng Math 115(1):157, 2019) have shown that a streamline based approach to modelling of steady state detonations can produce good results for rate laws which have maximal reaction at the shock. In this paper we consider a Variational Streamline Approximation (VSA) which introduces streamline curvature. Comparing results with Direct Numerical Simulations (DNS) and the existing Straight Streamline Approximation (SSA) model, we find that the VSA improves the predictive accurary of streamlines modelling compared to DNS calculations, capturing the shock front and sonic surfaces with greater accuracy than SSA.

Schematic of a steady-state detonation. In the frame of the shock the unreacted explosive flows into the shock wave in the negative y-direction. The streamlines intersect the shock at the point (x f , y f ). The flow speed in the shock frame increases as the chemical reaction proceeds and goes sonic at the sonic locus Due to the high cost of performing experiments, computational modelling provides an alternative to study the behaviour of explosives in complex geometries and can be used to characterise existing and novel explosives. The most common approach to modelling explosives is based upon the assumption that we can model them at the continuum level, with the reactive Euler equations describing the physics. The explosive material, and its reactants, are described by an equation of state (EOS) and a chemical reaction rate model [1]. Popular choices of EOS include Jones-Wilkins-Lee, Williamsburg, polytropic and finite-strain. The equations describing the chemical reactions are coupled to the shock strength by making the reaction rates depend upon the local thermodynamic variables. The Ignition & Growth Model [2] has reaction terms which are a function of the pressure. The CREST reactive burn model [3] has reaction terms which are function of the entropy of the shocked unreacted explosive. The EOS and reaction rate laws are empirical in nature and require calibration to experimental data for a given explosive [4][5][6].
The most accurate method of solving the problem is by using high resolution direct numerical simulations (DNS) [7] where the full time dependent equations are solved. For high resolution calculations, DNS is computationally expensive. When calibrating explosive models where hundreds, or maybe thousands, of simulations are required, a computationally cheap model is desirable. To this end, various simplified detonation models have been developed, which are significantly more efficient than DNS. The first such model was the Wood-Kirkwood model [8], but calibration to experimental data is required to obtain the relationship between the shock front curvature and the charge radius [5]. An improvement to the Wood-Kirkwood model came with development of the Detonation Shock Dynamics (DSD) model [9,10] which is based upon the assumption that the detonation velocity is close to the Chapman-Jouguet (CJ) velocity and that the shock front curvature is small. Although reasonably accurate for ideal explosives where the detonation velocity does not depart significantly from the CJ velocity, it tends to give too low detonation velocities for non-ideal explosives where detonation velocities can be significantly lower than CJ velocity. Watt et al. [11] developed a more sophisticated model, the Straight Streamline Approximation (SSA). This makes no assumptions about the shock front curvature, instead it assumes that the post-shock streamlines are straight. They showed that it can accurately reproduce the velocity-radius curves from high resolution DNS calculations for non-ideal explosives. A follow up study found that for detonation models whose reaction rate models did not have maximal reaction at the shock, the SSA proved to be less accurate [12].
In this paper we extend the streamline detonation model, relaxing the assumption that the streamlines are straight and allowing the streamlines to curve. We investigate a simple explosive model with the power-law reaction and a reaction rate with an induction zone.

Reactive Euler equations
The problem is described by the compressible Euler equations with reaction terms: where is the material derivative. Here ρ is the density, u the velocity, p the pressure, e the internal energy per unit mass and λ the chemical reaction progress variable (ranging from 0 for no reaction to 1 for complete reaction). W is the chemical reaction rate, which depends upon λ and local thermodynamic variables.

Streamline detonation model
The first step in deriving the equations for the streamline based model is to define a compressible streamfunction, ψ, which in slab geometry is where u and v are the velocities in the cartesian x and y directions respectively. Equations (3) and (4) give For steady detonations, the material derivative, Eq. (2), is The reactive Euler equations are transformed from (x, y) coordinates to streamline based coordinates (ψ, y) via application of Eqs. (3)-(6) to give

Master equation in v
The governing equations in streamline coordinates, (7)-(10) can be reduced to a master equation in the y velocity, v, coupled to the rate law. Following a similar approach outlined in a previous paper [12] the master equation in v is Here Q is the a thermicity parameter and is related to the equation of state of the detonation model via the relation The equation in v is coupled to the rate law via Eq. (11). Although not always required to solve the system, one can express equations for the remaining hydrodynamic flow variables in terms of (12), these are and

The variational streamline approximation
As the shape of the streamlines x(ψ, y) is unknown, part of the solution process is to determine their shape. Upstream of the shock, the streamlines are straight and parallel. The oblique shock relations show that the streamlines are deflected at the shock [13]. The primary assumption of the SSA was to assume that streamline curvature was negligible, such that the streamlines remain straight downstream of the shock. The VSA relaxes this assumption, allowing the streamlines to curve in the Detonation Driving Zone (DDZ). In general, we consider a streamline to be described by where F(ψ) is the streamline slope immediately behind the shock. Here G(ψ, y) is a term which represents the curvature of the streamlines and the form chosen will depend upon the reaction rate used in the explosive model. Note if G = 0 then (16) reduces to the SSA model. At this point we can state the following general results: and Integration of Eq. (4) shows that, ahead of the shock, the stream function satisfies where ρ 0 is the ambient density of the explosive, x f is the x-coordinate where the streamline intersects the shock and D is the detonation velocity in the shock attached frame. The shock surface can be parameterised as y = y f (x) i.e.
where s is the slope of the shock and K is the rate of change of the shock slope. From this we can obtain the result We can then apply the chain rule to relate derivatives in ψ to derivatives in terms of the shock front parameterisation. For example, the derivative in y f is Solving the full detonation problem requires an explosive model, comprising of EOS and reaction rate equations, and a suitable choice of the streamline shape ansatz.

Equation of state
The polytropic equation of state is used to model the explosive and reaction products. It is Here q represents the heat release, per unit mass, of the chemical reaction. We use γ = 3 throughout the paper. The polytropic EOS was chosen as it is commonly used to investigate and validate detonation models [11,[14][15][16] which is the primary interest of this paper. The EOS is non-dimensionalised by assuming the characteristic density to be that of the unshocked explosive, and the characteristic speed D C J to be the one-dimensional detonation velocity. The initial density is ρ 0 = 1 and we use strong shock assumption ( p 0 = 0). With D C J = 1 then q = 1/(2(γ 2 − 1)). The thermicity parameter, Q, is found by considering the first law of thermodynamics. Writing the first law (10) in differential form Collecting terms in dρ and rearranging Noting that for the polytropic EOS c 2 = ∂ p ∂ρ then d p − qρ (γ + 1) dλ = c 2 dρ.
Comparing this expression with Eq. (13) it immediately follows that

Reaction rate models
We investigate two reaction rate models with the VSA. First, a power law type which has maximal reaction at the shock, where the SSA has been shown to be predict diameter effect curves with good accuracy compared to DNS [11]. Secondly, the induction type which has zero reaction at the shock and reaches a maximum reaction downstream, this class reaction rate model is more representative of reaction rate models typically used; the SSA was shown to be less accurate compared to DNS for reaction rates of this form [12].

Power law
The power law reaction rate is extensively used when validating detonation models. It is a single term (i.e. reaction progress variable) model and is described by where m, n and α are constants. The choice of α was made such that, for a given choice of n and m, λ = 1/2 coincided with y = 1/2 in the solution to the one dimensional detonation problem. Pressure dependent rate laws are commonly used to develop detonation models. The choice of x(ψ, y) was guided by the following: exothermic heat release has the effect of curving the streamlines inwards [9]; the power law reaction model has maximal reaction at the shock, therefore its effect on the streamline curvature will be greatest near the shock, and reduce as reaction proceeds; in [11,12] it was shown that for such reaction rate models the streamlines curve immediately behind the shock and then straighten through the reaction zone so any suitable ansatz should have maximal curvature at the shock and straighten through the reaction zone. A convenient form of ansatz which will reproduce those features is where is the second derivate of the streamline shape with respect to y, evaluated at the shock. For γ = 3, F(s) = s/(1+2s 2 ) .τ is a free parameter which controls the rate of decay of the curvature term and therefore how much the streamlines are allowed to curve. We require τ > 0 for the streamlines to straighten through the reaction zone. We then have and which match the analytical values obtained from the oblique shock relations at the shock y = y f . In [12] it was shown that, for γ = 3, Note that Eq. (30) is not needed in the SSA approximation since the terms in k are ignored. Setting k = 0 in Eq.
(30) gives a differential equation which, subject to boundary conditions, could be used to solve the shock shape. However, this would be a bad approximation since the equation cannot be solved if there is no reaction at the shock and we are assuming that k = 0, which is almost never true.
Since F, k are functions of ψ, only, we take the derivative of Eq. (26) with respect to ψ To find the derivative of the streamline slope F we first note that F is a function of the shock slope s(ψ). With this we have dF(s) dψ = dF ds For the polytropic EOS with γ = 3 As for the streamline curvature derivative, we know from Eq. (30) that k is a function of s and K . We therefore have where the term on the RHS contains the third derivative y f of the shock front. Here ∂k ∂s K and ∂k ∂ K s are evaluated from Eq. (30). With the ansatz, EOS and reaction rate specified, Eqs. (26)-(34) can be substituted into the master equation (12) and, coupled with (11), form a set of ODEs which can be solved along a streamline ψ = constant. The oblique shock equations provide the state immediately behind the shock for a given shock position (x f , y f ). The flow must also satisfy the generalised CJ condition that the right hand side of equation (12) is zero at the sonic point (u 2 + v 2 = c 2 ). For a given streamline, where x f , y f , y f , y f are known, there will be a unique value of y f which solves the master equation subject to the boundary conditions described, we can therefore write where the function f has a unique solution on each streamline. The problem is subject to the boundary conditions on axis and, for an unconfined rate-stick the requirement that the post shock flow must be exactly sonic requires at the charge edge. Integration from the charge axis to the boundary condition gives the complete shock shape for a particular detonation velocity D and hence the charge width. Diameter effect curves are generated by varying the magnitude of D. On the axis, y f = 0 by symmmetry, so we can solve the master equation for y f to obtain its value on the axis. Note that if we ignore terms in k and its derivatives, then Eq. (26) reduces to the ansatz for the SSA model.

Induction
A reaction rate model with an induction zone (where the reaction starts at zero and increases to a maximum) was chosen to provide a test of the VSA for a reaction model for which the SSA has been shown to be less accurate than for the power-law reaction rate [12]. The equations for the model are where λ 1 represents a fast reaction and λ represents the global burn of the explosive and is coupled to the heat release term in the EOS. For convenience β = 1 and α was chosen such that for a one-dimensional detonation with n = 0, λ = 1/2 when the distance behind the shock was 1/2 also. The variation of λ through the reaction zone has similar properties to the Ignition & Growth reaction rate model commonly used to model explosives [2,17]. We note that the induction model does not include an ignition term.
In the absence of reaction, Eq. (30) shows that the streamlines curve away from the axis immediately behind the shock. Inspection of DNS shows that the streamlines curve inwards at the point where the reaction rate is maximal. Therefore, for rate models which have zero reaction at the shock, the streamline curvature at the shock does not provide any physical insight into the behaviour of the streamlines downstream. A streamline ansatz of the form of (26) is not appropriate for this type of reaction rate. We must consider an alternative approach to modelling the curvature of the streamline. The model should satisfy the following constraints: (i) The curvature of the streamline should be zero on axis (ii) The streamline curvature at the Crocco point should be zero, as the Crocco point is where the streamlines transition from diverging immediately behind the shock to converging.
Analysis of the streamlines from the DNS showed that a quadratic curve gave a reasonably good fit for the streamline shape in the DDZ, both close to the axis and near the charge edge. The streamline ansatz was where We and where s crocco is the shock slope at the Crocco point, and τ is the shape parameter. In addition to satisfying the conditions outlined above, the choice of κ(ψ) has the property that the magnitude of κ grows with the streamline slope, which reflects the streamline shapes we see in the DNS. From (39) we have and from (40) The highest order term in Eqs. (41)-(44) is y f , therefore for a streamline ψ = constant we are required to solve for the induction model. The boundary conditions and solution procedure for the induction model are the same as those described for the power law reaction. The primary difference is that we solve for y f on each streamline, rather than y f .

Error minimisation
Each VSA model has a free parameter, τ , in the streamline ansatz. For each detonation velocity we determine the optimal value of τ by solving the system for multiple values of τ and computing a global error E(τ ). The optimal solution is the one which minimises E(τ ). The corresponding charge half-width is a product of the error minimisation routine.
The global error in the DDZ was calculated by evaluating the local error in the radial momentum and integrating over the entire DDZ domain. At any data point we assume the error is given by the extent to which the radial momentum equation is not satisfied, i.e.
At a given grid point (x i , y i ) where u i , p i etc. were obtained from the VSA solution data. The VSA data was non-linearly spaced and required interpolation onto a uniform grid. Integrating over the entire domain of the DDZ gives the global error, The optimised solution is the value of τ which minimises E. The DNS data is not used as part of the error minimisation routine.

Model accuracy: DNS simulations
As in [12], high resolution simulations of the rate-stick problem were carried out using Falle's MG code. MG uses a second order, conservative, upwind scheme to solve the full time-dependent Euler equations. The HLLC approximate Riemann solver [18] was used to compute the fluxes. MG uses Adaptive Mesh Refinement so that the grid is only refined where necessary. The numerical domain was 0 ≤ x ≤ 3, 0 ≤ y ≤ y u in slab (x, y) coordinates with the explosive in 0 < x < 1 and an inert, low density, confiner in 1 < x < 3. y u was such that the shock remained in the interior of the domain and therefore depended on the particular case. The density of the confiner was ρ c = 0.81. The boundary conditions were symmetric on x = 0 and zero gradient on the other boundaries. In order to prevent reaction ahead of the shock front, the reaction source term was set to zero for p < 0.05.
The radius of the explosive was varied by multiplying the rate constant in the source term by an appropriate factor. The detonation was initiated by imposing a high pressure region in 0 ≤ x ≤ 1, 0 ≤ y ≤ 2 and the detonation was then allowed to propagate until it became steady. The y coordinate of the shock front on the axis was tracked as a function of time and the resulting y − t data used to compute the detonation velocity.
All calculations performed were unconfined. The numerical resolution of DNS was increased until the detonation velocity was converged to three significant figures. In all cases, this ensured that at least 50 cells were in the reaction zone, which has been shown to be required for convergence [19].
In calculations where the detonation velocity was significantly lower than D C J it was found that the subsonic region extended significantly further downstream of the shock than for higher detonation velocities, with some streamlines having more than one transonic point. For cases where the sonic point was far downstream of the shock we have excluded this data from the figures. See appendix A for a screenshot from the DNS and some further discussion. Figure 2 shows diameter effect curves for VSA, SSA and DNS (where computed) for different rate law parameters, where the VSA data points correspond to the minimum error solution. We see that the VSA and SSA are in good agreement for rate law parameters for which the reaction rate is less sensitive to the thermodynamic state (n = 0, 1). Indeed, in Fig. 2b the SSA, VSA and DNS data are almost indistinguishable for larger charge widths. The VSA slightly underpredicts the charge width at lower detonation velocities.

Power law reaction rate without Induction zone
For the most state sensitive rate law (m = 1/2, n = 2), Fig. 2e shows that the VSA offers a clear improvement in accuracy compared with SSA; for a given detonation velocity the VSA is much closer to the DNS data point than SSA. The VSA systematically underpredicts the charge width, for a given detonation velocity, the VSA data points lie to the right of the DNS points. Figure 2e suggests that SSA is in better agreement with DNS near the failure point. Figure 3 shows the DDZ for DNS, VSA and SSA for different detonation velocities. It is clear that the VSA is capturing the shock front accurately compared to DNS up to a charge width of 90% of the width of the explosive, which is a dramatic improvement on the SSA model. Close to the charge edge the VSA overpredicts the shock front curvature which results in the half-width being lower than DNS and explains why the VSA systematically underpredicts the half-width. In Fig. 3a and b, for D = 0.93 and D = 0.859 respectively, the VSA captures the sonic surface in good agreement with DNS, except near the charge edge. In this region curvature and edge effects will be greatest, and the assumed streamline shape is not a good approximation. In Fig. 3c it can be seen that the VSA fails to capture the sonic surface accurately for D = 0.73. We note that the sonic point near axis the is far downstream of the shock and is not included in the figure. Figure 4 shows the DDZ for m = 1, n = 1, and again shows that the VSA is more accurate at predicting the shape of shock and sonic surfaces compared to SSA for velocities close to D C J and for velocities much lower than D C J . The VSA is in excellent agreement with DNS, only diverging at the edge of the explosive. In Fig. 3a, which is D = 0.907, it can be seen that the shock surface predicted by VSA is in excellent agreement with DNS, except near the edge. The sonic surface is also in good agreement. In Fig. 3b and c there is a similar pattern, the shock surface is in good agreement but the sonic surface becomes less accurate as the detonation velocity reduces. Figure 5 shows how the global error of the DDZ for the VSA calculation varies as a function of the shape parameter τ for the different reaction rate parameters . In all of the models considered, we see that 0.3 < τ min < 0.9, which indicates that there is a non-trivial curvature immediately downstream of the shock. Figure 5a-d indicate that as the charge gets smaller the magnitude of τ min (corresponding to the optimal solution) increases, suggesting that curvature effects are greater for smaller charges. Note that larger values of τ result in straighter streamlines. However, in Fig. 5e, which corresponds to n = 2 and the most state sensitive rate law, we see the opposite behaviour, that for lower detonation velocities τ min decreases. We note that Fig. 5c is the only parameter-set that exhibits a turning point in the diameter effect curve. Formally n > 3/2 is required for a turning point [15]. Figure 6a-c show the streamline shapes for different VSA solutions for m = 1/2, n = 2. From visual inspection one can see that the streamlines do exhibit a small, but not insignificant, curvature towards the axis. This agrees with DNS results previously published for power law reaction rates [11,12]. Table 1 shows the computation time required to calculate a single point on a diameter-curve for the different methods. It is clear that for both SSA and VSA, solutions can be computed in a short time. The primary reason that the VSA takes longer than SSA is because solutions for mutliple values of τ must be computed. For DNS, the computational demands are far greater, approximately two orders of magnitude greater than VSA. In the time it takes for DNS to produce a diameter effect curve for a single set of model parameters, the SSA/VSA models could compute diameter effect curves for hundreds of different model parameters.  -c shows that the VSA is able to predict diameter effect curves more accurately than SSA. As for the power law reaction rate, we observe that VSA systematically underpredicts the half-width for a given detonation velocity. For n = 0, Fig. 7a shows that the VSA greatly improves the accuracy of the diameter curve compared to the SSA. For n = 1, 2, Fig. 7b and c shows a small improvement in accuracy. For the stateless reaction (n = 0) Fig. 8a-c shows that the VSA predicts the shape of the shock surface with good accuracy compared to the DNS. However, except near failure, the VSA significantly overpredicts the charge half-width. The VSA and DNS sonic surfaces match well for D = 0.938 and D = 0.875, except near the charge edge, where the VSA is not able to capture the turning point. This is likely due to the ansatz not accounting for the curvature effects at the charge edge, which are more significant at lower detonation velocities. For all detonation velocities, the VSA is more accurate than SSA.
When state sensitivity is introduced, the VSA is less accurate. Figure 9 shows the DDZ for different detonation velocities with n = 1. Figure 9a shows that the VSA captures the shock surface with good accuracy, but for the sonic surface VSA and DNS begin to diverge close to the axis. Figure 9b and c shows that the VSA fails to predict the shape of the sonic surface near the charge edge. For n = 2, which is a more sensitive rate law, it can be seen in Fig. 10a-c that the VSA systematically underpredicts the shock curvature and the sonic surface diverges from DNS not too far from the axis. For these model parameters, the VSA offers little improvement over the SSA. Further investigation would be required to understand whether the accuracy of the VSA could be improved by choosing a different ansatz.

Conclusions
We have demonstrated that introducing curvature effects into a streamline model (VSA) of unconfined steady state detonations offers a marked improvement in accuracy compared to the simpler model where the streamlines were assumed to be straight (SSA). With a polytropic EOS and power law reaction rate, where reaction is maximal at the shock, the VSA is able to capture the shock front and sonic surfaces accurately compared to DNS, only diverging near the charge edge. The VSA is significant improvement over the SSA, which systematically under predicts the shock front curvature and fails to capture the sonic surface compared to DNS.
For reaction rate models which have an induction zone, where the reaction is maximal downstream of the shock, the VSA is more accurate than the SSA in predicting the shape of the shock front and sonic surface when compared with DNS. However, compared to results for the power law reaction rate, the VSA does not perform as well near the charge edge. This could be due to the chosen ansatz only having a single degree of freedom. Adding further parameters to the VSA ansatz may yield more accurate solutions. The idea of approximate detonation models is that they are simple and computationally cheap compared to DNS. Introducing more free parameters to the VSA would increase the computational cost and at some point, the advantages of the approximate model over DNS would be lost. For both power law and induction reaction models, the VSA has been shown to produce more accurate diameter effect curves than the SSA.
Our results suggest that this approximate detonation model could be used for efficient characterisation of explosive models, where we are primarily interested in comparing the diameter effect curve, generated for a set of rate-law parameters, to experimental data. For reaction rates which have multiple free parameters and potentially require hundreds of simulations, the VSA offers improved accuracy over the SSA, whilst still being orders of magnitude more efficient than DNS. Indeed, one could use the VSA to constrain model parameters and then use DNS to fine tune those parameters. Additionally, for explosive models for which the reaction rate is maximal at the shock, the VSA can be used to gain a very good insight into the 2D structure of the DDZ. This value is independent of the detonation velocity, reaction rate and radius of the explosive. Therefore it must depend upon the equation of state only. A good approximation is c λ=1 = γ 2 (γ + 1) 1 2 , which, for γ = 3 gives c = 0.61237.