The response of high explosives (HEs), due to mechanical and/or thermal insults, is of great importance for both safety and performance. A major component of how an HE responds to these stimuli stems from its reactant equation of state (EOS). Here, the tri-amino-tri-nitro-benzene based explosive PBX 9502 is investigated by examining recent experiments. Furthermore, a complete thermal EOS is calibrated based on the functional form devised by Wescott, Stewart, and Davis [J. Appl. Phys. **98**, 053514 (2005)]. It is found, by comparing to earlier calibrations, that a variety of thermodynamic data are needed to sufficiently constrain the EOS response over a wide range of thermodynamic state space. Included in the calibration presented here is the specific heat as a function of temperature, isobaric thermal expansion, and shock Hugoniot response. As validation of the resulting model, isothermal compression and isentropic compression are compared with recent experiments.

## I. INTRODUCTION AND MOTIVATION

The thermodynamic response of high explosives (HEs) is of paramount importance for determining the response to mechanical and thermal stimuli. During impact of projectiles/fragments, the pressures and temperatures achieved in the HE are determined solely by its reactant equation of state (EOS). Likewise, the EOS is responsible for determining the expansion/contraction during heating/cooling of the HE. This expansion/contraction leads to changes in the initial energy density of the HE, and thus the Chapman-Jouguet detonation velocity upon detonation. Other observable quantities such as the shock-polar matching between the HE and surrounding inert material rely on knowing the HE reactant EOS,^{2–5} as well as the process of shock transmission across the interface between the HE and a second material. The thermodynamic response is critical for future, higher fidelity, reactive flow modeling that includes effects due to changes of the initial HE temperature^{6,7} and in pseudo-reaction zone modeling of HE.^{8}

In a series of papers,^{1,9,10} a complete thermodynamic model of reactant EOS was developed. In those papers, a procedure for calibration and an initial set of the relevant thermodynamic parameters for PBX 9502 [95% tri-amino-tri-nitro-benzene (TATB), 5% polymeric binder Kel-F 800] were determined utilizing the data available at the time.^{1} The focus of the present study is to revisit the role of each of the model parameters to obtain a wider range of thermodynamic validity by incorporating experimental data not included in the above described calibration.

The outline of the remainder of the paper is as follows: Sec. II details the analytical model reactant equation of state along with a few newly derived relationships between model parameters, thermodynamic variables, and measurable quantities observable from experiments. This is then followed by Sec. III, which outlines a calibration procedure for systematically determining EOS parameters. Some derived properties of the model are presented in Sec. IV. Section V presents validation tests of the calibrated model via isentropic and isothermal compression experiments. Finally, Sec. VI gives conclusions and possible improvements for future work.

## II. DAVIS EQUATION OF STATE FOR HE REACTANTS

A thermally complete Mie-Grüneisen equation of state based on a reference isentrope has been devised, and has been utilized in reactive flow models of HE.^{1} Following Wescott *et al.*, for a reactant EOS, the relations between specific internal energy, *e _{r}*, pressure,

*p*, and density,

_{r}*ρ*, for the reactants (thus the subscript

*r*) are given by the standard Mie-Grüneisen form

and

where Γ_{r} is the Grüneisen gamma, here taken to be only a function of *ρ*. The reference isentrope, denoted by a superscript *s*, of the reactants is given by the following piecewise functional form:

where *y* = 1 – *ρ*_{0}/*ρ* and $p\u0302=\rho 0A2/4B$. Here, *A*, *B*, and *C* are material dependent parameters and *ρ*_{0} is the reference ambient density. Note that the above two functional forms are continuous up to *y*^{3} at *ρ* = *ρ*_{0}. Furthermore, the above $prs(\rho )$ form is monotonically increasing with respect to *ρ*, assuming that the parameters *A*, *B*, and *C* are positive. Thus, along the principal isentrope $dprsd\rho >0$, and a real sound speed will result, even when the material is distended. Note that even though the isentropic bulk modulus is always positive, the isothermal bulk modulus becomes negative at sufficiently low enough density, which is a common issue with solid EOS models.^{11} If no switch is employed in Eq. (3) (and *C* > 0), and the second term is used for all *ρ* values, then there would always come a point for low enough density where $dprsd\rho <0$, implying an imaginary sound speed. Given this fact, the switch in the above equation is required to maintain reasonable mathematical behavior when *ρ* < *ρ*_{0}. An example plot of both functions is shown in Fig. 1.

Two other properties of the above fitting form are worth noting. First, by taking the derivative of $prs(\rho )$ with respect to *ρ*, it is easy to verify that the parameter *A* is the ambient sound speed of the material at *p* = 0 and *ρ* = *ρ*_{0}. The isentropic bulk modulus is defined as

along the principal isentrope. The bulk sound speed is related to the ambient isentropic bulk modulus, *K*_{0}, via *K*_{0} = *ρ*_{0}*A*^{2}. Furthermore, by taking the derivative of the isentropic bulk modulus with respect to pressure along the principal isentrope, the following is obtained:

Evaluating this function at *ρ* = *ρ*_{0}, it is found that

The slope, *s*, of the Hugoniot in shock velocity, *U _{s}*, and particle velocity,

*U*, space at

_{p}*U*= 0 is given by

_{p}^{11}

Thus, when *B* ≫ 1, the initial slope of the Hugoniot in *U _{s}* –

*U*space is

_{p}*s*≈

*B*. The

*C*term, in Eq. (3), was added to modify the Hugoniot at high pressures, to ensure that the reactant Hugoniot does not cross the product Hugoniot.

^{1}

From the thermodynamic identity *de *=* Tds* – *pdv*, for an isentropic process, simply *de* = –*pdv* is attained. Upon substitution of specific volume, *v*, in favor of density, the following is obtained: $de=p\rho 2d\rho $. The above reference pressure along the isentrope can then be integrated to obtain the reference energy along the isentrope

where *E*_{0} is the integration constant. This integration constant is often taken to be the stored chemical potential energy of the explosive, but simply represents an offset in examining just the reactant EOS. For inert solid EOS, this integration constant is often taken to be zero, yielding zero energy at the reference state.

The Grüneisen parameter is taken to be

where *Z* is a constant used to describe the changes to Γ_{r} with respect to density. Note that the *ρ* < *ρ*_{0} branch is not explicitly mentioned in the original Ref. 1, but is often implemented in practice. Note also that *y* = 0 at *ρ* = *ρ*_{0}, *y* → 1 as *ρ* → *∞* and *y* → –*∞* as *ρ* → 0. Considering the fact that the sound speed^{11} contains derivatives of Γ_{r}(*ρ*), there will be a discontinuity in the sound speed at *ρ* = *ρ*_{0} for nonzero pressures and *Z *≠ 0.

The reactant temperature, *T _{r}*, after a series of manipulations,

^{1}can be obtained as a function of energy and density

where $Trs(\rho )$ is the temperature along the reference isentrope

Here, $Cvr0$ is the reactant specific heat at constant volume at the reference temperature, *T*_{0}. The parameter *α _{st}* determines how the specific heat changes with respect to temperature. In particular, Eq. (10) can be inverted to determine

*e*(

*ρ*,

*T*), then holding

*ρ*fixed and taking the derivative with respect to

*T*yields the specific heat at constant volume

Note that at the reference density, $Cvr(\rho 0,T)=Cvr0(TT0)\alpha st$. This functional form for the specific heat at constant volume is a power law, and is not fundamental, but rather a consequence of choices made in assuming the specific heat at constant volume was a linear function of entropy. Other choices, in general, can be made.^{12}

It is instructive to obtain pressure as an explicit function of density and temperature. Upon substituting energy in favor of pressure from Eq. (1) into Eq. (10), and solving for pressure, one obtains

From this formula, isochoric pressurization can be obtained while keeping the density fixed. This is the pressurization due to heating at a fixed volume.

Similarly, the pressure along a specified temperature may be obtained by varying the density. This is often accomplished in diamond anvil cells for HE constituents, and is discussed in Sec. V B.

## III. CALIBRATION PROCEDURE

For the reactants EOS, the following parameters are required to define the model: *A*, *B*, *C*, *ρ*_{0}, *Z*, $\Gamma r0$, *E*_{0}, $Cvr0$, *α _{st}*, and

*T*

_{0}. Since

*E*

_{0}is associated with the stored chemical energy of the explosive, it is not investigated here. The specified initial pressing density range of PBX 9502 (Ref. 13) is from 1.885 g/cm

^{3}to 1.895 g/cm

^{3}, so a nominal initial density will be taken as

*ρ*

_{0}= 1.890 g/cm

^{3}. Here, the reference temperature is taken to be

*T*

_{0}= 297 K.

### A. Specific heat

It is difficult to directly measure the specific heat of a condensed phase material at constant volume over a wide temperature range. This is especially true for high explosives^{14} that decompose at temperatures from 465 to 600 K, yet there is interest in knowing thermal properties such as the shock temperatures under detonation conditions. These temperatures may reach 2000 K under shock loading, even prior to decomposition into detonation products. Considering that most high explosive molecules are rather large, and thus have many vibrational modes, the specific heat is expected to change significantly over these large temperature ranges. So, the best that can be done is to model^{15} the specific vibrational modes and frequencies^{16} to build a thermal model of *C _{vr}*. Here, the parameters $Cvr0$ and

*α*in Eq. (12) are fit to the model presented by Menikoff.

_{st}^{15}The resulting parameters, determined by least squares fitting over 0 <

*T*< 1800 K, are presented in Table I. A comparison between Menikoff's determined

*C*, the newly calibrated model, and the original model

_{vr}^{1}is shown in Fig. 2. Note that the original calibration matches Menikoff's work at relatively low temperatures, but the current calibration is reasonable over a wider range of temperatures; all models agree near ambient conditions. It should be pointed out that the functional form in Eq. (12) does not lend itself to simultaneously having the proper asymptotic behavior of $dCvrdT=0$ in the limits of

*T*→ 0 and

*T*→

*∞*. The best fit of Eq. (12) typically will overshoot a Debye model at both high and low temperatures, while being too low at some intermediate temperature range. This is seen in Fig. 2.

Parameter . | Value . | Unit . |
---|---|---|

A | 1.80 | mm/μs |

B | 4.6 | |

C | 0.34 | |

ρ_{0} | 1.890 | g/cm^{3} |

Z | 0.0 | |

$\Gamma r0$ | 0.56 | |

$Cvr0$ | 0.001074 | kJ/g K |

α _{ST} | 0.4265 | |

T_{0} | 297 | K |

Parameter . | Value . | Unit . |
---|---|---|

A | 1.80 | mm/μs |

B | 4.6 | |

C | 0.34 | |

ρ_{0} | 1.890 | g/cm^{3} |

Z | 0.0 | |

$\Gamma r0$ | 0.56 | |

$Cvr0$ | 0.001074 | kJ/g K |

α _{ST} | 0.4265 | |

T_{0} | 297 | K |

### B. Hugoniot and isobaric thermal expansion

The Hugoniot of many HE, including PBX 9502,^{17,18} has nonlinear relations between the shock velocity, *U _{s}* and particle velocity,

*U*. There is also a fair amount of uncertainty in the Hugoniot as shock velocities increase to the point where chemical reactions begin to take place. Furthermore, plastic bonded explosives are heterogenous and somewhat porous, which further complicates continuum modeling. As shown earlier, the parameters

_{p}*A*and

*B*can be related to the sound speed and the initial slope of the Hugoniot curve in

*U*–

_{s}*U*space, respectively. From previous studies,

_{p}^{1,17,19}there is a range of sound speeds from 1.75 mm/

*μ*s to above 2 mm/

*μ*s. Here, an intermediate value of

*A*= 1.8 mm/

*μ*s was chosen. Likewise, the initial slope of the Hugoniot curve varies significantly from 3.15 to 5.2 in the references. Again, an intermediate value of

*B*= 4.6 was chosen. Note that even in isothermal compression experiments, performed in diamond anvil cells, high uncertainty in the isothermal bulk modulus (30% variability) and its derivative with respect to pressure (100% variability)

^{20}are observed.

The isobaric volumetric compression/expansion as a function of temperature is governed by a variety of thermodynamic parameters. Once the specific heat parameters, $Cvr0$ and *α _{st}*, and the principal isentrope parameters,

*A*and

*B*, have been set, only the Grüneisen parameter is left to determine how the density varies with temperature at zero pressure (the

*C*parameter plays a minor role, and only at the cold temperature extremes). Namely, the thermodynamic relation

can be used to estimate thermal expansion and contraction. Here, *β* is the volumetric thermal expansion coefficient and *C _{p}* is the specific heat at constant pressure, which is typically only a few percent higher than

*C*for condensed phase materials. To go beyond estimation, one can compute the isobaric thermal expansion analytically by noticing that the term $(e\u2212ers(\rho ))$ in Eq. (10) becomes simply $\u2212prs(\rho )/\rho \Gamma r(\rho )$ when

_{vr}*p*= 0, from examining Eq. (1). This leads to a rather simple piecewise expression (depending on whether in contraction or expansion) for the temperature as a function of density. This can be inverted to obtain density as a function of temperature, which can be compared to experimental observations. Equivalently, one can set the pressure on the left side of Eq. (2) to zero, and solve for

*T*(

*ρ*).

There have been several measurements of thermal expansion/contraction of PBX 9502 in the literature.^{21–26} TATB, the main component of PBX 9502, is quite anisotropic in its thermal expansion, which can lead to uncertainty in the bulk thermal expansion and contraction. And although the HE molding powder has random initial crystal orientation, once the material is pressed into PBX 9502, the crystals tend to reorient depending on the applied pressure direction during pressing. This leads to pressed parts which are then anisotropic. Another contributing factor to the uncertainty of the experimental measurements is the fact that PBX 9502 can experience “ratchet growth” upon cyclic heating and cooling,^{24} which leads to hysteresis. These complications are obviously beyond the scope of this fluid-like continuum EOS modeling, but should be kept in mind as possible sources of uncertainty and are worthy of future research.

To match the experimental thermal expansion, $\Gamma r0=0.56$, is determined. The variation of the Grüneisen parameter with density is not used in this calibration, and so *Z* = 0 was taken. Taking *Z* = 0 also avoids the discontinuity of sound speed at the reference density. The isobaric thermal expansion for the original calibration, current calibration, and experimental data are shown in Fig. 3. For temperatures below 350 K, the current and original calibrations agree within about 1% in density. But, as the temperature increases beyond 350 K, significant differences between the original calibration and the experimental data are seen. For example, at 525 K, the measured data^{22} indicate a density of 1.7 g/cm^{3}, whereas the original calibration^{1} is below 0.2 g/cm^{3} at this temperature. These differences will significantly affect the HE energy density, which will in turn affect the detonation velocity and pressure. Likewise, the Hugoniot response due to flyer impact will drastically be affected by such large density inconsistencies. Also, for safety modeling, the isochoric pressurization due to confined heating of HE charges will be rather different between these two models.

The last parameter required for the EOS model is *C*. Once the other parameters are set, *C* determines the behavior at high pressures. In particular, it was originally introduced to keep the reactant Hugoniot from crossing the product Hugoniot in overdriven detonations; a crossing of Hugoniots implies that the products are more dense than the reactants, which is not expected in slightly overdriven detonations. From examination of the available Hugoniot data,^{17,18,27} a value of *C* = 0.34 is determined to adequately fit the observed Hugoniot. Figure 4 shows the original model, current model, and available experimentally measured Hugoniot points. Although the original and current models have quite different specific heats and thermal expansion properties, the Hugoniots are rather similar. This indicates that fitting to Hugoniot data, although very important, does not adequately constrain the complete thermodynamic EOS.

The calibrated set of parameters for PBX 9502 are summarized in Table I.

## IV. PROPERTIES OF THE MODEL

It is worth exploring other thermodynamic properties of the current model, so as to compare with experimental data as it becomes available for validation.

From the isobaric thermal expansion, the densities at a variety of initial temperatures are obtained and listed in Table II. The reason for listing these specific temperatures is due to the variety of dynamic experiments performed at these temperatures.^{21–23,28–31} For three initial temperatures, the *U _{s}* –

*U*relations are shown in Fig. 5. As expected, the colder material has a stiffer response than the warmer material. The differences between the ambient and cold material in the

_{p}*U*–

_{s}*U*space are hard to measure experimentally,

_{p}^{28}and the results here are not in contradiction with the available data.

Temperature (K) . | Temperature (°C) . | Density (g/cm^{3})
. |
---|---|---|

77 | –196.15 | 1.942 |

218.15 | –55 | 1.914 |

297 | 23.85 | 1.890 |

348.15 | 75 | 1.869 |

403.15 | 130 | 1.840 |

525 | 251.85 | 1.697 |

Temperature (K) . | Temperature (°C) . | Density (g/cm^{3})
. |
---|---|---|

77 | –196.15 | 1.942 |

218.15 | –55 | 1.914 |

297 | 23.85 | 1.890 |

348.15 | 75 | 1.869 |

403.15 | 130 | 1.840 |

525 | 251.85 | 1.697 |

It is also worth examining the temperature along the Hugoniots for variations in the initial temperature. It is known that the HE sensitivity to shock depends quite strongly on the HE's initial temperature.^{22,28–31} Figure 6 shows the predicted shock temperature along the Hugoniot as a function of pressure for material at three different initial temperatures. It is observed that the temperature variation, between the initially cold material and hot material, increases with shock pressure. Future experimental efforts are being directed to measure the temperature of shocked PBX 9502.^{32}

Another interesting property of the model is to compute the isochoric pressure increase due to thermal heating. HEs generally expand thermally much more than metals. As such, for the HE encased in metal without sufficient ullage, it is an important safety characteristic to know how the HE will pressurize internally. A rough approximation is to examine the isochoric pressurization of the HE from Eq. (2). This is shown in Fig. 7 for the current model. As seen, there can be significant pressurization of the HE, even prior to decomposition. The pressures reached after heating to 500 K are comparable to yield strengths in many metals.

## V. VALIDATION OF THE MODEL

In this section, two validation tests are performed on the current model. Comparisons with isentropic and isothermal compression experiments are explored.

### A. Isentropic compression

By loading an HE sample via a continuous pressure pulse, isentropic compression of the material can be achieved, at least for a limited thickness of the sample. The Sandia National Laboratory Z-machine^{33} is one such facility to isentropically compress solid materials to high pressures, and provides useful data for validating the current EOS model. Comparisons with reactive flow models of two such isentropic compression experiments have previously been made.^{34} The computational results presented here follow the Lagrangian methodology outlined previously.^{35,36}

Shot Z1489 (Ref. 37) compressed PBX 9502 up to ∼5 GPa. A schematic of the experiment is shown in Fig. 8. The main diagnostic is the measured interfacial velocity between the sample HE and a LiF window. Figure 9 shows the results of the experiment and the present model. It is seen that the model “drive” matches the measured drive by construction.^{38} For the various PBX 9502 sample thicknesses, the model does a reasonable job predicting the experimentally observed traces, although the model is a bit too stiff at these low pressures, resulting in early interface motion. This stiffness at relatively low pressures is also observed in the Hugoniot in the range 0.2 mm/*μ*s < *U _{p}* < 0.6 mm/

*μ*s in Fig. 4. It is difficult for an analytical EOS to reproduce the “kink” observed in the Hugoniot from Fig. 4, without over- and under-shoots near the kink.

Shot Z1981 compressed PBX 9502 to ∼20 GPa. See Fig. 10 for comparison between the experimental and model interface velocities at various HE sample thicknesses. The current model matches the experiment quite well in both timing and amplitude of the waves. Again, it should be noted that the model results presented for these isentropic compression experiments are predictions and were not part of the calibration procedure.

It should be noted in both Z1489 and Z1981 that there is early motion not well modeled by the simple fluid EOS utilized here. In particular, at very low pressures of less than 0.1 GPa, and particle velocities less than 0.025 mm/*μ*s, the Hugoniot elastic limit^{17} has not been reached and strength effects of the HE are important. In reality, the composite and porous nature of the material at these low pressures are also important aspects and such modeling is outside the scope of the present study.

### B. Isothermal compression

Isothermal compression of the PBX 9502 constituents has been investigated.^{20,39,40} Earlier works are cited in Refs. 20 and 40. Although discrete isothermal compression data points on TATB and Kel-F have been taken, there is no straightforward method to combine these discrete data to obtain the mixture's isothermal response. A simple methodology entails fitting the discrete data with an analytical fitting form for each of the constituents and determining the response to the overall mixture. One particularly convenient analytical form is from Murnaghan^{41}

where *K _{T}* and $KT\u2032$ are the isothermal bulk modulus at zero pressure and its derivative with respect to pressure, respectively. This particular fitting form allows analytical inversion to obtain

*ρ*(

*p*), which is convenient for obtaining the mixture response. From the data, for TATB, it is found that

*K*= 5.82 GPa and $KT\u2032=1.85$ with

_{T}*ρ*

_{0}= 1.937 g/cm

^{3}. For Kel-F 800, the parameters

*K*= 5.77 GPa and $KT\u2032=3.26$ with

_{T}*ρ*

_{0}= 1.998 g/cm

^{3}are found. Note that the bulk modulus and its derivative are quite sensitive to the particular fitting form used,

^{20}and these parameters represent a global fit to the available data. See the blue curve in Fig. 11 for the resulting mixture PBX 9502 isothermal compression based on the Murnaghan fit to the experimental constituent data. Note that at

*p*= 0, the mixture is assumed to be at theoretical maximum density, which in this case is ∼1.94 g/cm

^{3}.

For the current Davis EOS model calibrated in Sec. III, the isothermal compression curve, from Eq. (2), at *T *=* T*_{0} is shown as the red curve in Fig. 11. At *p* = 0, the density is 1.890 g/cm^{3}. The discrepancy in initial density is due to the fact that PBX 9502 is porous, with ∼2.5% void space. Upon compression above the yield point, one expects the voids to be removed, and should be comparable to the initially void-free theoretical maximum density mixture. As can be seen in the figure, the isothermal compression of the current model quickly rises in density under applied compression, and then tapers off to a more gradual rise in density. The current model and Murnaghan mixture of constituents agree in density within ∼2% for the currently available data up to ∼45 GPa, which is not unreasonable given the uncertainties in the experimental data, Murnaghan fit to that data and the limitations of the Davis analytical EOS model.

## VI. CONCLUSIONS AND FUTURE WORK

Although the current model is quite adept at modeling a wide range of thermodynamics, there is room for further improvements. Specifically, this model (and many others) could benefit from a better underlying functional form for the specific heat at constant volume as a function of temperature, as well as the functional form of the Grüneisen parameter. Also, it would be beneficial to model porosity explicitly, as it is a key driver in shock initiation of HE. One may also wish to examine both from a modeling and experimental point of view any possible phase changes under applied pressure/temperature loading of TATB.

Experimentally, measuring shock temperatures, isentropic and isothermal compression at different initial temperatures, and thus effectively measuring isochoric pressurization are worthy pursuits for validation and testing of the model.

## ACKNOWLEDGMENTS

The author is grateful to a number of experimentalists who readily shared their results, specifically, Raja Chellappa, Dana Dattelbaum, Richard Gustavsen, Larry Hill, Shawn McGrane, and Darla Thompson. The author would also like to acknowledge John Bdzil, Carlos Chiquete, William Davis, Nicolas Desbiens, Ralph Menikoff, James Quirk, Mark Short, Remy Sorin, D. Scott Stewart, and Bradley Wescott for fruitful discussions on the modeling aspects of this work. John Bdzil, Ralph Menikoff, and Matthew Price are acknowledged for proofreading the manuscript and making many useful suggestions. This work was performed under the auspices of the U.S. Department of Energy National Nuclear Security Administration under Contract No. DE-AC52-06NA25396.

## References

^{∘}C