Uncertainty quantification of material parameters in modeling coupled metal and high explosive experiments

Accurate calibration and assessment of associated uncertainty for equation of state (EOS) models related to energetic material detonation is of significant interest to the defense community for a wide variety of applications. The main challenge of calibrating these models emanates from complex physical chemistry underlying the chemical phase change occurring during the detonation process, particularly in high explosives (HE).¹ Performing focused experiments often requires a coupling between the HE and a metal; here, the metal facilitates experimental measurement of the free surface velocity, which aids in the modeling of the HE.^1–3 When modeling such experiments, this coupling between the two material types can potentially introduce additional uncertainty given the complex loading conditions imparted on the metal and the underlying complexity of accurately modeling metal strength. The constitutive model related to the metal strength will often require formulations that include rate, temperature, and pressure dependence as well as porosity nucleation and growth.

Our goal is to quantify the response uncertainty of the coupled metal/HE experiments. Of particular interest is the degree to which uncertainty in the metal strength model affects the uncertainty of the experimental response as well as the inferred HE parameters. We employ Monte Carlo error propagation to assess the impact of metal strength uncertainty on the experimental response uncertainty, and we use sensitivity analysis to identify the most significant parameters in the models of coupled experiments. Both the metal strength parameters and the HE parameters of such models are estimated using Bayesian inference, which naturally quantifies the associated uncertainty.

We consider two types of coupled metal/HE experiments: small plate tests and cylinder tests. We consider small plate tests using either LX-14³ or LX-17 HE² as well as a cylinder test using LX-17 HE.¹ LX-14 and LX-17 are two commonly used HE materials, and analyzing experiments with both types of HE provides a comparison between sensitive and insensitive HE, respectively. To leverage previous work in the statistical modeling of material strength, we focus on using tantalum as the metal component in each of these experiments. While tantalum has been well studied, modeling it across a wide range of pressure, temperature, strain rate, and strain conditions has proven to be quite challenging.^4–6 

Recently, Schill et al.⁷ studied tantalum strength using Bayesian analysis and compared pulsed power and gas gun experiments.

For a discussion of specific physical effects that accompany release of material from the compressed state, see Refs. 8, 9, 10, and 11.

For a more general application of calibration methods in problems of dynamic compression of condensed matter, we suggest Refs. 12–20. More broadly, similar methodology has been applied in high energy density physics^21–27 as well as computational plasticity.^28,29

Sensitivity analysis of experiments involving tantalum have also been previously studied.^7,30,31 References 7 and 31 specifically consider sensitivity of tantalum strength parameters. As shown in Ref. 31, the strength parameter sensitivity differs by experimental regime.

The literature available on HE model calibration and uncertainty using statistical inference is limited, but it does provide some useful resources. Higdon et al.³² calibrated energy-related parameters for an HE of interest in a cylinder expansion test using Bayesian inference. They found that the determination of the HE model parameters was challenging with limited experiments. In Lee et al.,³³ they calibrated EOS models for a variety of HEs to rate-stick experiments. In Andrews and Fraser,³⁴ they employed a Bayesian approach to verify that an HE EOS could be recovered using synthetic data. This work provides a basis to assess the effects of a metal constitutive model on our ability to infer HE model parameters. Additionally, Andrews et al.³⁵ applied Bayesian calibration to the Davis Products EOS model for PBX 9501, while Lindquist et al.³⁶ calibrated the reactant and product forms of the Davis EOS to a variety of conventional HE materials.

We leverage this existing work to analyze the uncertainty of coupled metal/HE experiments and unlike previous studies we include uncertainty of the metal strength model parameters. The outline of this study is as follows. In Sec. II, we describe the metal/HE experimental configurations. The modeling of these experiments is discussed in Sec. III. Section IV gives an overview of the statistical methodology. Implementation of these methods and their results are presented in Sec. V. In Sec. VI, we provide conclusions, including experimental recommendations based on our analysis.

In this section, we describe the small plate and cylinder experimental configurations. In both cases, tantalum is chosen for the metal component. Data collected from the experimental setups described below can be found in Refs. 1–3.

The first experimental configuration considered in this investigation is the small plate test. In this experiment, a thin disk of high explosive material is pressed against a thin plate of metal (see Fig. 1). We modeled two types of HE: sensitive (LX-14) and insensitive (LX-17). In either case, the detonation of the chosen HE pushes a tantalum plate, and the free surface velocity is measured using Fabry–Pérot interferometry.

The cylinder test is a calibrated pipe bomb designed to give detonation energies by measuring the free-surface velocity. The HE is ignited at one end of the cylinder with the velocity measurement taking place at approximately 2/3 of the total length assuming steady detonation has been achieved. In this case, the velocimetry is performed using Fabry–Pérot interferometry. The expansion process of the metal cylinder is visualized in Fig. 2.

A summary of the experiments and the salient details can be found in Table I. The small plate experiments were chosen due to the relatively thick tantalum plates, which increases the likelihood that spall will convolute the ability to infer HE model parameters. Additionally, the small plate configuration are experiments with large HE and metal diameters, which allows us to neglect 2D effects until relatively late times. Lastly, the tests are historical in the sense they are experiments that were performed using methods and data reporting that was not as detailed as it would be with modern experimental methods. The following analysis will allow us to assess the efficacy of historical experiments when calibrating modern HE models by assessing the effects of variance on as-reported measurements, notably the charge thickness in plate experiments.

Due to the complexity of coupled metal/HE experiments, modeling is achieved via multiphysics simulations. In this section, we highlight components of these simulations. Constitutive models for the metal component are detailed in Sec. III A. A brief description of the HE model is given in Sec. III B. The simulations were performed using the LLNL hydrocode, Ares.³⁷ The small plate and cylinder experiments were modeled with the metal being Lagrangian and the high explosive leveraging the code’s arbitrary Lagrangian–Eulerian (ALE) capability to allow the mesh to relax at late times. The configuration was modeled as 1D. The cylinder test was modeled with a 2D axisymmetric representation given the rotational symmetry about the center axis of the cylinder. A large number of simulations were run in order to train the Gaussian process used to emulate the computationally intensive multi-physics code. See Sec. IV A for more details related to the emulator.

The experimental configurations simulated here involve HE pushing metal in two different configurations. This requires that each material has its own constitutive model with the HE model simulating the detonation process related to a phase and thermochemical change from an unreacted solid to a gaseous product and condensed detonation soot. The metal, on the other hand, is modeled with what is colloquially known as a strength model. Here, the strength model comprises a combination of an equation of state (EOS), flow stress model, and a porosity model for smplt-LX-14 and cyl-LX-17 experiments and only the Preston–Tonks–Wallace (PTW) flow stress model for smplt-LX-17 due to challenges with coupling a porosity model and the soot parameter capability. The PTW³⁸ model, which is commonly used in hydrocodes, serves as our flow stress model. When porosity is included, the PTW model is coupled to porosity models, which help account for the wide range of strain rates, pressure, and temperature conditions experienced by the metal. For the following simulations, two different porosity models are used; they are described in Sec. III A 2. The EOS for the metal is Mie–Grüneisen in form and believed to be well understood and calibrated in the range of conditions in this study, therefore, it is not included in the calibration.

A useful feature of Bayesian statistical approaches is that they make use of subject matter expert knowledge. One way that knowledge comes into play is via selection of suitable ranges for parameter values. The PTW model parameters and their expert-determined bounds for tantalum are given in Table II. Parameters that are not specified in the table are fixed at their nominal values from Ref. 38 rather than inferred. Note that $γ$ is sampled on the log (base 10) scale when ensembles of simulation inputs are employed to cover the admissible parameter space for emulator training.

This porosity model is coupled to the PTW flow stress model (see Sec. III A 1), which also includes material softening due to the increase in porosity.⁴⁴ A subset of these porosity parameters are inputs to the coupled metal/HE experimental simulations; these parameters are noted in Table III. The following additional parameters were included in the emulator training simulations and calibration: volume fraction available for stress nucleation $f n α 0$⁠, calibration parameter for the strength degradation function $a 1$⁠, and the average pressure at which pores will nucleate $σ h M α$⁠. A more complete discussion of this porosity model along with the nucleation parameters can found in Qamar et al..⁴⁵ For evaluating mesh resolution effects, further work in the area could involve V&V based approaches.⁴⁶ 

For the HE model, we used the Cheetah thermochemical code to model the detonation of the two HE materials. Cheetah controls the conversion of unreacted HE materials to gaseous product species through kinetic rate laws as well as solve sets of thermodynamics equations to maintain instantaneous chemical equilibrium between possible product species. This allows us to estimate the pressure–volume relationship necessary for the HE EOS. For a more in-depth discussion of Cheetah, look to Ref. 48.

In Table V, we present the Cheetah model parameters and bounds used in the emulator training. The different experimental configurations modeled used different combinations of parameters. Choice of these parameters is detailed below.

For the LX-14 small plate configuration discussed in Sec. V A, we utilized two user adjustable parameters recommended for fine tuning the HE models, namely, the energy shift (⁠ $E shift)$ and chemistry freeze-out pressure (⁠ $P freeze)$ parameters. The energy shift (⁠ $E shift)$ parameter will adjust the energy density by specified amounts. The chemistry freeze-out pressure (⁠ $P freeze)$ represents pressure below which the Cheetah assumes no further instantaneous chemical equilibrium. Additionally, the charge thickness (⁠ $chrgthk)$—or, the thickness of the HE disk—is used as a parameter with appropriate bounds because the reported value was an average of multiple measurements on a nonuniform disk. For the LX-17 small plate configuration discussed in Sec. V B, an additional parameter controlling the molar density of soot formed through detonation⁴⁹ was included, which is denoted as $V m$⁠. Given that the LX-17 small plate data are also from historical experiments, we opt to include a charge thickness parameter for the same reasons as in the LX-14 small plate calibration. Lastly, we model the LX-17 cylinder test discussed in Sec. V C with energy shift (⁠ $E shift)$ and chemistry freeze-out pressure (⁠ $P freeze)$ parameters.

Here, we introduce the statistical methodology applied to the models given in Sec. III. In particular, we describe how a Gaussian process can be used to emulate the complex multi-physics models used to characterize our metal/HE experiments. We also highlight sensitivity analysis, which identifies the most influential model parameters, and Bayesian inference, which is used to calibrate model parameters and quantify the associated uncertainty.

Recall that our integrated metal/HE experiments are modeled with computationally expensive multi-physics models. For many statistical methods, a large number of model evaluations are required. With such an expensive model, this is not computationally feasible. In these cases, we instead utilize an emulator—which closely approximates an expensive model while also being much cheaper to evaluate—in place of the original model.

In particular, we employ a Gaussian process regression (GPR) model as our emulator for each of our three experimental configurations (small plate tests with both LX-14 and LX-17 and a cylinder test with LX-17). GPR is an attractive choice for an emulator due to the built-in quantification of uncertainty introduced due to emulator approximation. Each Gaussian process was implemented using the GPy package.⁵⁰ We chose a Matéernel with $ν=5/2$⁠.

To train each Gaussian process emulator, we performed a large number of simulations (see Table VI for more details). These simulations were run using sets of parameters (appropriate to the particular experiment) obtained via a Latin hypercube design bounded by the ranges given in Sec. III. These simulations were divided into training and testing sets. The training set, consisting of parameters and simulation outputs, was used to train the Gaussian process. The testing set was used to assess generalization beyond the training data. The simulations conducted for emulator training consumed a non-trivial amount of computational resources. The 1D LX-14 simulations utilized one node (with 36 cores) for 3 min of wall-clock run time for a total of 60 k core-hours, and the 1D LX-17 simulations used a similar allocation. The more computationally intensive 2D cylinder calculations utilized 168 cores with a run time of $∼2$ h each, for a total of $∼1.4× 10 6$ core-hours. In Fig. 3, the behavior of the emulator is compared to the original simulation for the LX-14 small plate configuration. This experimental condition was chosen as an exemplar for assessing emulator performance. The desired behavior in such a plot is the line of slope one passing through the origin, which would indicate a perfect match between the model and the emulator. The tighter the grouping around this line, the better the emulator is at mimicking the behavior of the original model. Note in Fig. 3 that the training set and the testing set exhibit similar grouping, indicating that we have not overfit. The gray region indicates the estimated experimental error from Fabry–Pérot interferometery.

The number of simulations necessary for a well-trained emulator was first determine for the LX-14 small plate experiment since had the highest dimensional parameter space and then replicated for the LX-17 small plate configuration. The size of the training test was determined by progressively increasing the number of simulations until the emulator training and testing performance was acceptable. Due to the similarities in the velocimetry, we are emulating the dimensionality was the dominate consideration when determining the size of the training set. A reduced training set was used for the cylinder test configuration given the previously noted computationally intensive calculations. Given the emulator performance was equivalent to the other configurations the number of training calculations was determined to be sufficient.

In sensitivity analysis, the variation (or, uncertainty) in a model response is proportioned to variation in each of the input variables (i.e., parameters). The goal of this type of analysis is to ascertain which parameters have the greatest impact on the model response. Sensitivity analysis techniques can be categorized as either local or global. Local sensitivity analysis quantifies the change in model output when the parameters are perturbed about a nominal value—e.g., using partial derivatives evaluated at a particular point. Global sensitivity analysis characterizes parameter influence throughout the entire range of possible parameter values (i.e., the admissible parameter space). We employ Sobol’ indices, a measure of global sensitivity, for our analysis. See Ref. 19 for more information of Sobol’ indices and their computation.

Standard Sobol’ indices, as described in Ref. 19, are used for models with scalar responses. For each of our experimental conditions, our model response is functional—namely, measured velocity over time. We instead employ a functional version of Sobol’ indices as in Refs. 30 and 31. This allows us to assess how the parameter sensitivity evolves over time. Note that the calculation of Sobol’ indices requires a large number of model evaluations, so we use an emulator, as described in Sec. IV A, for computational efficiency.

In addition to the experimental data, information encoded in the prior distribution $π 0$ also affects the parameter inference from (8). The prior distribution is chosen in such a way to reflect any previous knowledge about the model parameters (e.g., physical constraints, expert knowledge, or previous experiments). Note that information from previous experiments can only be used in the prior if you do not reuse the previous data in the current calibration (via the likelihood function).

Experimental data are incorporated in Bayesian inference through the likelihood function, denoted as $L$ in (8). Here, $L( y| θ)$ quantifies the likelihood of observing measurements $y$ given model parameters $θ$⁠. An appropriate likelihood function is derived from a combination of the physics model and an error model for the statistical noise. See Ref. 19 for more details.

As stated above, the goal of Bayesian calibration is to obtain the posterior density defined by (8). Despite the seemingly straightforward formulation, computation of the posterior is often difficult in practice. Due to the normalizing integral in the denominator, an analytic form of the posterior exists only in special cases.¹⁹ We instead utilize Markov chain Monte Carlo (MCMC) to obtain a posterior distribution. Since we cannot analytically represent $π( θ| y)$⁠, Markov chain Monte Carlo (MCMC) is employed to generate a sequence of samples from the posterior distribution. In particular, we use pyMC 2.3.8 to facilitate adaptive Metropolis MCMC sampling. The PTW model includes constraints on parameters as defined in Ref. 46. For the Bayesian calibration, constraints were enforced in the likelihood function. In practice, this was done in the code by assigning very small likelihood outputs to parameter inputs violating the constraints. This will lead to them being rejected during the acceptance/rejection stage of the Metropolis algorithm and not being included in the posterior.

To perform Bayesian calibration via MCMC, a large number of model evaluations are required. The exact number required is problem dependent, and in this case required $1.2 6$ model calls. If this were performed through hydrocode simulations it would take over 60 000 wall-clock hours, or $∼7$ years, for the small plate and over $2× 10 6$ wall-clock hours, or 273 years, for the cylinder tests. By leveraging the emulator, we can perform the same calibration in less than 2 h. Thus, we employ a Gaussian process emulator as described in Section IV A.

While Sec. IV introduced the statistical concepts, here, we give details on the implementation of our analysis. We first train a Gaussian process emulator as described in Sec. IV A. Next, we employ this emulator for Monte Carlo error propagation to ascertain how the metal flow stress affects the experimental response. We then perform sensitivity analysis followed by Bayesian inference.

Monte Carlo error propagation allows us to qualitatively assess the effect of parameter uncertainty on the uncertainty of the model response. We employ a distribution, corresponding to model input uncertainty, from which we draw parameters. We draw a large number of parameter values and use each of these resulting parameter sets to evaluate the model, giving an ensemble of model predictions. Due to the computational expense, we utilize our Gaussian process emulator. The width of the ensemble spread characterizes the model uncertainty attributed to the input uncertainty.

To assess the impact of well-characterized, the PTW strength model parameters on the model response uncertainty, we performed Monte Carlo error propagation with all three experimental configurations. By well-characterized, we mean a robust calibration to low-cost experiments such as quasi-static and Kolsky bar experiments as well as possibly an additional high strain rate, large strain experiment such as Taylor cylinder test. Focusing on the flow stress component of the model, we fixed all other parameters to the mean values of the samples used for training the Gaussian process emulator. For all three metal/HE experiments, we employed 500 parameters draws from two different distributions for the eight PTW parameters: a truncated uniform distribution (with each parameter sampled independently) bounded by the ranges in Table II and a posterior distribution obtained from a previous calibration of the PTW model to mechanical testing data. The calibration data set included stress-strain measurements from Ref. 51 (corrected in Ref. 52) and Taylor cylinder impact data from Ref. 52. Note that this posterior distribution is narrower than the bounded flat distribution; moreover, the parameters drawn from the posterior are not necessarily independent.

After ascertaining the qualitative influence of the metal flow stress parameters, we employed sensitivity analysis to quantify the relative influence of the full parameter set in all three experimental models. As introduced in Sec. IV B, we used Sobol’ indices to perform sensitivity analysis on models with functional output. Numerical calculation of the Sobol’ indices was done via Saltelli’s method (see Ref. 19 for more information). This technique involves a large number of model evaluations [in our case, $10 6( n p+2)$⁠, where $n p$ is the number of model parameters], so we again utilized Gaussian process emulation. We calculated first order Sobol’ indices over a dense grid of time points to obtain a picture of how the sensitivity changes over the time window of each experiment. While total sensitivity is often reported to provide a complete picture of parameter sensitivity, we found that sensitivity to higher-order parameter interactions was small for all experimental configurations. Hence, we limit our analysis to first order indices.

Finally, we performed Bayesian calibration for each experimental configuration. In each case, the model used for calibration is a Gaussian process emulator for the computationally expensive multi-physics code, and the likelihood function was derived using this model coupled with an error model of independent and identically distributed (iid) Gaussian noise. For all three setups, we performed two Bayesian calibrations, each with the same data and likelihood function but using two different prior distributions for the flow stress parameters. First, we employed a bounded flat prior. All PTW parameters were assigned a truncated uniform distribution (independent of all other parameters) reflecting the range of admissible values. These bounds are given in Table II. For our second prior, we used a collection of truncated normal distributions for the PTW parameters. The truncated normal distributions also used the same bounds as in the flat prior case. The mean and variance of the truncated normal distributions for the flow stress parameters were chosen to reflect information from a previous Bayesian calibration, which utilized the tantalum quasi-static, Kolsky bar, and Taylor cylinder data referenced above. The truncated normal distributions implemented in pyMC are such that they integrate to one. Priors pertaining to the non-PTW parameters did not change over the two calibrations.

Note that sensitivity analysis is often used to down-select the set of calibrated parameters. Insensitive parameters can be fixed rather than estimated with minimal effect to the calibration results. However, we opt to calibrate the full set of parameters in all cases, allowing us to use the resulting Bayesian posteriors to validate the sensitivity analysis results. In particular, the marginal posteriors for insensitive parameters will not be updated by the data and will appear very similar to the marginal prior.

We now present the results of our uncertainty analysis. Each of the following subsections pertains to a specific experimental configuration. In each subsection, we provide results of Monte Carlo error propagation, sensitivity analysis, and Bayesian calibration.

The first experimental configuration is a small plate test with LX-14 pushing a tantalum plate.³ The tantalum strength model is a combination of the PTW for the flow stress and the NP porosity model for the yield surface as described in Secs. III A 1 and III A 2, respectively. The HE calibration parameters include the initial energy shift (denoted as $E shift)$⁠, the freeze-out pressure for reactions (⁠ $P freeze)$⁠, and the charge thickness (⁠ $chrgthk)$⁠.

Figure 4 shows a comparison of the error propagation obtained with samples from the two different distributions. The response ensemble resulting from the bounded uniform distribution exhibits a wider spread, which corresponds to larger output uncertainty. The uncertainty resulting from the posterior draws is notably smaller with the spread of the responses visually approaching a single curve.

Figure 5 visually presents the first order Sobol’ indices as a function of time; the wider the colored band, the more sensitive the corresponding parameter. Recall that the Sobol’ indices partition the output variance to the variance of the parameters. Thus, the first order indices sum to one if there is no sensitivity to higher-order parameter interactions (e.g., the joint influence of two or more parameters).

For this experimental configuration, the HE parameter shift is obviously the most sensitive parameter. Other parameters with a small but visually observable level of sensitivity include PTW parameters $y 0$⁠, $−log⁡(γ)$⁠, $κ$⁠, and $s 0$⁠. The HE charge thickness parameter is similarly sensitive. The white dips into the plot near 5.0 and 5.4  $μ$s indicate some sensitivity to higher-order parameter interactions. In fact, of the three experimental configurations, the LX-14 small plate model showed the highest sensitivity to multiple-parameter interactions.

Recall that we are calibrating 16 parameters for the LX-14 small plate model: eight flow stress parameters and five porosity parameters as well as HE parameters $E shift$⁠, $P \,freeze$⁠, and charge thickness. We performed two Bayesian calibrations. First with a bounded flat prior for the PTW parameters and then with a truncated normal PTW prior. Priors for all other parameters were the same for both cases. With the exception of charge thickness, all parameters were assigned a truncated uniform prior (independent from all other parameters) reflecting the range of admissible values. Porosity parameter bounds are given in Table III. For the HE parameters, the bounds were defined for the energy shift and the freeze pressure to provide adequate coverage of the experimental data. Since the charge thickness is a well-understood value recorded at the time of the experiment, we chose a truncated normal prior with the mean equal to the reported value and with truncation bounds given by $±10%$ of the nominal value to account for any potential measurement or transcription errors; this prior gives more weight to the recorded charge thickness as compared to a relatively non-informative flat prior. We employed these two priors to determine if reducing the uncertainty in the flow stress parameters affects the inference of the other parameters. We now compare the results of the two calibrations. The resulting posteriors are plotted together in Fig. 17. Due to the size of the plot, these results were placed in  Appendix. When calibration is performed with flat priors for the PTW parameters, the resulting marginal flow stress posteriors, with the exception of those corresponding to $κ$ and $γ$⁠, remain relatively uninformed by the data—i.e., the posterior is very similar to the prior, which in this case is flat. This result agrees with the sensitivity analysis shown in Fig. 5. Similar to the majority of the flow stress parameters, the marginal posteriors for NP porosity remain largely unchanged from their priors. For the HE parameters, $E shift$—as expected from the sensitivity analysis—shows the greatest reduction in uncertainty when comparing the prior to the posterior. Calibration with the second prior results in narrower PTW marginal posteriors. This is expected since the truncated normal priors are more informative than the bounded flat priors. Note that the marginal posteriors of the non-PTW parameters appear to be unaffected by the change in the prior distribution. This suggests that prior knowledge about flow stress parameters does not aid in the inference of the HE and porosity parameters.

For a closer look at the HE parameter results, see Table VII and Fig. 6. In Table VII, we report for each HE parameter the maximum a posteriori (MAP) value as well as the posterior mean and standard deviation resulting from both calibrations. In Fig. 6, we show plots comparing the HE parameter posteriors for both priors. As noted above, all HE parameter posteriors are largely unchanged by the different flow stress priors. This implies that for the tantalum/LX-14 small plate experiment, HE model parameter uncertainty is largely insensitive to metal strength and, thus, cannot be reduced by well-characterized flow stress parameters.

To characterize the uncertainty in the model response based on the posterior parameter uncertainty, we utilized Monte Carlo error propagation. 8000 posterior samples from the truncated normal calibration were propagated through the emulator to obtain an ensemble of predicted velocity curves. Note that 95% Bayesian credible intervals provide a more quantitative alternative to Monte Carlo error propagation. In both cases, the resulting uncertainty envelope does not reflect the uncertainty due to measurement noise and, thus, is not expected to cover the data (as would be the case if 95% prediction intervals were employed). With this in mind, we validate our calibrated model by plotting the ensemble of predicted velocities with the experimental data in Fig. 7. Overall, the qualitative match to the data is good with a small excursion after the initial shock rise.

We again consider a small plate test but this time with LX-17 as the HE—that is, LX-17 is used to push a thin tantalum plate.² For this configuration, we neglect porosity given that results from previously conducted analysis suggested that this configuration was insensitive to such parameters. The PTW model is solely employed to characterize the metal’s flow strength. We also choose a different parameterization for the HE based on charge thickness and soot. The latter is related to the soot produced during the detonation process and is the molar volume used in the Murnaghan equation of state. This soot is an amorphous solid comprising carbon, nitrogen, and oxygen but not adhering to any specific molecular structure. The change in HE model construction with the addition of amorphous soot is reflective of recent soot recovery experiments and simulations.

A comparison of the resulting model response uncertainty from both error propagations is shown in Fig. 8. As with the previous experiment, the output uncertainty is greatly reduced—almost vanishingly small—when the samples are drawn from the narrower, previously calibrated posterior.

As in Sec. V 2, we obtain first order Sobol’ indices—this time for the LX-17 small plate model—as a function of time. The resulting Sobol’ indices are plotted in Fig. 9. The vast majority of the sensitivity lies with the HE soot parameter. Other parameters showing a discernible level of sensitivity include $y 0$⁠, $−log⁡(γ)$⁠, and the charge thickness. Note that the plot is nearly filled by the first order indices, indicating limited sensitivity due to parameter interaction. The exception to this trend is the early time regime, which suggests that parameter interaction may be influential before approximately 2.6   $μ$s.

Ten parameters, including eight PTW parameters and HE parameters soot and charge thickness, were used for calibration of our LX-17 small plate model. We compared calibration results obtained with two different prior distributions: one with bounded uniform priors for the PTW parameters and another with truncated normal priors for the PTW parameters. In both cases, we chose a bounded uniform prior for the soot parameter and a truncated normal prior for the charge thickness with the mean equal to the reported value and with truncation bounds given by $±10%$ of the nominal value; the soot parameter bounds were defined by the interval (4,8).

Utilizing these two priors allows us to explore the effect of better-characterized metal strength parameters on the calibration of the HE parameters. In Table VIII, the MAP, mean, and standard deviations for both calibrations are reported for the soot parameter and the charge thickness. The informed priors for the metal correspond to a modest reduction in uncertainty for the soot parameter, but the main point of interest is the effect of the shift in the mean on the predicted velocity. The shift in the two posterior distributions is evident in Fig. 10. While visual inspection does not suggest a radical difference in these posteriors, the HE modeling community considers a difference of $0.5%$ in the predicted velocity to be significant. Thus, we compared hydrocode calculations for two sets of HE parameters, corresponding to the posterior means from each choice of prior. The maximum velocities for the difference in soot parameter means is 0.31% and for charge thickness means it is 0.98%. This indicates that variations in charge thickness are something to consider and include for historical experiments given the value reported will be an average of multiple measurements on an nonuniform disk of HE.

As with the LX-14 small plate analysis, we characterized the model response uncertainty by performing an Monte Carlo error propagation. This was done by densely sampling the posteriors determined from the truncated normal PTW priors and propagating the resulting parameter sets through the Gaussian process emulator. In Fig. 11, the experimental data are plotted with the shaded area corresponding to the spread of the model responses. Recall from Sec. V 3 that the propagated error envelope need not cover the data. Nevertheless, we can validate the model by qualitatively assessing the goodness of fit to the data. Overall, the model and data match well with some late-time deviations. This discrepancy could be attributed to the 1D assumptions used in the computational model being no longer valid as 2D effects start to become dominate at late times.

Lastly, we analyze the most common experimental configuration for studying the detonation velocity of HE materials. In the cylinder test configuration, a collection of HE pellets are inserted into a metal cylinder and then detonated to expand the metal casing. Simulation requires a 2D axisymmetric computational model to accurately represent the experiment. The specific HE/metal configuration modeled here is LX-17 pushing a tantalum cylinder with the velocity being measured using Fabry–Pérot laser interferometry. As with the other experimental setups, we characterize tantalum flow stress with the PTW model. For porosity, we employ the classical Gurson porosity model as described in Sec. III A 2.

As with the small plate experiments, we propagate Monte Carlo samples of the PTW parameters from a bounded flat distribution and a previously calibrated posterior. With all other parameters fixed, we obtain the response ensemble shown in Fig. 12. Similarly to the other experimental configurations, drawing Monte Carlo samples out of a tighter distribution results in greater certainty with respect to the model response.

For our final sensitivity analysis, we calculate first order Sobol’ indices for the LX-17 cylinder configuration. The Sobol’ indices as a function of time are plotted in Fig. 13. As with the previous experimental types, the sensitivity is dominated by an HE parameter. Here, $E shift$ is by far the most sensitive parameter. Other parameters with visually discernible sensitivity include PTW parameters $κ$⁠, $−log(γ)$⁠, and $s 0$ as well as porosity parameter $q 2$ and HE parameter $P freeze$⁠. Moreover, some sensitivity to higher-order interactions is indicated for the time window before approximately 39 $μ$s.

Calibration of the LX-17 cylinder model involves 12 parameters: eight PTW parameters, two porosity parameters, and two HE parameters. The HE model used in this calibration is the same as the one used for the LX-14 configuration, employing energy density shift and chemical-freeze-out pressure parameters. We again carried out two calibrations, which employ priors encoding differing levels of information about the PTW parameters—specifically, a non-informative bounded uniform prior and a more-informed truncated normal prior. In both cases, the non-flow stress parameters were assigned bounded uniform priors (see Tables III and V for bounds).

We compare the resulting posteriors of the two calibrations to ascertain how knowledge of the metal strength can influence the posteriors of other parameters. Table IX provides a summary of the calibration for the HE and porosity parameters. Figures 14 and 15 show an overlay of both posteriors for the HE and porosity models, respectively. Of all of the configurations, the non-flow stress parameters of the LX-17 cylinder model showed the greatest sensitivity to the different PTW parameter priors. Modest differences are present in the HE posterior comparison with the posteriors from the more informative priors showing a slight reduction in uncertainty, but a notable effect is displayed with the porosity posteriors. Figure 15 shows a significant narrowing of the porosity model posteriors when a more informative prior is used for the PTW parameters. This suggests that incorporating flow-stress-related knowledge results in better-constrained porosity parameters.

Note the narrowing of the $q 2$ parameter, which controls the porosity evolution. Given that ductile failure of metals is related to the complexity of the loading condition and that the cylinder test will evolve non-uniform deformation as well as transients in strain rate, pressure, and temperature due to the complexity of HE loading, it is valuable to constrain porosity-related parameters for applications that extrapolate beyond the subset of focused experiments in this study.

In addition to the truncated normal priors for the PTW model, we include informed priors for the Gurson porosity model parameters with the mean being the average of uniform bounds and the standard deviation set to $±10%$ of the mean. This more informative prior will help us determine if a well-calibrated porosity model could further reduce the uncertainty of the HE parameters. Given the strong dependence on the global response of $E shift$ and $P freeze$ parameters, it is unsurprising that for a model calibration a tighter porosity prior results in a negligible difference in the HE posteriors.

As with the small plate configurations, Monte Carlo error propagation was used to study the effect of the parameter uncertainty on the model response uncertainty. We sampled parameter sets from the informed porosity posterior and propagated them through the Gaussian process emulator. Figure 16 shows the spread of the resulting response ensemble plotted with the experimental data. Visual inspection suggests that the model provides a good fit to the experimental data (Fig. 16).

Lastly, the posterior means and standard deviations. for the PTW parameters for all three configurations can be found in  Appendix, in Figs. 17–19, and compared to those used in the truncated normal priors in Tables X–XII. Notably, while there is some shifts and reduction in standard deviations, the differences are relatively small.

This article has presented a Bayesian analysis of high explosive-driven experiments, including plate and cylinder experiments of both LX-14 and LX-17. The goal of this work was to simultaneously assess the inelastic behavior of metals and the behavior of materials simultaneously from metal/HE experiments. We estimated the optimal parameters of Cheetah-based reactive-flow models of the energetic materials, the plasticity behavior, and the porosity behavior simultaneously. We used a relatively simple model for tantalum, the specific metal under consideration. Utilizing constraints from focused metal experiments, we were able to assess the effects of uncertainty in the metal material model on the metal/HE experiments. Note that the specific quantitative results described here for tantalum may not transfer to the use of other metals, such as copper, in coupled metal/HE experiments, but the methodology applied here is transferable to other HE/metal system.

Intriguingly, for our tantalum/HE models, significantly reducing uncertainty in the metal prior has minimal effect on the HE model parameters pertaining to the energy content of the material; however, there is a significant effect on other HE model parameters such as molar density of soot formed through detonation, $V m$⁠. Some energetic material parameters—such as $E s h i f t$ and $P f r e e z e$⁠, which are associated with the basic energy content of the HE—are strongly constrained by the experimental data and are relatively insensitive to the priors of the plasticity model. In contrast, inference of the soot density and charge thickness is sensitive to the choice of prior for the plasticity models. Analysis of the velocimetry is suggestive that the charge thickness appears to have small effect on the velocity prediction up to 0.9%. This implies that we may safely utilize historical HE data (with possibly high uncertainty in the reported charge thickness), but given the specificity of the models, data sets, etc. considered in our study, this should not be taken as absolute.

Based on our results, we offer recommendations for future statistical analysis and experimentation. When performing statistical inference, it would be beneficial to estimate charge dimensions when calibrating to historical data. Additionally, HE model parameters may benefit from relatively basic metal strength characterization when performing coupled metal/HE experiments for model calibration. Furthermore, we have observed some interesting interactions between the plasticity and the porosity model; prior knowledge of metal flow stress model parameters allows for improved calibration of metal porosity model parameters in cylinder tests. Hence, performing experiments with well-characterized metal is suggested. Next, the small plate calibration potentially indicates that insensitive HEs may exhibit more sensitivity to metal strength than conventional HEs. This indicates that when performing experiments with insensitive HEs it is especially important to use a well-characterized metal.

Lastly, HE lot to lot variability is not considered in this study. We close this study by offering some suggestions for future work in this area. First, we suspect that the decomposition of HE parameters into two classes—those having to do with the HE energy content and those which have to do with specific details of the products—may suggest future experiments and methods of analysis by which to learn about HEs; an obvious extension would be to explore the uncertainty related to HE density variations that are known to exist. Additionally, using a statistically motivated design of experiments strategy, additional experiments could be designed to maximize the information gained about parameters that are more challenging to define unambiguously when including the uncertainty arising from the metal. Second, we suggest that further work is needed to specifically examine the interplay between plasticity and porosity models; the correlations arising in our results among these parameters are consistent with the physics of the strength models. It will be important to utilize experiments in future analyses that are insensitive to the plasticity–porosity interaction, helping to disambiguate the causal behaviors. Finally, we endorse that analysis of model form error in the HE model as a valuable direction of inquiry and would likely have a profound impact on understanding uncertainty.

This work is performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DEAC52-07NA27344 (LLNL-JRNL-864794).

The authors have no conflicts to disclose.

Matthew Nelms: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). William Schill: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – review & editing (equal). I.-F. William Kuo: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal). Nathan Barton: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). Kathleen Schmidt: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

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

Figures 17–19 show the complete parameter posterior comparison for the LX-14 small plate configuration, complete parameter posterior comparison for the LX-17 small plate configuration, and complete parameter posterior comparison for the LX-17 cylinder test configuration, respectively. Tables X–XII show the PTW model parameters with a comparison of calibrations for smlplt-lx-14 and those used in the truncated normal prior distributions, PTW model parameters with a comparison of calibrations for smplt-lx-17 and those used in the truncated normal prior distributions, and PTW model parameters with a comparison of calibrations for cyl-lx-17 and those used in the truncated normal prior distributions, respectively.