Reactive burn models and ignition & growth concept

Plastic-bonded explosives are heterogeneous materials. Experimentally, shock initiation is sensitive to small amounts of porosity, due to the formation of hot spots (small localized regions of high temperature). This leads to the Ignition & Growth concept, introduced by Lee and Tarver in 1980, as the basis for reactive burn models. A homogenized burn rate needs to account for three meso-scale physical effects: (i) the density of active hot spots or burn centers; (ii) the growth of the burn fronts triggered by the burn centers; (iii) a geometric factor that accounts for the overlap of deflagration wavelets from adjacent burn centers. These effects can be combined and the burn model defined by specifying the reaction progress variable λ = g(s) as a function of a dimensionless reaction length s(t) = rbc/ bc, rather than by specifying an explicit burn rate. The length scale bc(Ps) = [Nbc(Ps)]−1/3 is the average distance between burn centers, where Nbc is the number density of burn centers activated by the lead shock. The reaction length rbc(t) = ∫ t 0 D(P(t′))dt′ is the distance the burn front propagates from a single burn center, where D(P) is the deflagration speed as a function of the local pressure and t is the time since the shock arrival. A key implementation issue is how to determine the lead shock strength in conjunction with a shock capturing scheme. We have developed a robust algorithm for this purpose based on the Hugoniot jump condition for the energy. The algorithm utilizes the time dependence of density, pressure and energy within each cell. The method is independent of the numerical dissipation used for shock capturing. It is local and can be used in one or more space dimensions. The burn model has a small number of parameters which can be calibrated to fit velocity gauge data from shock initiation experiments.


Introduction
Plastic-bonded explosives (PBX) are heterogeneous materials consisting of explosive grains of various sizes and shapes, a polymeric binder to hold them together and a small amount of porosity; see Fig. 1.Due to the heterogeneities, shock compression gives rise to localized spatial regions of high temperature known as hot spots.Since chemical reaction rates are strongly temperature dependent, reaction in and around hot spots dominate the burning in a shock-to-detonation transition.
Pore collapse is a key hot-spot generation mechanism for shock initiation.Experiments have shown that initiation is sensitive to the meso-scale structure and in particular to porosity.The scale of the heterogeneities in a PBX (tens of microns) is smaller than the cell size that can reasonably be used to simulate a typical high explosive (HE) application.Therefore, a homogenized or volume averaged burn rate is needed to account for the sub-grid burn physics.
The underlying hot-spot reaction physics is best described by the Ignition & Growth concept reviewed in sec.2. We then present a new formulation for a burn rate in sec.3 that better captures the underlying physics of this concept.A key idea is to base the burn model on the strength of the lead shock.In sec. 4 we present a robust algorithm based on the Hugoniot jump condition for detecting the lead shock.A comparison of numerical results with velocity gauge data for shock initiation experiments of PBX 9502 are shown in sec. 5. Final remarks and extensions of the model formulation are briefly described in sec.6.

Ignition & growth concept
Currently available burn models are empirical.They work well only for the limited regime in which they are calibrated.The more physical hot-spot models are based, at least heuristically, on the Ignition & Growth concept introduced by Lee and Tarver in 1980.The ignition phase corresponds to the formation of hot spots, their reaction and thermal runaway on a fast time scale.The growth phase corresponds to the subsequent burning from reactive wavelets triggered by the burned hot spots.Typically, hot spots represent a small fraction of the HE volume and almost all of the burning occurs in the growth phase.
Leading examples of this type of model, or models that can be interpreted with the ignition and growth concept, along with the form for the burn rate are listed below: i Ignition & Growth model, Lee and Tarver (1980).
Rate is a function of local pressure and reaction progress variable.The functional dependence on the reaction progress variable can be interpreted as proportional to the average burn front area of the reactive wavelets emanating from the hot spots.ii DAGMAR, Wackerle et al. (1978) and Anderson et al. (1981).
Rate is a function of lead shock pressure, equilibrium mixture temperature and first order in reaction progress variable.Model parameters were calibrated to embedded Manganin stress gauge data.iii JTF model, Johnson et al. (1985).
Multi-step reaction with the rates a function of the lead shock pressure, equilibrium mixture pressure and first order in reaction progress variable.Noteworthy for comparison with embedded magnetic velocity gauge data; see Vorthman et al. (1985).iv WSD model, Wescott et al. (2005).This is a generalization of model i above.The shock density is used to switch between rates calibrated to the ignition and propagation regimes.In each regime the rate depends on the local pressure and the reaction progress variable.

00003-p.2
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter v CREST model, Handley (2006).
Rate is a function of reactant entropy for two-phase model with pressure but not temperature equilibrium.Two reaction progress variables are used to tune reaction profile based on analysis of embedded magnetic velocity gauge data; see James and Lambourn (2006).
At the time that the DAGMAR and JTF models were developed, artificial viscosity was used for shock dissipation.The lead shock was determined by scanning the viscous pressure profile for the first local minimum.This algorithm had difficulty when a compressive wave followed the lead shock.Due in part to the numerical issue with robustness, these models never gained favor.
The implementation of the WSD model in (Stewart et al. 2007, § 2.4) detects the time of arrival of the lead shock based on the maximum ∂ t P for each cell.Then a level set algorithm applied to the shock arrival time is used to determine the shock speed, which determines the shock state.This may not be too accurate but is good enough to use as a switch on the initiation regime or growth regime.
For the CREST model, the solid entropy can be thought of as a measure of the lead shock strength since the chemical energy released is apportioned to the products.There are issues with shock capturing in two-phase models since the component energy equations are not in conservation form; see Petitpas et al. (2009) for details.
Burn rates based on the lead shock strength naturally account for the phenomena of shock desensitization.With the possible exception of the CREST model, different rate parameters are needed to account for the initial temperature or initial porosity of a PBX.In effect, a PBX at a given initial condition is treated as a distinct explosive.To overcome this and other deficiencies, a burn model needs to better account for the sub-grid or meso-scale burn physics.
Based on the ignition and growth concept, three key physical effects need to be incorporated in a burn model: (i) the density of active hot spots or burn centers, which depends on the lead shock strength; (ii) the growth of burn fronts triggered by burn centers, which depends on the local deflagration speed; (iii) a geometric factor that accounts for the overlap of deflagration wavelets from adjacent burn centers, which depends on material heterogeneities that determine the distribution of hot spots.

New formulation
For our formulation, Scaled Unified Reactive Front (SURF) model, the volume average of the hot-spot reaction is accounted for by directly specifying the reaction progress variable λ as a function of time.The effective burn rate is then determined implicitly as R = (d/dt)λ. (1) The key point is that burn centers are formed and reaction begins only after the lead shock pass over potential hot spot sites; such as a void or an embedded high density particle.Moreover, the average amount of burn depends on the distribution of hot spots, whether a hot spot evolves into a burn center and how a burn center interacts with its neighbors.First, let us define the shock functions t s (x) = time the lead shock arrives at point x, P s (x) = pressure of lead shock at point x.We express the reaction progress variable as (2) where P(x, t) = local pressure, g(s) = reaction scale function, s(t, Ps, P) = scaled burn center reaction length = r bc (t, P)/ bc (Ps), r bc (t, P) = reaction length for single burn center, bc (Ps) = distance between burn centers.

00003-p.3 EPJ Web of Conferences
For now we neglect the initial hot-spot volume, and express the reaction length for a single burn center as which depends on the local deflagration speed, D(P).The average distance between burn centers can be written bc where N bc (P s ) is the number density of burn centers or those hot spots activated by the lead shock that burn on a fast time scale.We expect the reaction scale function g to depend on the statistical distribution of hot spots.It accounts for the interaction of burn centers and the geometric factor for the burn front area due to the overlap of burn wavelets from each burn center.For a given PBX composition, g would depend on the structure of heterogeneities -grain size distribution, morphology, imperfections, impurities, voidsthat give rise to hot spots.
To summarize the SURF model: it is determined by specifying three functions g(s), N bc (P s ) and D(P) that separate the sub-grid burning into distinct parts, each of which has a clear physical interpretation.

Fitting forms
For independent randomly distributed point burn centers in three dimensions, the reaction scale function can be shown to be g(s) = 1 − exp(−s 3 ); see Nichols and Tarver (2002) and Hill et al. (2009).This includes the geometric effects of outward hot-spot burning and inward grain burning similar to the functional dependence of the burn rate on the reaction progress variable, λ 2/3 (1 − λ) 2/3 , in the original ignition & growth model of Lee and Tarver (1980).
The geometry of burn wavelets triggered by the burn centers depends on the dissipative mechanism generating the hot spots.We use a form corresponding to cylindrical burn centers g(s) = 1 − exp(−s 2 ). (5) Measurements on the regression rate for propellants and higher pressure measurements with diamond anvil cells, see for example Esposito et al. (2003), suggest that the deflagration speed is proportional to the pressure.We express it as The sensitivity of distance-to-detonation with shock pressure suggests that the number of burn centers rises rapidly with shock strength.We use a shock strength factor of the form where P ref and t ref are units for pressure and time.The scaled reaction length is then expresses as For these choices, the SURF model has only 4 parameters (A, B, P s,max , n).The model can be calibrated to shock initiation experiments by fitting to embedded velocity gauge data; i.e., Lagrangian time histories of the particle velocity for several points along a shock-to-detonation trajectory; for details see Shaw and Menikoff (2010).We note that s = 1 corresponds to 1 e-fold of reaction.The time scale for this to occur is roughly τ = 1/ f (P s ), and the corresponding burn center length scale bc = D(P s )/ f (P s ).For PBX 9502, used to illustrate the model in sec.5, the reaction time scale for a shock-to-detonation transition varies between 10 ns and 1 μs.Assuming D(P s ) is on the order of 0.1 mm/μs, the burn center length scale would be between 1 and 100 microns.This is comparable to the grain size for the PBX.
The effective burn rate for the fitting form corresponding to Eq. ( 5) can be expressed as since s 2 = − ln (1 − λ).We note that the λ dependence of the rate is very similar to that in the CREST model; see (Handley and Lambourn 2009, Eq. (2)).Moreover, Eqs. ( 5) and ( 8) are similar to the history variable reactive burn model in the CTH code; see description in (Starkenberg 2002).

Lead shock detection algorithm
A numerical shock profile mimics that of a steady-state solution to the Euler equations with viscosity and heat conduction.In the (V, P)-plane, the continuum profile lies between the Rayleigh line and the Hugoniot locus, as shown in Fig. 2; see Gilbarg (1951).The distance between a point on a shock profile and the Hugoniot locus is roughly proportional to the Hugoniot function, whose zero level set defines the Hugoniot locus.The subscript '0' denotes the state ahead of the shock.
For the lead shock, the ahead state is well defined; the ambient or initial state of the HE.For a given cell in a numerical simulations, the shock state can be determined by calculating the time dependent Hugoniot function, h(t) = h (V(t), e(t) , (11) 00003-p.5 and looking for either a zero crossing or a local minimum.The former case occurs when the shock is followed by a compression wave and the latter when the shock is followed by an expansion wave.Following waves would cause the shock strength to change and the wave front to either accelerate or deaccelerate.

EPJ Web of Conferences
The ability of the algorithm to detection the lead shock is shown in Fig. 3.In this example, a piston is used to drive a shock in non-reacting PBX 9502.An accelerating piston strengthens the lead shock, while a deaccelerating piston causes the shock to decay.The algorithm works equally well for the case in which reaction accelerates the lead shock.Possibly the error in one of the other shock jump conditions, such as, δ = (u s − u 0 ) 2 − (P s − P 0 )(V 0 − V s ), can be used to estimate the error in the shock state.
Three key features of the shock detection algorithm are: (i) It is independent of the numerical dissipation and can be used with a shock capturing algorithm based on either artificial viscosity or a Riemann solver, i.e., Godunov scheme.(ii) It is local to a cell and independent of the shock direction, hence can be applied in one or more space dimensions.(iii) It determines time of shock arrival along with the shock state.Consequently, the shock detection algorithm is well suited for our model formulation.

Numerical results for PBX 9502
To illustrate the effectiveness of the SURF model for shock initiation we compare simulated results with data from gas gun experiments.A series of shock-to-detonation experiments were carried out on PBX 9502 by Gustavsen et al. (2006) in which the initial shock pressure was varied.The gas gun drives a well supported planar shock in the PBX.Hence the flow is one dimensional.Embedded magnetic velocity gauges were used to measure Lagrangian time histories for up to ten particles.This provides information on the flow behind the lead shock as well as the trajectory of the lead wave.
Numerical simulations were performed utilizing the Amrita framework Quirk (1998a,b).The solver algorithm is a Godunov scheme on a 1-D Lagrangian mesh.Both the projectile from the gas gun and the PBX were included in the calculation.The measured projectile velocity was used as the initial condition to set the initial shock pressure.
Adaptive mesh refinement allowed shock fronts and the reacting regions to be well resolved.A noteworthy point of the shock detection algorithm is that the Hugoniot function variable h enabled the refinement criterion to select the finest mesh for the lead shock, and to prevent reaction from occurring in the numerical shock profile.Coarsening the mesh by 1 level behind the lead shock is a smoothing operation that damps out numerical variations in the shock state due to the discretization of the shock profile.

00003-p.6
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter A comparison of the velocity gauge data for two experiments and numerical results are shown in Fig. 4. The initial shock pressure in shot 2s57 is 13.5 GPa and in shot 2s68 is 11.5 GPa.We note that initiation is sensitive to the shock pressure.The relatively small pressure difference (15%), changes the transition time from 1.3 μs to greater than 2.0 μs (50%).In addition, the timing of the jump off velocity of the gauges depends on the shock velocity.Hence accurate equations of state are needed for the reactants and products of the PBX.The one we are using is described in (Menikoff 2009).
A noteworthy feature of the data is displayed by the first gauge which is at the PBX interface.After the shock, the particle velocity is nearly constant.This is compatible with a small initial hot-spot volume and burning resulting from the subsequent growth of the hot-spot reaction length; i.e., r bc (t, P).This is qualitatively different then the original Ignition & Growth model in which the burn rate is a function of the local pressure.We expect that the different behavior for the rate behind the lead shock would have a large effect on short-shock initiation experiments.
Further details on the SURF model and comparisons with additional experiments can be found in a companion paper of ours (Shaw and Menikoff 2010).

Final remarks
We have focused on shock initiation.But a burn model needs to deal with two additional detonation phenomenon: (i) The curvature effect for propagating detonation waves.(ii) Change of the ignition sensitivity with the initial state of a PBX.We conclude with a brief comments on incorporating these properties within the context of a hot-spot model.
The detonation speed and state of an underdriven planar detonation wave is uniquely determined by the CJ condition and the equation of state of the HE.In general, however, the detonation speed and state depends on the local curvature of the detonation front.This is known as the curvature effect.Due to the jump conditions for a quasi-steady detonation (Menikoff et al. 1996), the curvature effect is sensitive to the reaction zone width.This can lead numerical simulations to have a mesh dependent curvature effect.Typically, the reaction zone width decreases with a finer mesh and results in a smaller numerical curvature effect.
Alternatively, the reaction rate of a burn model can be adjusted to tune the reaction zone width in order to fit the curvature effect.In fact, Wescott et al. (2005) fit both regimes for PBX 9502 using a variation of the original Lee & Tarver model; terms were added to the rate which effectively switch from initiation parameters to propagation parameters based on the shock density; see (Wescott et al. 2005, Eqs. (44) and (45 We note that the reaction zone for a propagating detonation wave depends on the rate for the lead shock pressures greater than P CJ , while shock initiation depends on lead shock pressures less than P CJ .Thus for our formulation, it would be natural to fit the shock strength factor f (P s ) separately for the two disjoint regimes.Simulations would need the reaction zone to be sufficiently resolved.Very likely adaptive mesh refinement of the detonation front would be needed for efficiency.
Initiation sensitivity of a PBX depends on its initial state.Lower initial density increases the shock sensitivity.This is compatible with a higher porosity resulting in an increased number of potential sites for hot spots.In addition, a PBX is more sensitive when hot then when cold.This can be attributed to varying the number of hot spots that become burn centers.
Dissipative mechanisms that localize energy in a shocked PBX give rise to a temperature distribution; see for example Menikoff (2004).Hot spots would correspond to the high temperature tail of the distribution.Due to the temperature sensitivity of a chemical rate, burn centers (in which thermal runaway occurs on a fast time scale) can be approximated as those hot spots above a cutoff temperature.A higher initial temperature would shift the temperature distribution resulting in a greater density of burn centers.For our formulation of a hot-spot burn model, the change in sensitivity could be taken into account by modifying the shock strength factor to depend on the initial state as well as the shock pressure; i.e., f (P s , ρ 0 , T 0 ).
Finally, we note that all the burn models mentioned in sec. 2 have been calibrated to fit some experimental data.To assess the relative merits of different models one needs to compare each model with a large number of experiments covering a wide range of detonation phenomenon.We think that 00003-p.7 EPJ Web of Conferences a burn model that incorporates the underlying physics of hot spots is likely to do better than ad hoc fitting forms for the reaction rate.

Fig. 4 .
Fig. 4. Comparison of embedded velocity gauge data for shock initiation experiments in PBX 9502; shots 2s57 (13.5 GPa) and 2s68 (11.5 GPa) from Gustavsen et al. (2006).Solid lines are experimental data and dashed lines are simulations.Gague positions are shown in the legend of each plot.