AWSD reactive flow model for PBX 9404 Open Access

PBX 9404 is an HMX-based high explosive (HE) consisting of 94% HMX, 3% tris- $β$-chlorethylphosphate (CEF), 2.9% nitrocellulose (NC), and 0.1% diphenylamine (DPA) by weight. It is no longer in active use but still plays an important role in the development, assessment, and validation of theoretical models applied to legacy HE data. The results of those theoretical models can be used to make comparisons to other HMX-based HEs in order to understand the performance of different HE compositions.¹ PBX 9404 has a high energy density, so it is a powerful and brisant HE.² Nitrocellulose—the binder in PBX 9404—is chemically unstable, which raises questions about how the safety and performance of PBX 9404 change over long time-horizons, specifically years to decades. In the 1960s, PBX 9404 was one of the PBXs with the highest energy density that could be obtained with the state-of-the-art HE formulation and fabrication methods.² Over time, PBX 9404 has been replaced with better and safer HEs, such as PBX 9501 and PBX 9502. However, because of its historical significance in HE research and its usefulness in HE aging studies,³ examining the behavior of PBX 9404 by developing theoretical models to describe its behavior is currently an important research area.

There is a significant amount of experimental data available for PBX 9404, and these data can be used to calibrate and validate theoretical models. Svingala et al. measured in situ particle velocity wave profiles using gas-gun plate-impact experiments.¹ They specifically used an embedded magnetic particle velocity gauge technique to obtain particle velocity histories and shock wave times of arrival for PBX 9404.¹ Hugoniot data for PBX 9404 are available from multiple sources. Marsh,⁴ Gibbs and Popolato,⁵ and Svingala et al.¹ have all reported Hugoniot data for PBX 9404 in the reactants state. Green et al. and Kineke et al. have reported overdriven Hugoniot data for the PBX 9404 products state.^6,7 These experimental data can be supplemented and supported using thermochemical codes. For example, codes, such as MAGPIE, can be used to predict thermodynamic properties of PBX 9404 under various conditions.⁸ 

Models of shock-induced chemical reactions and high explosive detonations play an important role in predicting the performance of HEs and energetic materials.^9–16 For example, reactive flow models of HEs are used to understand and predict the performance of energetic materials and HEs under various shock conditions. One emerging reactive flow model is the Arrhenius–Wescott–Stewart–Davis (AWSD) model.^12,17 The AWSD model consists of a Davis reactants equation of state (EOS), a Davis products EOS,^18–20 and an Ignition and Growth type of model^21,22 that describes the time evolution of the reaction progress with modified rates having an Arrhenius form^12,16,23,24 with the temperature being the shock temperature. AWSD models have been successfully applied to model shock-to-detonation (SDT) behavior and other thermodynamic and mechanical properties in several HEs, including PBX 9501, LX-14, TNT, and others.^17,25–27 One of the principal limitations of the AWSD model is that, without modifications, it is unable to describe thermal nonequilibrium chemical kinetics.^28–33 Nevertheless, the AWSD model has been shown to be a powerful theoretical tool in HE modeling, often able to capture reactive flow of an HE with more accuracy than other reactive flow models (for example, compare the results in Ref. 26 with those in Refs. 34–36).

In this work, we present an AWSD reactive flow model for PBX 9404. Specifically, we calibrate an AWSD model to experimental data from a number of different sources and to calculated thermochemical results. After performing the calibration, we find that the AWSD model is in excellent agreement with previous experimental results, including data that were not used in the calibration procedure. Other computational models have been developed to predict the detonation propagation of PBX 9404.^37–40 These models can be used to examine the reactive behavior of the HE but do not provide the robust description of the detonation of PBX 9404 given here. The results presented in this work illustrate the ability of the AWSD model to accurately describe the reactive burn of PBX 9404. This work provides a model based on data from a legacy HE that can be used for comparison to future HMX-based HE formulations.

The remainder of this article is organized as follows: The mathematical details and results of the AWSD calibration are presented in Sec. II. Section II A contains the details of the reactants EOS calibration. Section II B contains the details of the products EOS calibration. In Sec. II C, the kinetic parameters in the AWSD model are calibrated to experimental gas-gun wavefront data and to experimental diameter effect results. In Sec. III, the performance of the model is validated by comparing to experimental gas-gun and cylinder expansion test (CYLEX) data. Conclusions and implications of this work are discussed in Sec. IV.

We follow the mathematical exposition of the AWSD model given in Refs. 12 and 17. The calibration procedure we use generally follows Ref. 25. For completeness, some of the exposition and explanation of the model given in Ref. 26 is repeated here.

The Davis reactants EOS has ten parameters: $A$⁠, $B$⁠, $C$⁠, $C v ( 0 )$⁠, $Γ 0$⁠, $Z$⁠, $α s t$⁠, $ρ 0$⁠, $T 0$⁠, and $E det$⁠. Values for these parameters can be obtained by calibrating the EOS model to a combination of experimental data, thermochemical calculations, and simulation data.

The parameter $A$ in the Davis EOS model is the ambient state sound speed of the HE. For PBX 9404, values for $A$ are available from several sources. The value given by Marsh is $A=2.26km/s$⁠.⁴ The value taken from Ramsay and Popolato is $A=2.43km/s$⁠.⁴² Here, we use a value of $A=2.18km/s$ as we found that this value resulted in a calibrated EOS that was in the best agreement with the experimental data; specifically, the SDT results presented later were better fit using this value.

The parameter $B$ determines the slope of the Hugoniot in the shock velocity (⁠ $U S$⁠)–particle velocity (⁠ $U P$⁠) plane at small $U P$⁠. Here, we used the value of $B=3.3$⁠, which is the AWSD value for PBX 9501 given by Aslam et al.²⁵ We used this value because we found that calibrating $B$ to the available $U S- U P$ data for PBX 9404 resulted in highly stiff EOS prediction, exceeding the stiffness of HMX.

The $C$ parameter was calibrated by applying the Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimization algorithm⁴⁷ to minimize the error with both the Hugoniot data from Ref. 1 and single crystal HMX data from Ref. 4. In the optimization, the HMX $U S$ data were shifted down by $150m/s$⁠. The rationale for this shift is that we estimate the porosity of PBX 9404 to be $≈2%$ through a comparison of the nominal density $ρ 0=1.841g/ cm 3$ and the theoretical maximum density $ρ max=1.878g/ cm 3$⁠. The initial guess for $C$ was taken to be the value for PBX 9501 from Ref. 25. This procedure resulted in a value $C=0.5329$ for PBX 9404. Compare this value to values of $C=0.2$ for PBX 9501 and $C=4.591$ for LX-14.^25,26 A full listing of the calibrated Davis reactants EOS parameters for PBX 9404 is given in Table I. The value for $E det$ is determined from the calibrated products EOS described later.

The results of the reactants EOS calibration are shown in Fig. 2. The $U S- U P$ results are shown in Fig. 2(a). The calibrated PBX 9404 EOS is shown as a solid blue curve. The circular markers are data points to multiple data sources. The calibrated EOS is in strong agreement with the experimental data across all $U P$ values. The dashed black curve is the PBX 9501, and the dashed orange curve is the result for HMX. At larger $U P$⁠, the PBX 9404 EOS value lies between HMX and PBX 9501 because the PBX 9404 EOS is more stiff than PBX 9501 but not as stiff as HMX due to factors, such as porosity.⁴⁸ The PBX 9501 and PBX 9404 Hugoniot curves are similar in magnitude and functional behavior due to their close chemical relationship. The Hugoniot curves in the $P-V$ plane, where $P$ is the pressure and $V$ is the specific volume, are shown in Fig. 2(b). As the specific volume becomes smaller, the PBX 9404 pressure exceeds the PBX 9501 pressure, a result of the greater stiffness of the PBX 9404 EOS in comparison with PBX 9501. For specific volume values larger than $≈0.4 cm 3/g$⁠, the PBX 9501 pressure is slightly greater in relative terms than the PBX 9404 pressure.

There are seven parameters $a$⁠, $b$⁠, $k$⁠, $n$⁠, $p c$⁠, $v c$⁠, and $C v p$ in the Davis products EOS model (Table II). The values of these parameters are generally determined by fitting the EOS model to available data sources. Here, to fit the products EOS, we used overdriven Hugoniot data from Green et al. and Kineke et al. for the PBX 9404 product state supplemented with data generated from thermochemical calculations performed using the code MAGPIE.^6–8 

The results of calibration are shown in Fig. 3. Figure 3(a) shows the products Hugoniot curve in the $U S- U P$ plane. The circular markers in the figure are the experimental results, and the blue solid curve is the Davis EOS optimized to fit the experimental data and the calculated values. Good agreement is observed between the calibrated EOS and the experimental results. The results of the PBX 9501 Davis products EOS, shown as a dashed black curve, are shown for comparison to PBX 9404. The PBX 9404 EOS predicts a slightly higher $U S$ value than PBX 9501 as $U P$ becomes larger. The PBX 9404 CJ state given by Gibbs is shown as a square orange marker, and the corresponding state calculated using the calibrated Davis products EOS is shown as a square blue marker. The detonation velocity from Gibbs is $D CJ Gibbs$ $=8.803km/s$⁠, and the optimized value is $D CJ=8.802km/s$⁠;⁵ therefore, our model is in excellent agreement with the known detonation velocity. Figure 3(b) shows the products Hugoniot curves in the $P-V$ plane, again with good agreement observed between the calibrated EOS and the experimental data. The value for the pressure at the CJ state given by Gibbs is $P CJ Gibbs=36.8GPa$⁠, and the calibrated value is $P CJ=34.9GPa$⁠.

The Hugoniot $U S- U P$ curves for the PBX 9404 products EOS and reactants EOS are both shown in Fig. 4. This plot illustrates that Hugoniots from the two calibrated equations of state do not cross. Similarly, although not shown, we have confirmed that the two EOSs do not cross in the $P-V$ plane. This is an important observation because the detonation products EOS being above the reactants EOS may avoid computational problems when the EOS models are used in multiscale simulations and hydrocodes due to the unique identification of the system thermodynamic state for each EOS.³⁶ 

Svingala et al. performed gas-gun experiments to measure the wave front propagation of PBX 9404 at various pressures.¹ They specifically performed embedded particle velocity gauge measurements of PBX 9404 at different pressures using gas-gun plate-impact techniques. Five different shots were performed, and a list of those shots is given in Table III. Gibbs and Popolato reported diameter effect data for PBX 9404.⁵ We use these datasets to calibrate the AWSD reaction rate parameters.

We calibrated the rate law parameters in the AWSD model by fitting the results of Eq. (20) to four of the shots of Svingala et al., specifically 1S-1665, 1S-1666, 1S-1668, and 2S-1079. The data from shot 1S-1667 were not used in the calibration so that it could be used to validate the developed model. The mesh size used in the simulations was $10μm$⁠, and the simulations were performed using the multimaterial multiphysics hydrocode FLAG.^50,51 The gas-gun impactors were either Sapphire or Kel-F81, depending on the shot as shown in Table III. Both materials were modeled with a linear $U s$- $U p$ EOS, with the Sapphire parameters taken from Ref. 52 and the Kel-F parameters taken from Ref. 53.

There are several steps in the rate law parameter calibration. They are listed sequentially:

1. A Nelder–Mead simplex optimization algorithm⁵⁴ was used to obtain an initial set of optimized parameter values. This procedure used the experimental gas-gun data from shots 1S-1665, 1S-1666, 1S-1668, and 2S-1079 as target data. The objective function in this calibration was the sum of the time-averaged absolute difference between the experiment and the simulation for each embedded gauge signal. The simplex optimization procedure was employed to obtain the parameter values that minimize this objective function. More specifically, the Nelder–Mead optimization was performed starting with the initial guess for the parameter values being the AWSD parameter values for PBX 9501.²⁵ 400 iterations of the simplex optimization were performed, and we found that this resulted in a converged value for the model parameters. The error metric used was the $L 1$ error, which was calculated at each time point for each gauge and over each shot by comparing the experimental and simulation data. Every time point in this optimization received equal weight in the error metric. The error value at each point for each gauge in the selected shots was then combined into a single total error metric, which was minimized using the Nelder–Mead procedure.

2. Rate-stick experiments for PBX 9404 were simulated using a shock-fitting code.^55,56 We used diameter effect data taken from Ref. 5 to perform the calibration. The parameters $k 2$ and $b 2$ were calibrated in this step, while all other parameters were held constant. A variational Bayesian calibration method was used in this step.⁵⁷ In this calibration, the objective function was used to maximize the posterior probability of the AWSD model parameters given the diameter effect data. The likelihood of the diameter effect data was modeled as a normal distribution with a variance given by the experimentally measured uncertainties in detonation speed.

3. A Nelder–Mead optimization procedure for the gas gun data was performed again to obtain an updated set of parameter values calibrated to experimental results. The overall procedure was similar to step 1.

4. We hand-tuned the ambient sound speed (the parameter $A$ in the reactants EOS) in order to get better agreement with the SDT experimental results. After setting $A$⁠, the parameter $C$ in the reactants EOS was again calibrated to experimental EOS data using a BFGS optimization algorithm.

5. Another Nelder–Mead optimization procedure over the gas-gun data was performed to obtain an updated set of AWSD parameter values using a similar procedure as in steps 1 and 3.

6. The parameters $k 2$ and $b 2$ in $F 2$ (the high temperature rate) were recalibrated to the available diameter effect data.

Figure 5 shows a comparison between the SDT results predicted by the calibrated PBX 9404 AWSD model and experimental data from the four shots used in the calibration. The principal overarching observation is that the AWSD model is in good overall agreement with the experimental results for each shot. Before the fitting was performed, we trimmed noisy tails of some gauges, for example, gauges 4 and 5 in shot 1S-1665 and gauge 7 in shot 1S-1666 so that the AWSD optimization would better target the initial wavefront shape of the experimental results. At later times, gauges can break and or bend yielding non-physical velocities. Gauge 8 in shot 2S-1079, the cyan curve from position $5.89mm$ in Fig. 5(d), had a significant fallout, so the data from that gauge were not used in the calibration.

In Fig. 5(a), the results for shot 1S-1665 are shown. The initial pressure of 1S-1665 was $5.27GPa$⁠. For this shot, the results of the calibrated model are in good overall agreement with the experimental data. The shock arrival times, maximum velocities, and wavefront shapes are predicted well by the model. For gauges further from the shock initiation, the model predicts later arrival times than the experimental values. Overall, the calibrated model performs well on this shot.

The results for shot 1S-1666, which had an initial pressure of $4.18GPa$⁠, are shown in Fig. 5(b). The wavefront shapes, initial shock front particle velocities, shock arrival times, and maximum velocity magnitudes are all in very strong agreement with the data for each gauge. Overall, this is the shot on which the AWSD model performed the best.

The SDT results for shot 1S-1668 are shown in Fig. 5(c). This was the shot with the lowest initial pressure (⁠ $2.98GPa$⁠) in the suite of PBX 9404 shots of Svingala et al. The AWSD model and other reactive flow calibration models often do not accurately capture the shock wavefront in low pressure regimes. However, in this case, the wavefront shapes and maximum velocities are in good agreement with the data. We observe excellent agreement with the wavefront shapes for each gauge. The shock arrival times predicted by the AWSD model are slightly early for all gauges.

The results for the two-stage shot 2S-1079 are shown in Fig. 5(d). The initial pressure of this shot was $6.66GPa$⁠. Again, the AWSD model is in strong agreement with the experimental results. The wavefront shapes are generated well by the model. The AWSD shock arrival times are also in good agreement with the experimental results, particularly further from the shock initiation.

Figure 6 shows the diameter effect results for PBX 9404. The red markers are the results of the AWSD model, and the solid red curve is a fit to that data that we include to guide the eye. The green triangular markers are experimentally measured values for PBX 9404. For large diameters (small inverse radii), the AWSD model captures the experimental behavior of the diameter effect well. This agreement is continued for small diameters (large inverse radii); however, visible deviations between the experimental datapoint and the model are observed for the experimental points with $D 0≈7.9mm/μs$ and $D 0≈7.3mm/μs$⁠. Note that compared to simulations by Tarver in Ref. 59, the AWSD model is generally in better agreement with experiments for small inverse radii and worse agreement for inverse radii above approximately $1.5 mm − 1$⁠. Overall, the experimental PBX 9404 rate-stick data are captured well by the AWSD model. For comparison, experimental data for PBX 9501 are shown using gray triangular markers. The experimental velocities for PBX 9404 are generally higher than PBX 9501, and that behavior is captured in the AWSD model.

The final parameter values used in the AWSD rate law are given in Table IV.

We did not use the data from shot 1S-1667 in the model calibration so that it could be used for model validation. The results for this shot are shown in Fig. 7. The reasons we withheld this shot are (a) it is an intermediate pressure regime relative to the other shots and (b) the first two gauges failed in this experiment, so we concluded that it should be used for validation as opposed to calibration. The AWSD results and the experimental data are in overall good agreement. The predicted wavefront shapes are generally very close to the experimental shapes, and the maximum velocities are in reasonable agreement with the experimental values. The shock arrival times are excellent early in the shot, but the AWSD model predicts slightly later arrival times than the experimental results further from the shock initiation. Overall, these validation results illustrate that the calibrated AWSD model can be used to predict the SDT performance of PBX 9404 in regimes and conditions that were not included in the data used to calibrate parameters in the model.

The developed model was also validated against CYLEX data. The standard CYLEX tests consist of a nominally 1.0-in. diameter by a 12-in. long explosive charge surrounded by a 0.1-in. thick copper tube and initiated at one end. The expanding copper wall velocity is measured with Photon Doppler Velocimetry (PDV) probes after a steady detonation wave has developed. Simulation of a 1.0-in. CYLEX experiment for PBX 9404 was run in FLAG using the calibrated AWSD parameters. A short 0.5-in. long section of programmed burn PBX 9404 was used to initiate the AWSD reactive burn HE region. The simulations were 2D axially symmetric with an initial zone size of 40  $μ$m. The copper was modeled with a Sesame EOS and PTW strength model.⁶⁰ Wall velocities in the simulations are compared with PDV probe data from CYLEX experiments with PBX 9501⁶¹ and PBX 9404³ in Fig. 8. PBX 9501 and PBX 9404 are both HMX-based HEs with similar detonation performance, so wall velocities are similar between the experiments. However, there is an ongoing debate about the effect of aging on PBX 9404 performance. In these recent PBX 9404 CYLEX tests, the HE was over 46 years old and is believed to have a slight decrease in detonation energy when compared with less aged material.³ Overall, the AWSD calibration provides a good agreement with the experiment and is only slightly below the PDV data during the early wall expansion.

We have calibrated an AWSD reactive flow model for the high explosive PBX 9404. The calibrated model is able to accurately predict the wavefront shape and shock arrival times of gas-gun experimental data and is in strong agreement with the experimental diameter effect and CYLEX results. The model’s performance was validated on experimental data that were not used to calibrate parameter values, with good agreement observed between that unseen experimental data and the model prediction. The primary implications of this work are that it can be used to (a) compare the performance of a legacy HE (PBX 9404) to other HEs and (b) examine how aging in PBX 9404 affects performance. Specifically, in future work, the AWSD model for PBX 9404 developed here can be compared against reactive flow models for other HMX-based high explosives in order to further understand the performance of different HE compositions. It can also be compared to AWSD models developed for aged PBX 9404 in order to determine the specific ramifications and implications of aging.

This work was supported by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218NCA000001).

The authors have no conflicts to disclose.

Galen T. Craven: Investigation (equal). Matthew A. Price: Investigation (equal). Stephen A. Andrews: Investigation (equal). Kirill A. Velizhanin: Investigation (equal). Tariq D. Aslam: Investigation (equal). Jeffery A. Leiding: Investigation (equal). Christopher Ticknor: Investigation (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request and institutional approval.