Compaction of granular HMX : P-α porosity model in CTH hydrocode

Compaction waves traveling through porous cyclotetramethylene-tetranitramine (HMX) are computationally modeled using the Eulerian hydrocode CTH and validated with gas gun experimental data. The method employed use of a newly generated set of P-α parameters for granular HMX in a Mie-Gruneisen equation of state. The P-α model adds a separate parameter to differentiate between the volume changes of a solid material due to compression from the volume change due to compaction, void collapse in a granular material. Computational results are compared via five validation schema for two different initial-porosity experiments. These schema include stress measurements, velocity rise times and arrival times, elastic sound speeds though the material and final compaction densities for a series of two different percent Theoretical Maximum Density (TMD) HMX sets of experimental data. There is a good agreement between the simulations and the experimental gas gun data with the largest source of error being an 11% overestim...


I. INTRODUCTION
Energetic formulations are commonly comprised of granular explosives which have complicated heterogeneous microstructures.Shock waves traveling through this material compact the inherent voids and generate local temperature increases which are known as "hot spots".These hot spots are of interest as they increase the sensitivity of the energetic composite and may cause detonationintentional or otherwise.Porosity in granular explosives is one of several factors that govern reactivity thresholds and therefore influences the shock and impact sensitivity of the explosive.Accurate modeling of porosity is consequently important to predict the explosive reactivity under shock and impact conditions.This work summarizes an approach to model the change in density, and thus porosity, of the explosive cyclo-tetramethylene-tetranitramine (HMX) due to an imparted shock wave.Simulations are compared to existing Los Alamos National Laboratory (LANL) experimental gas gun data which captured velocity and stress as measured by magnetic particle velocity and polyvinylidene diflouride (PVDF) gauges.The compaction wave experimental setup and results are described in Sheffield, et al. 1 Previous computational simulations of similar experiments detailed in Ref. effects from elastic-plastic interactions as a function of yield strength of HMX grains in mesoscale calculations.Menikoff's strength based model is analogous to the P-α model in that it relates the average material stress to the stress of the solid matrix material. 3The goal herein is to simulate the Sheffield LANL gas gun experiments in the hydrocode CTH with macroscale porous (P-α model) granular HMX parameters to validate the modeling methodology employed in the present work.An outcome of this work is the generation of a working P-α parameter set for granular HMX in CTH.
Computational model setup is illustrated in FIG. 1 for the simulated LANL gas gun experimental configuration with 65% theoretical maximum density (TMD) HMX impacted by a polychlorotrifluoroethylene (Kel-F) flyer.The granular HMX bed is positioned between a Kel-F front plate and polymethyl-methacrylate (PMMA) back plate.The aforementioned magnetic particle velocity and PVDF instruments are positioned at the front plate to HMX interface, as well as the HMX to back plate interface to capture velocity and stress data.Numerous experimental trials were repeated to obtain complete data sets as individual shots were instrumented with either magnetic particle velocity or PDVF gauges, not both simultaneously.
Results are presented for one particular test case series using fine grain HMX with a mean particle size of 10-15 µm.The density is 1.24 g/cm 3 or 65% TMD, corresponding to a porosity of 35%.The Kel-F flyer velocity is 279 m/s.However, there exists a comparable set of LANL gas gun test data with granular HMX density of 1.4 g/cm 3 , 74% TMD or 26% porosity, and a flyer velocity of 270 m/s.These 1.4 g/cm 3 simulations were likewise repeated but are not discussed in the body of the paper.Results for the 1.4 g/cm 3 simulations are summarized in the concluding sections in conjunction with the 1.24 g/cm 3 HMX simulations.

A. Computational Model and Parameters
CTH version 10.2 is utilized for the present analyses.CTH is a three-dimensional multi-material Eulerian hydrocode capable of modeling high strain rates characterized by high velocity impact, shock wave transmission through dissimilar materials and shock wave coalescence.All solid materials in this analysis are modeled with a pressure dependent Mie-Gruneisen Equation of State (EOS) formulated from Hugoniot curve parameters.
To capture these effects for the LANL gas gun experiment in CTH it was necessary to formulate an inert Mie-Gruneisen user defined EOS, described in greater detail in Ref. 4, for granular HMX from multiple references with additional P-α porosity model parameters.The P-α model was originally formulated by Hermann 5 and was intended to accurately model porous compaction in a thermodynamically consistent manner at high stresses and approximate compaction at lower stresses.In Herrmann's porosity model the distention parameter α = v/v s , is defined as the ratio of specific volume of a porous material, v, to the specific volume of a solid matrix material, v s .
A fundamental assumption of Herrmann's model is that the specific internal energy is equivalent for the porous material and solid matrix material when at the same temperature and pressure conditions.The primary function of α is to distinguish the volume change of the solid matrix material due to material compression from the volume change of the porous material due to void collapse (compaction). 5Herrmann's P-α model recasts an EOS where pressure is a function of specific volume of a solid matrix material and internal energy to include porous effects by making pressure a function of porous material specific volume divided by α as well as internal energy, shown in Eqn.(1).
Herrmann separates the effects of compaction into two regimes: elastic and plastic.Assuming that shear strength is negligible and consequently sound speed can be appropriately represented by only the bulk sound speed, the porosity term α p can be expressed in the form below, Eqn.(2).In this equation from Herrmann's original 1969 publication, the subscript p denotes the plastic deformation regime, p s is the pressure at which compaction is complete where α = 1, p e is the elastic compaction region pressure threshold, and α 0 is the initial porous material distention parameter.
In 1971, Carroll and Holt of Lawrence Livermore National Laboratory published a notable augmentation of the original formulation of the P-α model.Carroll and Holt modified the pressure dependent EOS relation in Eqn.(2) to include a factor of 1/α, as displayed in Eqn.(3), to correct for the observation that the pressure in the porous material is approximately equivalent to 1/α times the solid matrix material average pressure value. 6 Mie-Gruneisen EOS with P-α porosity parameters implemented in CTH to model porous granular HMX was based on iteratively determined parameters as well as a cross-referenced compilation of parameters published in Refs. 2 and 7-10 formatted in accordance with the CTH input instructions. 11P-α model parameters for granular HMX are listed in Table I.Mie-Gruneisen parameters for Kel-F were obtained from Menikoff 2 and Mie-Gruneisen parameters for PMMA were obtained from the CTH EOS library.CTH simulations included air in the regions of the computational domain unoccupied by solid materials.Strength of all materials was represented with Elastic-Perfectly-Plastic strength models.
CTH parameter sensitivities were evaluated throughout this investigation.By iterative computational analysis it was determined that the minimum pressure at which compaction begins parameter, p e , influences the wave front shape.A p e parameter value of 0.1 MPa results in a planar wave front whereas higher values of this parameter p e result in an increasingly pronounced non-planar curved wave front.CTH simulation results were quite sensitive to the value assigned to the pressure at which  12 Both Lagrangian and Eulerian discrete user defined tracer points were defined in the CTH porosity model validation simulations to record local pressure (stress), velocity, and density data as a function of time.Similar to the configuration modeled in Menikoff, 2 the present simulations utilized rows of approximately 50 tracers one baseline coarse mesh cell distance into the interface boundaries spaced 0.01 cm apart linearly.This tracer configuration facilitates comparison of hydrocode results to experimental data in Sheffield, et al. 1 Specific variable results were averaged at each time, recorded in 0.1 microsecond increments.

B. Validation Schema
In order to demonstrate that CTH is modeling the pertinent physics correctly, validation of the P-α porosity model in CTH followed a five-pronged multi-parameter validation approach.This validation schema included comparison of experimental data or analytical calculations (where noted) to CTH in terms of: Each of these schema are addressed in the section below, the compendium of which supports a robust modeling technique and further application of the iteratively determined P-α parameters.

C. CTH Simulation Summary
Prior to delving into results and validation of the aforementioned Approaches 1-5, it is useful to first review time sequences of pressure contour plots to examine the simulation behavior.FIG. 2 contains CTH simulation results for a simulated LANL gas gun experiment with granular HMX density of 1.24 g/cm 3 , 65% TMD, and piston velocity of 279 m/s with a fine-mesh spacing of 10 µm.
Progress of the compaction wave through granular HMX is evident in the images from three to six microseconds given in one microsecond intervals.In this configuration the planar wave reaches the beginning of the HMX sample at 2.8 µs.Piston impact imparts a plastic compaction planar wave of approximately 230 MPa into the granular HMX.Plastic compaction wave pressure, as well as wave front velocity, decreases as a function of time as the compaction wave travels through the granular energetic material sample.An elastic precursor compaction wave traveling at approximately the sound speed of granular HMX leads the planar compaction wave.Edge effects minimally influence the wave shape in HMX at later times, shown below at six microseconds by rounding of the planar wave edges due to rarefaction waves impinging at the granular HMX to air interface.The thick black line near y = 2 cm denotes tracers positioned within HMX in close proximity to the back interface.Edge effects are more pronounced in the region of the Kel-F front flyer and plate.Finite plastic compaction wave thickness shown in the pressure contour plots supports Sheffield's assertion that "porous materials do not propagate sharp shock waves.Instead, the waves are spatially and temporally diffuse or spread out 1 ".
LANL experimental gas gun data for both "coarse" 120 µm mean particle size and "fine" 10-15 µm mean particle size tests were compared to CTH simulation results on a broad range of mesh spacing, from 10 to 200 micrometers, in order to determine which grain size distribution CTH more closely replicates.CTH simulations with appropriately tuned P-α inputs compare very well to back plate experimental velocity versus time traces for fine mean particle size HMX.CTH significantly underestimates back plate velocity for experimental data collected with coarse grain HMX particle distributions.The P-α porosity model assumes a uniform distribution of porosity throughout a matrix material and does not consider statistical particle distribution and packing.Uniform porous distribution is therefore more closely replicated in samples with consistent and small interstitial spaces, referred to as void regions in the study of granular energetic materials, which is more characteristic of fine grain samples.Large variations in void size throughout a coarse grain sample do not match the simplifying assumptions inherent in the P-α porosity model.Subsequent validation cases are only compared to fine grain, 10-15 µm mean particle size, experimental granular HMX data.CTH simulation results included in figures and corresponding discussions are for simulations performed on the finest mesh spacing of 10 µm.FIG. 4 in the subsequent Approaches 2 and 3 results discussion contains velocity profiles obtained on a 10 µm computational mesh to demonstrate the importance of fine mesh spacing in low velocity impact porous compaction simulations.

II. GAS GUN COMPUTATIONAL COMPARISON
The following section contains comparisons of stress and velocity experimental results for LANL gas gun cases with fine grain HMX at a density of 1.24 g/cm 3 to CTH simulations.In order to lend confidence to the presented p-α parameters, the conclusion contains tabulated simulation data for comparison to simulations with HMX at a density of 1.4 g/cm 3 .

A. Validation Approach 1: Stress
Validation Approach 1 compares the Sheffield experimental PDVF gauge stress data at the Kel-F front plate to approximately 50 CTH Lagrangian tracer values spaced 0.01 cm apart linearly and positioned 200 µm from the interface boundary to avoid shared volume fraction calculation and impedance issues.FIG. 3 plots experimental PDVF Kel-F front plate input wave pressure versus time along with corresponding CTH results post-processed in Matlab.The CTH results have been shifted in time to facilitate comparison with experimental results.Time zero occurs when the input wave reaches the Lagrangian tracers of interest.This time shift is consistently employed in subsequent plots.
Sheffield 1 data plotted with black filled circles in FIG. 3 correspond to experimental shot number 2477 with a piston velocity of 285 m/s, coarse grain HMX, and a TPX back plate.This is the only input experimental stress data provided in Sheffield.CTH simulations replicated Sheffield experiments with fine grain distribution 1.24 g/cm 3 HMX samples with a piston velocity of 279 m/s and a PMMA back CTH results and experimental data in FIG. 3 compare well near time zero when the input wave reaches the granular HMX sample after traveling though the Kel-F front plate.The slight depression at the beginning of the experimental data set is not evident in CTH results.This depression may be due to impedance mismatch at the gauge interface position between the Kel-F front plate and HMX sample.CTH results deviate from the experimental trend at 2.2 µs, likely due to impingement of rarefaction waves at tracers positioned closet to the edge.At the crest of the initial empirical depression in digitized results, occurring at 0.39 µs, CTH results at the nearest data capture point of 0.4 µs are within 11% of the experimental data.Part of the difference in CTH and empirical input wave data can be attributed to slight differences in simulation versus experimental conditions.The CTH simulations were intended to replicate fine grain HMX back plate velocity data.Available front plate data piston velocities are 2% higher (285 m/s for coarse HMX Shot 2477 versus 279 m/s for fine grain HMX) leading to slightly higher input wave pressure.

B. Validation Approaches 2 and 3: Velocity and Arrival Time
Validation Approach 2 concerns the comparison of velocity versus time traces at the front Kel-F plate and back PMMA plate for the same CTH simulation scenario of fine grain HMX at a density of 1.24 g/cm 3 and a piston velocity of 279 m/s.As with front plate input wave PDVF stress gauge data, magnetic particle velocity data at the Kel-F to granular HMX interface were only available for Sheffield Shot 912 configured with 1.24 g/cm 3 coarse grain HMX with a piston velocity of 285 m/s and a TPX back plate.Though it is appropriate to compare upstream input wave front plate profiles despite the differences in HMX coarse and fine grain samples, slightly higher input wave velocity profiles are expected due to the 2% higher piston velocity in available input wave front plate velocity data.
FIG. 4 comprises subplots of front and back plate velocity versus time profiles.Subplot (a) contains front tracer data compared to experimental results 1 from Shot 912 plotted with black filled circles.CTH results on the 10 µm mesh compare quite well to digitized experimental data at early times in the simulation.A break away in data trends is evident again around 2.2 µs.Prior to this trend line deviation, which is likely due to rarefaction wave impingement, CTH results on the two finer meshes are within 6% of experimental results from 0 to 2.2 µs.A portion of the CTH under-predicted velocities can be attributed to the 2% lower piston velocity in CTH.Velocity spikes at the wave front in the computational results were also observed in yield strength based porosity model hydrocode calculated in Ref. 3. Menikoff and Kober 3 attribute this velocity spike to the "blow-off velocity at the free surface and is a consequence of the pores between grains."They further note that the velocity Once the input wave begins transmission through and compaction of granular HMX two separate waves propagate through the CTH simulated HMX sample.This phenomenon is consistent with the distinction of elastic and plastic compaction regimes in the P-α model. 5The concept of a discontinuity propagating through a solid material at a speed less than the sound speed is explained in detail in Ref. 13. Powers et al. 13 assert that either subsonic or supersonic compaction waves can exist, depending on piston impact velocity, and a critical piston velocity threshold separates the two regimes.Subsonic compaction waves are "characterized by a smooth rise in pressure from ambient to a higher pressure equal to the static pore collapse stress level" while supersonic compaction wave behavior is that of "a discontinuous shock leads a relaxation zone where the pressure adjusts to its equilibrium static pore collapse value."Subsonic compaction has been experimentally documented while supersonic compaction waves have yet to be experimentally observed. 13In the piston velocity range of interest for this study, approximately 280 m/s, subsonic compaction waves occur.The average compaction wave velocity is substantially lower than the bulk sound speed of HMX, 2,740 m/s, according to Refs.7 and 10.Note that the cited sound speed value is for solid matrix HMX crystals at a theoretical maximum density of 1.903 g/cm 3 .Porosity in a granular energetic material is postulated to lower the bulk sound speed because the wave is transmitted through contacting grains and therefore has a longer path to travel in granular materials. 3alidation Approaches 2 and 3 distinguish between the rise time, more generally the velocity versus time curve profile, and the arrival time of the transmitted plastic wave through the granular HMX as two separate validation criteria.This distinction has been made to differentiate between the finite rise time at a single geometric location characteristic of subsonic compaction waves and the arrival time of the wave.Shock velocity, U s , changes as a compaction wave traverses a granular material sample and thus obtaining an accurate arrival time represents accurate modeling of the average plastic compaction wave transmission throughout the thickness of granular HMX.
Subplot (b) in FIG. 4 is a comparison of empirical magnetic particle velocity gauge data, plotted with black filled circles, for a 1.24 g/cm 3 fine grain HMX experimental configuration with piston impact velocity of 279 m/s and a PMMA back plate. 1 CTH simulation input conditions replicate the known experimental conditions.Lagrangian tracers are utilized to capture data adjacent to the HMX to PMMA back plate interface.A smooth transition and 0.1 µs rise time is observed in the experimental data. 1 CTH results on the 10 µm mesh contain a slope change 0.1 µs into the wave arrival.Measuring the rise time from only the single sharp slope region the calculated rise time is 0.2 µs.
These results are indicative of the need for very fine mesh spacing to accurately resolve inert porous granular energetic material behavior.Again, one source of error may be due to the date write interval of 0.1 µs.Back plate velocity values at the peak of the arrival wave are within +5% on the 10 µm mesh.For comparison, 200 µm mesh results completely diffuse the arrival wave profile, with a rise time of approximately 2 µs, and underestimate the back plate interface velocity by nearly 50%.
In validation Approach 3 the transmitted wave arrival times are of primary focus.As previously stated, obtaining comparable arrival times implies that compaction behavior throughout the granular sample has been accurately approximated due to the degradation of shock velocity as a function penetration distance into the granular compact.Transmitted wave arrival occurs at 5.12 microseconds in the digitized experimental data.Past 30 µm mesh results have an arrival time of 4.8 µs, within 6% of experiment, as compared to the start of an upward velocity trend in the right subplot CTH results of FIG. 4 which indicate that the 10 µm mesh calculates an arrival time of 4.9 µs, within 4% of empirical data.

C. Validation Approach 4: Elastic Compaction wave
Validation Approach 4 compares published HMX longitudinal sound speed to the elastic precursor wave speed calculated in CTH.The elastic wave speed is calculated at the same tracer locations as the back plate velocity profiles to avoid impedance matching issues with the adjacent plate.For a granular HMX sample thickness of 0.39 cm, the elastic precursor arrival time in CTH is 1.4 µs, corresponding to a shock velocity of 2,785 m/s.For the previously cited longitudinal sound speed value of 2,740 m/s for crystalline HMX, a calculation error of 2% is obtained.These results further validate the excellent agreement between CTH and the P-α porosity model.Lower sound speeds theorized to exist for granular energetic materials are not observed due to simplifying assumptions inherent in the P-α porosity model.

D. Validation Approach Five: Hugoniot Density
The final validation point, Approach 5, is a comparison of analytical final compaction density calculations as well as additional experimental results to CTH results at the center thickness of the HMX sample.A plot of density versus time of an Eulerian fixed tracer particle positioned at the granular HMX center is provided in FIG. 5. Final compaction density can be derived from Hugoniot jump conditions for pressure and density as a function of shock and particle velocity, as detailed in Ref. 14.The resultant non-dimensional equation is given below in Eqn.(4).In this equation, ϕ is the percent density ratio equal to 100 × (ρ/T M D).Subscript 0 denotes initial granular material state and FIG. 5. Geometric Center HMX Density vs Time of LANL Gas Gun Simulation at ρ = 1.24 g/cm 3 .subscript h denotes the Hugoniot jump condition state.TMD is assumed to be the previously cited HMX crystal density of 1.903 g/cm 3 .Shock velocity, U s , changes as a function of time and position in the granular material during compaction.The particle velocity, u p , levels off to a steady-state value at a discrete location following attenuation of the input shock from transmission through the front plate.
Shock velocity is calculated via the arrival time of the plastic compaction wave at the granular HMX center thickness, 0.39 cm / 2 = 0.195 cm.For an arrival time of 2.5 µs, the corresponding shock velocity is 780 m/s.Particle velocity is initially equal to the piston velocity of 279 m/s and decreases as the compaction wave traverses the granular material.Particle velocity at the HMX sample center is measured from a Lagrangian tracer and is equal to 220 m/s.This particle velocity compares well with the steady state experimental front plate velocity oscillations in the range of 230-235 m/s in the left subplot of FIG. 4. Further particle velocity degradation is expected by the time the plastic compaction wave reaches the granular HMX center.With the known values of initial density as well as shock and particle velocity, the corresponding final compaction density is calculated via Eqn.(4) as 1.73 g/cm 3 , 91% TMD, and plotted on FIG. 5 with black filled circles to facilitate comparison to CTH results.Excellent agreement, within 1%, is achieved between the Hugoniot theoretical compaction density and CTH results.
A brief literature review was also conducted to determine if supporting experimental results intended to measure final compaction density have been published.According to Ref. 14 there are four different experimental techniques to test dynamic compaction of granular energetic material.These four experimental methods include: gas driven, piston driven, ramp loaded, and shock driven compaction experiments.Results with the ramp loaded technique are reported in Ref. 15 for 64.6% TMD, or 1.23 g/cm 3 .Ramp loaded experimental results are reported in the final compaction range of 86.5%-97.3%TMD for 64% HMX at maximum stresses in the range of 18.7 to 227 MPa.Variations in results as compared to the present gas driven compaction mechanism are expected given the differences in experimental approach.

III. SUMMARY OF RESULTS AND CONCLUDING REMARKS
To extend the results further, a simulated LANL gas gun experiment with granular HMX density of 1.4 g/cm 3 , 74% TMD, a piston velocity of 270 m/s, and PMMA back plate simulations were likewise completed with various mesh spacing.The simulations follow the same trend as the 1.24 g/cm 3 HMX samples.CTH simulation results for HMX densities of 1.24 and 1.4 g/cm 3 are summarized in Table II with simulation versus experimental conditions listed in Table III  a The 100% relative error is indicative of the simulation data-write interval of 0.1 µs, thus the CTH rise times are not resolved on a sufficiently fine interval.This is an area for future refinement.

Approach 1 :
PDVF gauge stress data at the Kel-F front plate Approach 2: Magnetic particle velocity data velocity vs time traces at the front plate and back plate, specifically focused on rise time and wave profile Approach 3: Magnetic particle velocity data arrival time of the plastic compaction wave, specifically instantaneous velocity at the back plate.Approach 4: Analytical comparison of elastic compaction wave speed at the back plate to longitudinal sound speed Approach 5: Analytical Hugoniot comparison to final compacted density at the average shock velocity

FIG. 4 .
FIG. 4. LANL Velocity Magnitude in Gas Gun Simulations with Fine HMX ρ = 1.24 g/cm 3 and V = 279 m/s.(a) Velocity Profile at the Kel-F Front Plate.(b) Velocity Profile at the PMMA Back Plate.

TABLE I .
Granular HMX P-α Parameters.Pressure at which plastic compaction is complete, p s 300 MPa Pressure at which elastic compaction is complete, p e 0.10 MPa compaction is complete, at 100% TMD.Mesh independent CTH results for the LANL gas gun experimental configuration showed that decreasing the pressure at which compaction is complete increased the final compacted TMD percentage, increased transmission time of the compaction wave through the granular energetic material and lowered back plate average velocity, and vice versa for increasing p s .Iterative parameter determination initially began with p s equal to 230 MPa based on simulation results in Menikoff, 2 specifically from the input wave pressure in similar Sheffield et al. experiments where Menikoff predicts near complete compaction.Setting p s to 300 MPa produces CTH results with excellent agreement with experimental data and is consistent with calculations in Baer & Nunziato.
. Table entries containing N/A indicate that no published experimental data was available for comparison.

TABLE II .
Comparison and Relative Error of Validation Approaches.