Numerical modelling of steady detonations with the CREST reactive burn model

Watt et al. [J Eng Math 75(1):1–14, 2012] have shown that one can obtain good results for the propagation of detonation waves in cylindrical charges by assuming that the post-shock streamlines are straight. In this paper, we compare this Straight Streamline Approximation (SSA) to high-resolution Direct Numerical Simulations (DNS) for different models of explosives. We find that the SSA is less accurate for realistic explosion models than it is for polytropic equations of state with power-law reaction rates.

Schematic of detonation in a cylindrical charge. The shock wave (red line) is propagating from right to left at velocity D into the explosive. The flow behind the shock wave is expanding upper limit on the detonation velocity, the Chapman-Jouget (CJ) detonation velocity. For small radii, the detonation will fail to propagate; the maximum radius at which failure occurs is called the critical radius.
Confinement of an explosive enables the detonation to occur at radii below the critical radius for an unconfined charge. This is because the confiner reduces the lateral expansion of the gaseous products, thus maintaining a higher pressure in the detonation reaction zone than in the unconfined case [2].

Modelling
One of the primary motivations for detonation modelling is the high cost of performing experiments. Such models can provide understanding as well as predictions of the behaviour of detonations in more complex geometries than a simple cylindrical charge. This means that the models can be calibrated with simple experiments and then applied to cases that are relevant to applications. This can reduce the need for expensive experiments.
Since detonation in heterogeneous explosives is complex, current models are phenomenological and attempt to capture the burning of the explosive at the continuum level. At this scale, the inviscid reactive Euler equations are assumed to provide an accurate description of the physics. These need to be supplemented by an equation of state for the explosive and its products together with equations describing the chemical reactions [3]. Equations of state are usually multiphase with separate equations for each phase. The phases are assumed to be in pressure (mechanical) equilibrium and moving with the same velocity. Popular choices of equations of state include Jones-Wilkins-Lee (JWL), Williamsburg, polytropic and finite-strain. Apart from the polytropic EOS, these equations of state are empirical and must be calibrated with experimental data. The equations governing the chemical reaction rate are also empirical and are typically calibrated with ignition timing and size-effect data [4,5]. The reaction rates are generally related to the shock strength via a dependence upon one of the thermodynamic variables, such as pressure [6], temperature [7] and entropy [8].
The problem can be solved by numerical integration of the full time dependent Euler equations in a Direct Numerical Simulation (DNS) [9]. Although this produces an accurate solution provided the numerical resolution is high enough, it is computationally expensive. A number of simplified detonation models have therefore been developed, which have much lower computational cost. The first such model was the Wood-Kirkwood model [10], but this requires experimental data to calibrate 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 [11,12] which, at leading order, avoids the need for such calibration. It is based upon the assumption that the detonation velocity is close to the 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 the detonation velocity can be as low as 50% of the CJ velocity [13].
Watt et al. [14] have developed a more sophisticated model, the Straight Streamline Approximation (SSA) model. 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. However, they only considered the simple case of a polytropic equation of state and power-law reaction rate. It is yet to be determined whether the SSA model produces accurate results for different equations of state and forms of the reaction rate. In this paper, we explore the accuracy of the SSA model with a variety of equations of state and expressions for the reaction rate.

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 the depends upon λ and local thermodynamic variables.

Explosive model
CREST is a model for condensed phase explosives with a two-phase equation of state and an entropy dependent reaction rate. The use of an entropy based reaction rate is unique to this model: it being more usual to use pressure or temperature [5,6,11,15]. It has been shown that CREST can accurately reproduce run-time to detonation and radius effect curves [16]. DNS have been carried out using this model [8,[16][17][18], but it has not yet been combined with any steady state detonation models such as the Wood-Kirkwood model or DSD. The two phases in the equation of state are an unreacted solid phase and a fully reacted gaseous phase. The equation of state for each phase can be written in the form e k = e k ( p k , ρ k ) where subscript k denotes the phase, s for the solid and g for the gas.
For the solid phase, where p i (ρ s ) and e i (ρ s ) are the pressure and internal energy, respectively, on the principal isentrope that have been fitted to experimental data, Γ (ρ s ) is the Grüneisen parameter-complete definitions can be found in [16]. The gas phase is modelled by the JWL [19] equation of state: where A, B, R 1 , R 2 , ρ 0s , and w are constants. The Isentropic Solid Equation (ISE) model is used to determine the equation of state of the mixture [20], which assumes that the two phases are in pressure equilibrium. The solid phase is assumed to be on an isentrope, which means that the first law of thermodynamics for an adiabatic process holds: The remaining mixture equations are where λ is the mass fraction of the gaseous phase.
The reaction rate model consists of a fast reaction associated with "hot-spots", with a reaction progress variable λ 1 , which activates a slow bulk reaction, described by λ 2 . The reaction rates [8] are where and Here the c j are constants and Z S is a function of the solid-phase entropy where τ (v s ) is related to the Grüneisen Γ by The reactions are one-way and the chemical energy is released implicitly by the transition from solid to gas. Differentiating equation (5) w.r.t. pressure at constant entropy and λ gives 1 where is the square of the frozen sound speed. The sound speeds of the phases are given by  The parameters are taken from [16] and are reproduced in Table 1 for reference.

Streamline detonation model
Watt et al. [14] have developed an approximate model for detonations in cylindrical charges which make an assumption about the shape of the streamlines and then solves the equations along the streamlines. The simplest approximation is to assume that the post-shock streamlines are straight, the straight streamline approximation (SSA). This gave results which were in good agreement with DNS for both confined and unconfined charges for a polytropic EOS and a power-law reaction rate. In particular, they obtained much more accurate diameter effect curves and failure diameters for both ideal and non-ideal explosives than the standard DSD method. However, although it works well for a polytropic EOS and power law reaction rate it is not known whether this is true for other EOS and reaction rates.

Governing equations in streamline coordinates
Since the method uses the equations along the streamlines, it is convenient to work with the compressible stream function, ψ. In cylindrical polar (r, z) coordinates, this satisfies where u and v are the velocities in the r and z directions, respectively. Equations (18) and (19) give i.e. curves of constant ψ are streamlines as they should be for a stream function.

Master equation in v
The governing equations in streamline coordinates can be reduced to a master equation in the z velocity, v. Downstream of the shock there is no dissipation and the flow is isentropic if there is no reaction. We must therefore have where Q is the thermicity parameter [7]. To determine Q in this case, we start with which follows from the first law of thermodynamics in the absence of dissipation. Since the solid is assumed to be on an isentrope, Eq. (4) gives The gas phase is not isentropic if λ is not constant, so we have dρ g can be found from the differential form of Eq. (5) Using (28), (29) and (30) in (27) gives Since the sound speed of the gas satisfies From Eq. (15), the term in the first set of square braces is just 1/ρ 2 c 2 , so comparing the above equation with (26) we finally get Reducing the governing equations to a master equation in the flow velocity involves successive substitutions of the energy and momentum equations into the continuity equation. We start by expanding the partial derivative w.r.t. z in the continuity equation (21) r ∂ ∂ψ Substituting for u/v using (20) gives This can be rearranged to give an equation for ∂ρ/∂z To get the z derivative of the pressure, we use (22) to eliminate ∂ p/∂ψ in (23) On substituting for u from (20), we get Expanding the final term and rearranging gives Using (33) and (34), the energy equation (26) becomes Rearranging this gives the master equation for the z velocity, v For a given streamline shape, Eq. (35) is an ordinary differential equation for v and can be solved together with Eqs. (25), (33) and (34) for ρ(z), p(z), u(z), v(z) given the boundary conditions at the shock. However, since the streamline shape is not known, we must make some assumptions about their shape.

Straight streamline approximation
The simplest assumption is that the streamlines are straight as shown in Figure 2. In the frame of the shock, the upstream streamlines are parallel to the axis, and the velocity is −D 0 , where D 0 is the detonation speed. At the point (r f , z f ) where the streamlines intersect the shock, they are deflected by the curved shock and are then assumed to remain straight until they reach the sonic locus.
Integrating (19) upstream of the shock gives

Fig. 2
Schematic of a steady-state detonation. In the frame of the shock, the unreacted explosive flows into the shock wave in the negative z-direction. The SSA assumes that the flow behind the shock has straight but diverging streamlines. The streamlines intersect the shock at the point (r f , z f ). The flow speed in the shock frame increases as the chemical reaction proceeds and goes sonic at the sonic locus from which we get A streamline at radius r f intersects the shock at z = z f as shown in Fig. 2. Since it is straight, we have where F(ψ) is the slope of the streamline. This is where u f and v f are the components of the flow velocity immediately behind the shock. We also have since the streamlines are straight. From these results, we can obtain (∂r/∂ψ) z and ∂ 2 r/∂ψ∂z, which are required by (35). From (38), we have where Using this and (37) in (41), we get Taking the derivative of Eq. (41) gives Equations (43) and (44) show that the equations along the streamline will depend upon its point of intersection with the shock, (r f , z f ), the shock slope z f and its second derivative, z f . Hence, for a given intersection point and slope, there will be a unique value of z f that solves the eigenvalue problem given by the master equation (35). We therefore have where the function f is given by the solution to the eigenvalue problem along the streamline. The solution of this equation with the boundary condition gives the shock shape and hence the complete solution.
Special care needs to be taken with the master equation (35) on the axis (r = 0). The first term on the right hand side is 1 r which is ill defined on the axis. Using L'Hopital's rule we have On the axis both (∂ 2 r/∂ψ∂z) and (∂r/∂ψ) −1 z are non-zero, therefore in the case where r = 0 this expression was used in place of (1/r )(∂r/∂z) ψ in the master equation.
The derivative dF/dz f in (44) is found by taking the derivative of (39). With CREST the shock relations must be solved numerically, thus dF/dz f must be computed numerically. A central difference approximation was used in this case.

Solving the system
To solve the system, we choose a detonation velocity D 0 and then compute the shock surface by integrating equation (45) from the axis to the charge edge using the boundary conditions (46). For an unconfined charge, the charge edge is determined by the condition that the post-shock flow be sonic [11]. If the charge is confined, then the condition at the charge edge must be determined by a shock polar analysis. In either case, this gives D 0 as a function of charge diameter.
As |z f | increases, both the post-shock Mach number and streamline deflection increase monotonically until the streamline deflection reaches a maximum at the Crocco point [21]. For larger |z f |, the streamlines converge behind the shock, which means that dF/dψ is negative and hence from (44) so is ∂ 2 r/∂ψ∂z. The master equation, (35), has no solution in this case so that one cannot integrate to the sonic point when the Crocco point occurs before the sonic point. This is the case for the CREST EOS, so we simply assumed that the charge edge was at the Crocco point for unconfined charges. The impact on the accuracy of the results is discussed in Sect. 5.3.

DSD model
We also compared our diameter effect curves with those obtained using a leading order DSD model [12,13,15], which is described and critically analysed in [15]. This is a one dimensional model for which the equations given in Appendix B in [12] dρv where −z is the post-shock distance along the shock normal, v is the velocity normal to the shock, D n is the normal velocity of the shock and κ is its local total curvature. Using (26), we can derive the master equation for this model These equations can be integrated from the shock downstream and for a given normal shock velocity D n , κ is determined by the condition that the solution passes smoothly through the sonic point v = c. This gives a relation D n (κ) between the detonation speed and the local shock curvature. The total shock front curvature and normal shock velocity are given by for axisymmetry [22]. Here s = z f (r ) is the shock slope and s = z f (r ) is its derivative. This provides a relation for z f of the same form as (45), which can be integrated to give the shock shape. For a given detonation velocity, the charge radius is determined by integrating along the shock front from the charge centre out to the edge. As for the SSA, the Crocco point was used as the boundary condition at the charge edge.
A solution can only be obtained for D n if κ is below a critical value, κ c , at which point D n = D c . For κ < κ c there are two solutions for D n , an upper stable branch with D n > D c and a lower unstable branch with D n < D c [23]. Since D n decreases as we integrate out from the axis, it may happen that it becomes equal to D c and we have to start using the value of κ from the unstable branch. Since this occurs near the failure radius, it is clear that this model is somewhat dubious for such cases.

Numerical method
In order to explore the accuracy of the SSA, high-resolution numerical calculations of the cylindrical charge problem were performed using Falle's MG code [24]. This uses a second-order, conservative, upwind scheme to solve the full time-dependent Euler equations. The HLLC approximate Riemann solver [25] was used to compute the fluxes since it maintains stationary contacts and slip lines. This is important for 2D calculations in which there is a sharp interface between the explosive and the external medium. For reasons to be explained later, it is necessary to use a reasonably diffusive Riemann solver.  N grid  levels, G 0 , G 1 , . . . , G N −1 with mesh spacings h n = h 0 /2 n , where h 0 is that on G 0 . G 0 and G 1 cover the whole domain, but finer grids need not do so. Refinement is on a cell-by-cell basis and is controlled by error estimates based on the difference between solutions on different grids i.e. the difference between the solutions on G n−1 and G n determines refinement to G n+1 . An additional refinement condition, based upon the chemical reaction rate, was included to ensure that the solution was always resolved to the finest level in the reaction zone.
Like all shock-capturing codes, MG works with the conserved variables, (e, ρ, λ), but the pressure and densities of the phases are needed to compute the fluxes and source terms. These were obtained from the pressure equilibrium condition (7) and mixture equations, (5), (6) using a secant method iterating on the solid density. To prevent reaction ahead of the shock front, the reaction source term was set to zero for p < 0.05 Mbar. For λ < 0.02, the solid-phase equation of state was used; otherwise, the mixture model was used.

One-dimensional calculations
The code was tested by comparing one-dimensional calculations with a steady-state ZND calculation. The detonation was initiated with a high -pressure region and allowed to propagate for a long time to ensure a steady state had been obtained. The results are summarised in Table 2, from which it can be seen that CJ detonation velocities are in excellent agreement. Figure 3 also shows that the agreement between the structure of the reaction zone from the DNS calculations and that of the steady-state solution for EDC37 is equally good.

Two-dimensional calculations
The numerical domain was 0 ≤ r ≤ 2, 0 ≤ z ≤ z u in cylindrical polar (r, z) coordinates with the explosive in 0 < r < 1 and an inert, low density, confiner in 1 < r < 2. z u was such that the shock was in the interior of the domain and therefore depended on the particular case. For the unconfined cases, the confiner density was ρ c = 0.9 g cm −3 for EDC37 and ρ c = 1.5 g cm −3 for PBX9502, which are small enough for the expansion fan at the charge edge to be transonic. The boundary conditions were symmetric on r = 0 and zero gradient on the other boundaries.
The radius of the explosive in physical units 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 ≤ r ≤ 1, 0 ≤ z ≤ 2 and the detonation was then allowed to propagate until it became steady. The z coordinate of the shock front on the axis was tracked as a function of time and the resulting z − t data used to compute the detonation velocity. Table 3 shows the results for EDC37 for different numerical resolutions. It can be seen that 320 cells/mm are required to obtain well-converged results for solid-phase entropy and the detonation velocity. This corresponds to 69 cells in the reaction zone which is consistent with previous work in which it was found that 50 cells were required for convergence [13].   The table shows the values of the reaction progress variable, λ shk at the end of the numerical shock structure. These show that somewhere between 160 cells/mm and 80 cells/mm there is a transition from pure solid phase in the post-shock state (λ shk ≤ 0.02) to one in which there is a mixture of gas and solid. In the CREST model, the reaction rate depends on the solid entropy, so that it is important that the post-shock state be pure solid with the correct entropy. It is also essential that the entropy varies smoothly along the shock in the two-dimensional calculations. This was achieved by adding extra dissipation to the HLLC solver as described in [26] in order to widen the shock structure. The accuracy, therefore, cannot be improved by using a method that gives a narrower shock structure, such as a higher-order scheme. Note that, in CREST, the extent of the reaction depends upon the solid entropy, and there is no guarantee that the reaction will go to completion, i.e. λ max = 1. It is for this reason that λ max varies significantly with resolution when there are less than 160 cells/mm. This is not a feature observed for reaction rates which are always positive for λ < 1. Figure 4a shows the diameter effect curve for EDC37 with the CREST reaction rate given by Eqs. (8) to (12). It can be seen that the DNS results lie between those from the DSD approximation and SSA. This is consistent with results obtained by [14], who found that the SSA model gives too large a detonation velocity. At small radii the disparity between the DNS results and the SSA is significant and SSA gives a failure radius of (∼ 0.83 mm) compared to (∼ 1.25 mm) for the DNS calculations. These differences are much larger than those found using a polytropic EOS Entropy refers to the unreacted equation of state, λ max is the maximum value of the reaction progress variable, D is the detonation velocity, CIRZ is the Cells in the Reaction Zone and is the number of cells on the finest level in which the reaction rate is significant and λ shk is the value of λ at the downstream end of the shock structure and power-law reaction rate, for which the difference between the failure radii was less than 5%. Figure 4b shows that similar results are obtained for PBX9502. The streamlines in the DNS calculation shown in Figure 5 are weakly curved throughout the reaction zone and they curve towards the axis except near the shock and the charge edge, where they curve away. As is shown in Appendix A and in [11], the shock curvature causes the streamlines immediately behind the shock to curve away from the axis, whereas the subsequent reaction curves them towards it [11]. As can be seen from Fig. 6, the reaction source term is small at the charge edge, while the shock curvature is significant. Figure 7 is a magnified view of the region near the charge edge from Figure 5 and it can be seen that the streamlines curve towards the axis at smaller radii where the chemical reaction is significant. The finite width of the numerical shock structure makes it very hard to see the streamline curvature at the shock, but it does affect the streamline slope downstream of the shock. At the charge edge the Prandtl-Meyer rarefaction curves the streamlines away from the axis.

Power-law reaction rate
To investigate the effects of a different forms of reaction, calculations were performed for the CREST EDC37 EOS with a simple power-law reaction rate of the following form: where α = 2.84 µs −1 to maintain consistency with the remaining equations in the model. A state insensitive equation was chosen so that the reaction rate would be independent of shock strength. The main difference between this reaction rate and that of the CREST model is that the reaction rate is a maximum at the shock and monotonically decreases, whereas for CREST the reaction has an induction zone immediately behind the shock before reaching its maximum value a significant distance behind the shock. Figure 8 shows that the SSA and DNS results are in excellent agreement for large charge radii and, although there is some disagreement for smaller charge radii, it is much smaller than for the CREST reaction rate. This is because, as we show in Appendix A, the large reaction rate at the shock now partly compensates for the streamline curvature due to the shock curvature, as can be seen from Figs. 9 and 10. As a result, the streamlines are straighter and their slope is quite close to that immediately behind the shock i.e. the assumptions of the SSA are more nearly satisfied. We would therefore expect SSA to work well for the widely used ignition and growth model (e.g. [27]) since this has a non-zero reaction at the shock.

Confined calculation
Calculations were also performed with the CREST EDC37 model with confinement. Confinement reduces shock curvature and also prevents the flow from being sonic at the charge edge. This enables detonation waves to propagate in charges with a much smaller radius than in the unconfined case [2]. The reduced curvature means that the streamlines are straighter and the SSA should be more accurate than in unconfined detonations. The polytropic equation of state was used for the confiner with a stiffened polytropic gamma (γ = 5) and initial density of ρ 0 = 7 g cm −3 . The explosive and confiner were initialised with the same pressure. The boundary angle used in the SSA was determined by a shock-polar match between the explosive and inert. Figure 11 shows that the agreement between SSA and DNS is much better than in the unconfined case shown in Fig. 4a. In contrast, the DSD is not significantly better than in the unconfined case.

Failure to integrate to the charge edge
In Sect. 3.4, we showed that it is not possible to integrate to the charge edge for unconfined charges when the Crocco point occurs before the sonic point. Figure 12 shows how the post-shock Mach number varies with streamline deflection for EDC37, ANFO and polytropic equations of state. The ANFO equation of state was a pseudo-polytropic equation of state in which gamma is a quadratic function of the density [15]. It is clear that the sonic point occurs after the Crocco point both for ANFO and the solid-phase equation used by CREST, whereas they coincide for a polytropic EOS [28]. Figure 12 suggests that for detonation modelling of unconfined explosives with the SSA Fig. 9 Subsonic region (blue) from DNS calculation with EDC37 EOS and single-term power-law reaction rate: charge radius R c = 2.5 cm with streamlines superimposed. The explosive-confiner interface is indicated with a bold line. The detonation velocity is D = 0.836 cm µs −1 and D C J = 0.877 cm µs −1 Fig. 10 Reaction source term from DNS calculation with EDC37 EOS and single-term power-law reaction rate: charge radius R c = 2.5 cm with streamlines superimposed. The explosive-confiner interface is indicated with a bold line. The detonation velocity is D = 0.836 cm µs −1 and D C J = 0.877 cm µs −1 the Crocco point will define the boundary condition at the edge of the rate stick because for any realistic EOS the Crocco point is at a shock angle smaller than the angle required for sonic flow post shock. Figure 13a, b shows the shock slope as a function of radius for D = 0.861 cm µs −1 (near CJ) and D = 0.760 cm µs −1 (near failure) for EDC37 with the CREST reaction rate. Both cases show the same general trend: the SSA and DNS are in good agreement for r < R/2 and thereafter the shock slope for SSA increases more rapidly than for DNS. Inspection of Fig. 6 in conjunction with Fig. 13b shows that the radius at which the SSA and DNS begin to diverge closely matches the radius at which the peak reaction rate reduces significantly compared to the rate on axis. As the reaction rate is smaller here, then curvature effects will more significantly affect the flow. The DNS shows that the curvature is large near the charge edge, which means that the Crocco point occurs at almost the same radius as the sonic point. The discrepancy between SSA and DNS in Fig. 4 must therefore be due to streamline curvature, rather than our inability to integrate past the Crocco point.  Fig. 12 Post-shock Mach number versus the streamline deflection shock polar for EDC37, polytropic (Poly) and polytropic with quadratic gamma (ANFO). D = 0.813 cm µs −1 was used for EDC37 and the remaining curves are independent of the detonation velocity. Note that for the polytropic EOS the sonic and Crocco points coincide Figure 14a, b is similar to Fig. 13a, b except that the reaction rate is a power-law one. For a detonation speed near to CJ (Fig. 14a), the SSA and DNS are in reasonable agreement, but for the lower detonation speed shown in Fig. 14b, the SSA shock slope decreases more rapidly than the DNS. This suggests that for a power-law reaction rate, curvature effects only become significant at lower detonation velocities, which is why the agreement between the SSA and DNS diameter effect curves in Fig. 8 is worse at lower detonation velocities. The main difference between Figs. 14a, b and 13a, b is that the shock curvature diverges between the SSA and DNS at a radius much closer to the charge edge with the power-law reaction compared to the standard CREST reaction model. This suggests that for modelling with CREST curvature effects are significant at all detonation velocities and in a greater portion of the reaction zone than for a power-law reaction.

Streamline curvature
Since the SSA assumes that the streamline curvature is zero between the shock and the sonic loci, we need to examine the validity of this approximation. First consider streamline curvature immediately behind the shock for a polytropic EOS with γ = 3. This is 2 3 , from Eq. (77) in Appendix A. Here s is the shock slope. From [11], we have s 2 = (γ − 1)/(γ + 1) = 1/2 at the sonic point, so s 2 < 1/2. Note that with W = 0, this agrees with the expression given in [21] which does not include reactions.
The first term in the square braces is proportional to the reaction rate and is strictly positive. The second term, which is proportional to the shockwave curvature κ, is of the opposite sign since κ < 0 and s 2 < 1/2. This is also true for the CREST solid-phase equation of state, although in this case the post-shock streamline curvature must be found numerically from Eqs. (68) and (55) in Appendix A. In either case, the effect of the reaction term is to reduce the streamline curvature.
where α, m and n are constants [14]. These show the post-shock streamline curvature as a function of shock slope for SSA calculations in slab geometry for different detonation velocities. It can be seen that in the absence of reaction at the shock, the streamline curvature is a positive and monotonically increasing function of shock slope, whereas it is non-monotonic and significantly smaller in magnitude when there is a reaction. A reaction rate that is a maximum at the shock is also smaller in the rest of the driving zone and this reduces the streamline curvature downstream of the shock. A careful comparison of the streamlines in the reaction zone for the CREST reaction rate in Fig. 6 and the power-law reaction rate in Fig. 10 shows that the curvature of the streamlines is greater for the former. The SSA assumption that the streamline slope is equal to that at the shock and is constant thereafter is therefore significantly more accurate for a power-law reaction rate. Comparing Fig. 15a, c shows that the magnitude of the streamline curvature increases signficantly at smaller detonation velocites; with reaction for a detonation speed D = 0.45, K = 0.6 at the edge of the rate stick whereas for D = 0.97, K = 0.16 at the edge.  Figure 16 shows a comparison between the post-shock streamline curvature without reaction for the CREST solid-phase equation of state and a polytropic equation of state. It is not possible to consider the effect of reaction on the streamline curvature at the shock for CREST as reaction is strictly zero here. For CREST the results were calculated numerically. In both cases, the streamline curvature has a maximum, but for a polytropic equation of state it is independent of shock velocity and becomes zero at the sonic point as expected from Eq. (77). In contrast, for CREST it depends on the shock velocity and becomes zero before the sonic point is reached. This coincides with the point where the streamlines begin to converge post-shock and is the maximum shock slope that can be reached with the SSA. The magnitude of the streamline curvature for the polytropic EOS is significantly higher than for CREST, approximately a factor of two at the peak. This suggests that streamline curvature effects would be more significant for polytropic EOS if there was no reaction at the shock.

Conclusions
We have compared the predictions of the SSA with DNS for the CREST model and have shown that it gives accurate values of the detonation velocity for unconfined cylindrical charges whose diameter is significantly larger than the failure diameter. However, for smaller diameters, it performs poorly compared with previous results for polytropic equations of state with a power law reaction rate [14].
There are two reasons for this. The first is that the CREST reaction rate is zero at the shock, whereas it is a maximum there for a power-law reaction rate. Since reaction at the shock reduces the curvature of the streamlines for a given shock curvature, this means that the assumption that the streamlines are straight with slope equal to that at the shock is less well satisfied for CREST. Secondly, for the CREST equation of state, the sonic point on the shock polar occurs after the Crocco point at which the streamline deflection is a maximum. Since the SSA requires  The results are much better for confined charges since the shock curvature is smaller and there is no difficulty in integrating to the charge edge in this case. Although the properties of explosives are usually determined from experiments with unconfined charges, the confined case is the most relevant to applications. The SSA approximation is therefore still usable for the CREST explosive models and others with similar properties.
The only way to improve the accuracy of SSA is add the effect of curvature. A possible approach is to introduce some curvature of the streamlines that depends upon a set of adjustable parameters. The complete solution can be constructed for some initial guess for the values of these parameters on each streamline. The extent to which this solution fails to satisfy the equations can then be used to construct a global error. Minimization of this error determines the optimal values of the parameters on each streamline. This Variational Streamline Approach (VSA) is clearly more expensive than SSA, but can be made significantly more accurate, particularly if there are several parameters on each streamline. Preliminary calculations with a single parameter determining the curvature near the shock show considerable promise.
Finally, we have been entirely concerned with the verification of the SSA with CREST by comparison with DNS, not the validation of the CREST model by its ability to predict the results of experiments. As with any model, CREST must be properly calibrated to other experimental data, and the SSA can be useful for this. Even for cases in which it is not very accurate, it can provide initial guesses for the parameters in explosion models which can then be refined with DNS calculations.
Since this a local analysis, it is convenient to use Cartesian coordinates (x, y) with the streamline given by x = x s (y) and the upstream flow in the negative y direction. The slope, M, and curvature, K , of the streamline are then Immediately behind the shock where u and v are the post-shock velocities in the x and y directions, respectively. Using (52) and (54) in Eq. (53) gives The partial derivatives on the RHS of this equation can be determined from the Euler equations: u · ∇ρ + ρ∇ · u = 0, ρu · u + ∇ p = 0, u · ∇e − p ρ 2 u · ∇ρ = 0, u · ∇λ = W.
In Cartesian coordinates, these are The energy and reaction equations, (59) and (60), can be combined with the equation of state e = e( p, ρ, λ) and the expression for the sound speed (16) to give To proceed further, we need the derivatives along the shock. If the shock position is given by y = y f (x), then its slope, s, and curvature, κ, are For an oblique shock, the post-shock flow variables are a function of the shock slope only. Therefore, for a flow variable g = g(s), we have dg = dg ds from Eq. (63). Along the shock, dg = (dx, dy) · ∇g = dx (1, s) · ∇g.
Equating (64) and (65) gives In Cartesian coordinates, the operator on g is Given g(s) for each of the variable u, v, p, ρ from the shock relations, this gives four equations for their partial derivatives with respect to x and y. These, together with Eqs. (56)-(58) and (61) give eight linear equations for the partial derivatives of u, v, p, ρ. These can be written in the following form: and For a given equation of state e = e( p, ρ, λ) and reaction rate W , this system can be solved for the vector z, and hence, the partial derivatives of the velocities required to determine the streamline curvature from Eq. (55).
For the CREST equation of state, it is not possible to obtain an explicit expression for the post-shock state variables in terms of the shock slope. However, it is possible for a polytropic equation of state of the form: The shock relations then give Here we have used the strong shock approximation that the upstream pressure vanishes and have also set the upstream density to unity. Using (72)-(75) in (70) gives Substituting the corresponding solution to (68) into (55) gives the streamline curvature for a polytropic equation of state K = s (γ − 1) (γ + 1) 3 1 + s 2 3 qW − 4κ D 3 (2γ − 1) 1 − γ + s 2 (γ + 1)