Free Energy Based Equation of State for Pentaerythritol Tetranitrate

Subscriber access provided by University of Otago Library

Article

Free Energy-Based Equation of State for Pentaerythritol Tetranitrate Marc J. Cawkwell, David S. Montgomery, Kyle J Ramos, and Cynthia A. Bolme J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b09284 • Publication Date (Web): 30 Nov 2016 Downloaded from http://pubs.acs.org on November 30, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Free Energy-based Equation of State for Pentaerythritol Tetranitrate M. J. Cawkwell,∗ D. S. Montgomery, K. J. Ramos, and C. A. Bolme Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States E-mail: [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract An equation of state for the energetic molecular crystal pentaerythritol tetranitrate (PETN) has been developed from a parameterized model for its Helmholtz free energy. The ion motion contribution to the free energy is represented by a sum of Debye models for the vibrational modes of mainly lattice phonon and intramolecular character. The dependence of the frequencies of the normal modes on density is captured using the quasi-harmonic approximation whereby the Debye temperatures for both populations of modes depends explicitly on specific volume. The dependence of the Debye temperatures on specific volume was parameterized to normal mode frequencies computed from solid state dispersion-corrected density functional theory. The model provides a good description of the thermophysical properties of PETN. The equation of state has been applied to the calculation of thermodynamic states along the principal Hugoniot of single crystal PETN.

Introduction Pentaerythritol tetranitrate, C5 H8 N4 O12 , (PETN) is a relatively sensitive secondary explosive. It adopts a body centered tetragonal unit cell containing two molecules with space group P ¯421 c under ambient conditions. 1 Another polymorph, PETN-II, has been observed at temperatures approaching melt at atmospheric pressure. The latter forms an orthorhombic unit cell containing four molecules with space group P cnb. 2 An orthorhombic structure with space group P 21 21 2 has been observed under static high pressure conditions around 6.3 GPa. 3 The structure of the PETN-III phase is closely related to the tetragonal polymorph that is stable under ambient conditions. Density functional theory calculations suggest PETN-III is only stable under non-hydrostatic pressure conditions. 3 Hence, it might be accessible in shock experiments on single crystals that are oriented appropriately. Strong impacts on energetic materials generate high pressure, high temperature conditions from which chemical reactions initiate. PETN is a particularly interesting material 2


Page 2 of 23

Page 3 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in this respect because its impact sensitivity is known to depend on the crystallographic orientation of the shock propagation direction. 4,5 The locus of thermodynamic states accessible by shock compression, the Hugoniot, is important for understanding the behavior of materials under impact. Pressure and density on the principal Hugoniot can be derived from experimental measurements of the shock wave speed, Us , for a given particle velocity, Up , through the Rankine-Hugoniot jump conditions. However, these measurements do not allow the temperature of the shocked state to be evaluated directly. Information on the temperature of shocked explosives is vital since chemical reactions in energetic materials are driven by combined high temperature, high pressure conditions. Continuum and mesoscale models for shock processes in materials typically rely on an equation of state to relate state variables. An accurate representation of the heat capacity by the equation of state is particularly important for organic materials since their heat capacities under ambient temperature and pressure are significantly smaller than the classical limit and are sensitive to temperature. For instance, models using a constant value for the heat capacity equal to the value at room temperature or the classical limit will tend to overor underestimate adiabatic heating, respectively. These sources of error in temperature estimates are quantified in Ref. 6 which provided a comparison of the shock heating attained in the energetic molecular crystal cyclotrimethylene trinitramine (RDX) from classical molecular dynamics simulations of uniaxial shock compression with an estimate derived from an equation of state with a constant, experimentally derived value for the heat capacity. A number of molecular simulation techniques have been applied to the calculation of Hugoniot states. These span the direct simulation of plate impacts 7 to molecular dynamics and Monte Carlo approaches to the solution of the Hugoniot relation, 1 E − E0 = (P + P0 )(V0 − V ), 2

(1)

where E is the total kinetic plus potential energy, P the pressure, and V the volume. 8–14 The

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

variables with subscripts denote the properties of the system in its pre-shocked state. The physical accuracy of the Hugoniot states derived from these methods depends sensitively on the quality of the interatomic potential that is used. We have developed an equation of state for PETN based on a parameterized model for its Helmholtz free energy that captures accurately the temperature dependence of its heat capacity. The Helmholtz free energy is represented using Debye theory within the quasiharmonic approximation. 15–17 Separate Debye models were constructed for the vibrational modes of mainly lattice phonon character and those of mainly intramolecular character. 18 The vibrational normal modes and their dependence on specific volume was determined using dispersion-corrected density functional theory calculations. The static, zero temperature lattice energy was parameterized such that complete model reproduces isothermal compression data from experiment. Physical properties derived from the free energy including the heat capacity, Gr¨ uneisen parameter, bulk modulus, and volumetric thermal expansion coefficient are in accord with values derived from experimental measurements. The free energy model has been applied to the computation of the principal Hugoniot up to 12 GPa where good accord with experimental data is found over the entire pressure range. We also present the dependence of temperature on the shock pressure.

Free energy model Our model for the Helmholtz free energy of a molecular crystal was described in detail in Ref. 18 Only the essential parts of the framework are outlined here. The total free energy, F (V, T ), is separated into a static lattice energy, Φ(V ), that depends only on the specific volume, V , and a term arising from the harmonic, thermal motion of the ions,

F (V, T ) = Φ(V ) + Fion (V, T ),

4


(2)

Page 5 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where T is the temperature. The ion motion depends on the spectrum of normal mode frequencies of the crystal and their dependence on the specific volume, 3(N −1) 3(N −1) X 1 X Fion (V, T ) = h ¯ ωi (V ) + kB T ln(1 − exp(−¯hωi (V )/kB T )) 2 i=1 i=1

(3)

where h ¯ is the reduced Planck constant, kB Boltzmann’s constant, N the total number of atoms, and ωi (V ) the angular frequency of normal mode i of a crystal with specific volume V . 15 Zerilli and Kuklja employed a similar approach for the development of a complete equation of state for β-cyclotetramethylene tetranitramine. 19 The full set of 3(N − 1) normal modes is partitioned into a set of mainly lattice phonon character and a set of intramolecular character. The former control the thermal pressure whereas the latter lead to the large and temperature-dependent heat capacity that is a characteristic of organic molecular crystals. The Q normal modes of lowest frequency are assigned to the population of lattice phonons by inspection of the vibrational density of states. The remaining 3(N − 1) − Q modes are assigned to the intramolecular vibrations. We replace the Fion terms for each population of normal modes by a Debye function,

FD = 3Nvib

"

1 kB θ D + kB T 8

T θD

3 Z

θD /T

2

z ln(1 − e )dz 0

−z

#

(4)

where θD the Debye temperature, and Nvib is the number of normal modes in Eq. 3 that the Debye function. In this work, Nvib = Q for the lattice phonons and 3(N − 1) − Q for the intramolecular modes. The Debye temperatures for the population of lattice phonons and intramolecular vibrations are determined numerically by minimizing the difference in the free energy from Eq. 3 and the corresponding Debye function, Eq. 4. Hence, the Debye temperature is that which minimizes, χ2 = (Fion (V, T ) − FD (V, T, θD ))2

5


(5)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

In Ref. 18 we determined the Debye temperatures by minimizing an objective function written in terms of the difference between the heat capacities calculated directly from the normal modes and the corresponding Debye models. The advantage of the new approach in Eq. 5 is that it is significantly faster to compute. The dependences of the Debye temperatures for the phonons and intramolecular vibrations on specific volume, θD,1 (V ) and θD,2 (V ), respectively, are represented analytically by,

θD (V ) = θD (V0 )

V V0

−γ ∞

exp −A

V V0

B −1 − 2

V V0

2

−1

!!

,

(6)

where γ ∞ = 2/3 is the Gr¨ uneisen parameter at infinite compression, V0 the specific volume of a reference state, and A and B are adjustable parameters. 16,17 The static lattice energy, or cold curve, Φ(V ), is represented by the Vinet equation of state, 20 Φ(V ) =

4V ∗ B ∗ [1 − (1 + X) exp(−X)] , (B1∗ − 1)2

(7)

where V ∗ , B ∗ , and B1∗ are adjustable parameters, and 3 X = (B1∗ − 1) 2

"

V V∗

1/3

−1

#

(8)

The cold curve is parameterized by requiring that the complete free energy model reproduces experimental isothermal compression data, Pref (V, T ). The contribution to the pressure from the cold curve is, P0 (V ) = Pref (V, T ) − PD,1 (V, T ) − PD,2 (V, T )

(9)

where PD,1 and PD,2 are the pressures derived from the Debye functions for the lattice phonons and intramolecular vibrations on the isotherm, respectively. The cold curve is then parameterized by a least squares fit of

P0 (V ) =

−4V ∗ B ∗ ∂X X exp(−X), (B1∗ − 1)2 ∂V 6


(10)

Page 7 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where, ∂X B1∗ − 1 = ∂V 2(V 2 V ∗ )1/3

(11)

to the right hand side of Eq. 9.

Application to PETN The vibrational normal modes of PETN for compressions 1 ≤ V /V0 ≤ 0.6 were calculated using density functional theory with the D3(BJ) dispersion correction of Grimme et al. 21 Empirical dispersion corrections provide additional weak cohesion to compensate for the underestimation of van der Waals-like interactions in regular density functional theory. As such, the use of dispersion corrections significantly improves the predicted equilibrium specific volume of molecular crystals. 22 All of the calculations reported here in addition used the generalized gradient approximation exchange-correlation functional of Perdew, Burke, and Ernzerhof (PBE), 23 the DZVP basis set, and Geodecker-Teter-Hutter pseudopotentials 24 as implemented in the cp2k package. 25 A plane wave cut-off of 1200 Ry and a relative cut-off of 60 Ry provided well converged total energies and a suitable distribution of Gaussian basis functions across the integration grids. Since the PETN unit cell is relatively small, containing only two molecules, all of our calculations were performed on a 2 × 2 × 2 supercell, containing 464 atoms, to reduce finite size effects. The crystal structure published by Conant et al. 1 was used as the starting configuration for the calculations. The internal coordinates and lattice parameters of the supercell were relaxed at zero external pressure until the force acting on any atom was less than 5.14 × 10−7 eV/˚ A. The predicted lattice parameters at zero pressure are a = b = 9.378 ˚ A and c = 6.574 ˚ A. These compare well with experimental measurements of a = b = 9.378 ˚ A and c = 6.708 ˚ A. 1 The relaxed supercell at zero pressure was compressed uniformly in 21 increments up to a compression of V /V0 = 0.6. The atomic positions in each compressed supercell were relaxed until the maximum force acting on any atom was less than 5.14 × 10−7 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

eV/˚ A. The vibrational normal modes were then calculated by numerically evaluating the force constant matrix within cp2k. The uniform compression of the tetragonal PETN supercell, where the a : c ratio of the PETN unit cell remains constant during compression, leads to non-hydrostatic stress conditions. This approach is required in the application of the general expansion of the thermoelastic free energy developed in Ref. 26 In that work, Luscher et al. decomposed the thermoelastic free energy into a term that depends only on the volumetric strain and a term that depends on coupled volumetric and deviatoric distortions. The former can be replaced by a complete, free energy-based equation of state for isotropic distortions whereas the latter accounts for anisotropic distortions and depends on tensors of second order elastic constants and anisotropic thermal expansivity. 26 Hence, including the effects of non-uniform distortions into the equation of state by using hydrostatic rather than isotropic compression in its construction would lead to small ’double counting’ errors when it is used in the complete, anisotropic free energy expansion. The vibrational density of states, n(f ), of the PETN supercell at zero pressure is presented in Fig. 1. The vibrational densities of state of the 21 compressed PETN supercells are provided in the Supplemental Information. The normal mode frequencies are in good agreement with those presented in Ref. 27 The inset of Fig. 1 highlights the continuous band of lattice phonons from zero frequency to about 150 cm−1 . In accord with our free energybased equation of state for RDX, we have assigned all normal modes with frequencies less than 210 cm−1 in the vibrational density of states evaluated at zero pressure and strain to the sub-set with mainly lattice phonon character. 18 The value of the cut-off frequency between the two populations is somewhat arbitrary owing to the coupled inter- and intramolecular nature of the normal modes of crystals of flexible molecules such PETN. The selection of the cut-off frequency was also guided by identifying a relatively wide gap in the computed vibrational density of states around the target frequency range. Based on these arguments, in this work we have assigned Q = 349 for the 2 × 2 × 2 PETN supercell (3N = 1392). The

8


Page 8 of 23

Page 9 of 23

349 modes of lowest frequency are assigned to the lattice phonons for all of the compressed supercells.

2 1

n(f)

-1

Density of states (modes/cm /unit cell)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 0

0 0

1000

0

100

200 -1 f (cm )

2000 -1 Frequency (cm )

300

400

3000

Figure 1: Phonon density of states at zero pressure calculated with dispersion-corrected density functional theory. The vertical broken line in the inset depicts the frequency cut-off between the modes of mainly lattice phonon character and intramolecular character. The dependences of the Debye temperatures on specific volume were calculated from the normal mode frequencies via Eq. 5 and are presented in Fig. 2. As expected, the Debye temperature for the low frequency lattice phonons depends sensitively on density whereas the Debye temperature for the intramolecular vibrational modes depends only weakly on density. In the parameterization of Eq. 6 we set the reference specific volume equal to that given by our density function theory calculation at zero pressure, i.e., V0 = 0.550683 cm3 /g. The remaining terms in Eq. 6, namely θD (V0 ), A, and B were determined by a least squares fit. These values are provided for both Debye temperatures, θD,1 and θD,2 , in Table 1. 9



2.0 θD(V) / θD(V0)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

Phonons

1.5

Intramolecular vibrations

1.0 0.40 0.45 0.50 3 Specific volume (cm /g)

0.35

0.55

Figure 2: Dependence of the Debye temperatures for the lattice phonons and intramolecular vibrations on specific volume.

Table 1: Parameterization of Eq. 6 for the dependence of the two Debye temperatures, θD,1 and θD,2 , on specific volume. θD,1 θD,2

θD (V0 ) (K) 211.531 2288.19

A 1.34571 -0.604418

10

B -0.480332 0.003960763


Page 11 of 23

The cold curve, Φ(V ), was parameterized to the isothermal compression data of Olinger, Halleck, and Cady 28 by a least squares fit of Eq. 10 after subtracting the contribution to the total pressure from the Debye functions via Eq. 9. The parameters for the Vinet equation of state, Eq. 7, are V ∗ = 0.530640 cm3 /g, B ∗ = 14.1217 GPa, and B1∗ = 9.19302. The room temperature isotherm of PETN computed from our free energy model and the experimental data of Olinger et al. is presented in Fig. 3

12 10 Pressure (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8 6 4 2 0 1.8

2.0 2.2 3 Density (g/cm )

2.4

Figure 3: Isothermal compression data for PETN at 300K. The solid line was derived from the equation of state and the symbols are experimental data from Ref. 28

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

PETN Properties and Hugoniot The heat capacity at constant volume, CV , is derived from the heat capacity at constant pressure, CP , via CV = CP − T V α 2 B T ,

(12)

where α is the coefficient of volumetric thermal expansion and BT the isothermal bulk modulus. There exists some uncertainty in the appropriate value of the isothermal bulk modulus of PETN. Thermodynamic stability requires that the isothermal bulk modulus is less than the isentropic bulk modulus, i.e., BT < BS . The isothermal bulk modulus derived from a least squares fit of the Vinet equation of state, Eq. 7, to the isothermal compression data of Olinger et al. yields BT = 10.41 GPa whereas the Voigt average isentropic Bulk modulus derived from the single crystal elastic constants tabulated by Sun et al. yields BS = 9.70 GPa. 29 Olinger et al. also report the results of earlier experiments by Vasil’ev that focus on the pressure range 0 - 2 GPa. 30 The latter give BT ≈ 8.3 GPa at zero pressure, in accord with the thermodynamic stability condition. Using the value for the isothermal bulk modulus derived from the measurements of Vasil’ev with experimental data for CP , 31 α, 32 and the specific volume under ambient conditions 1 gives CV = 325.5 J/mol/K. The heat capacity at constant volume in terms of the set of normal mode frequencies is, 3(N −1)

C V = kB

X (¯hωi /kB T )2 exp (−¯hωi /kB T ) . 2 (exp (−¯ h ω /k T ) − 1) i B i=1

(13)

Using the normal mode frequencies calculated at zero pressure from dispersion-corrected density functional theory and T = 300 K we find CV = 324.7 J/mol/K, which is in outstanding agreement with the value derived from experimental measurements. The total heat capacity at constant volume from the sum of Debye functions is,

CV = 3QkB

T θD,1

Z

θD,1 /T 0

z 4 ez dz + 3(3N − Q)kB (ez − 1)2

12

T θD,2


Z

θD,2 /T 0

z 4 ez dz. (14) (ez − 1)2

Page 13 of 23

At T = 300 K and the equilibrium density at zero pressure, the heat capacity from the Debye functions, Eq. 14, is CV = 265.2 J/mol/K. This value underestimates experiment and the value computed directly from the normal mode frequencies using Eq. 13 by about 19%. This discrepancy is consistent with that seen in our equation of state for RDX. The temperature dependences of the heat capacities derived from the frequencies of the normal modes, Eq. 13, and the Debye models, Eq. 14, are presented in Fig. 4. The heat capacity derived from a single Debye function fitted to the normal mode data using Eq. 5 is plotted in Fig. 4 for comparison with the extended model. The extended model with two Debye functions clearly provides a better description of the heat capacity of PETN at low temperatures than the model with a single Debye function.

800

600 CV (J/mol/K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


400

200

0 0

Exact θD = 1601.4 K θD,1 = 204.8 K, θD,2 = 2285.3 K Experiment

1000 Temperature (K)

2000

Figure 4: Heat capacity of PETN at constant volume. The exact curve was calculated directly from the set of normal modes at zero pressure using Eq. 13 whereas the broken lines were calculated using Debye models with one and two Debye functions. The value of CV derived from experimental data is also included. 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

The experimental value of the coefficient of volumetric thermal expansion, α, of PETN is 2.355 × 10−4 K−1 at room temperature. 32 The value derived from our free energy model, α = 1.547 × 10−4 K−1 , underestimates the experimental data by about 52%. Our free energy model for RDX gave an almost identical relative error in the thermal expansion coefficient. 18 The isothermal and isentropic bulk moduli derived from the free energy model both overestimate experimental data. At room temperature and zero pressure we find BT = 9.8 GPa and BS = 10.3 GPa. The bulk moduli are controlled mainly by the static lattice energy term in the free energy. The cold curve, Eq. 7, was parameterized to the isothermal compression data of Olinger et al. However, since the data from Ref. 28 is relatively sparse at small compression, the bulk modulus at ambient conditions is not captured adequately. Indeed, Olinger et al. used low pressure isothermal compression data from Ref. 30 to estimate BT . The Gr¨ uneisen parameter, γ=

V αBT , CV

(15)

of PETN derived from room temperature experimental data is γ = 1.07. 28 The Gr¨ uneisen parameter computed from our free energy model is γ = 1.03 and hence compares very favorably with experiment. Nevertheless, the good accuracy of the predicted value is somewhat fortuitous owing to a cancellation of errors between the predicted volumetric thermal expansion, isothermal bulk modulus, and heat capacity. The Gr¨ uneisen parameter at zero pressure can also be estimated from the dependence of the Debye temperatures on specific volume, Eq. 6, via

γ = γ ∞ + A + B.

(16)

The contribution to the total Gr¨ uniesen parameter from the lattice phonons is large, γ = 1.532, since these vibrational modes control the thermal pressure. Similarly, the intramolecular vibrations contribute weakly to the thermal pressure of the crystal and the corresponding

14


Page 15 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Gr¨ uneisen parameters is only γ = 0.0662. As reported in Ref., 18 the total Gr¨ uneisen parameter lies between the contributions from the phonons and intramolecular vibrations. Adiabatic shock heating on the principal, or single shock Hugoniot is an important quantity for energetic materials since the first reactions are thought to be driven mainly by thermal activation. Furthermore, the capability to compute states on the Hugoniot are useful for the design and interpretation of shock compression experiments. The principal Hugoniot of PETN under initial conditions of zero pressure, density ρ0 = 1.778 g/cm3 and T = 300 K has been computed from the free energy model. In Fig. 5 we present the theoretical Hugoniot together with experimental data for shocks on PETN single crystals from Ref. 33 The agreement between the theoretical and experimental data is very good over the pressure interval 0 < P < 12 GPa. Nadykto proposed that a compression-induced change in the electronic structure of PETN gives rise to a change in the gradient of the isotherm and Hugoniot at about 4.15 GPa. 34 The response of PETN to compression was modeled using different equations of state for conditions above and below the transition. However, the complete, free energy-based equation of state developed in this work and the incomplete equations of state developed in Refs. 35,36 demonstrate that the isotherm and Hugoniot of single crystal PETN can be described adequately by assuming a single phase. The good agreement between the predicted and experimental Hugoniot data in Fig. 5 allows us to apply our free energy model to shock problems with confidence. We provide theoretical predictions of the temperature along the PETN Hugoniot in the range 0 < P < 12 GPa in Fig. 6. Over the studied pressure range the dependence of the shock heating, ∆T , on pressure is described accurately by the polynomial ∆T = 0.552488P 2 + 24.1506P .

15



12 Free energy: ρ0 = 1.778 g/cm Single crystal experiments

10 Pressure (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

3

8 6 4 2 0

1.8

1.9

2.0 2.1 2.2 3 Density (g/cm )

2.3

Figure 5: Principal Hugoniot for PETN for ambient initial conditions (ρ0 = 1.778 g/cm3 , T = 300 K) computed from the free energy with experimental Hugoniot data for PETN single crystals from Ref. 33

16


Page 17 of 23

700

600 Temperature (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


500

400

300 0

2

4

8 6 Pressure (GPa)

10

12

Figure 6: Temperature as a function of shock pressure on the PETN principal Hugoniot calculated from the free energy model. The initial conditions were ρ0 = 1.778 g/cm3 and T = 300 K.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Summary A model for the Helmholtz free energy of energetic molecular crystal PETN has been developed based on Debye theory. Separate Debye functions for the lattice phonons and intramolecular vibrations have been parameterized to normal mode frequencies calculated using dispersion-corrected density functional theory. The dependence of the Debye temperatures on specific volume is captured with the quasi-harmonic approximation. The static lattice energy was parameterized to high quality isothermal compression data. The free energy model provides a good description of the thermophysical properties of PETN. The Gr¨ uneisen parameter and the temperature dependence of the heat capacity are captured particularly accurately. Finally, the free energy model was used to calculate the principal Hugoniot of PETN. Experimental data confirms the accuracy of the theoretical predictions. The shock heating along the principal Hugoniot of single crystal PETN was also obtained

Acknowledgments This work was supported by the National Nuclear Security Administration Science Campaign 2 and performed at Los Alamos National Laboratory under DE-AC52-06NA25396. M.J.C. thanks Jeff Leiding and Scott Crockett for helpful discussions.

References (1) Conant, J. W.; Cady, H. H.; Ryan, R. R.; Yarnell, J. L.; Newsam, J. M. The Atom Positions of Pentaerythritol Tetranitrate (PETN, C5H8N4O12) Determined by X-ray and by Neutron Diffraction (LA-7756-MS). 1979. (2) Cady, H. H.; Larson, A. C. Pentaerythritol tetranitrate II - Its crystal structure and transformation to PETN I - An algorithm for refinement of crystal structures with poor data. Acta Crystallogr. 1975, B41, 1864. 18


Page 18 of 23

Page 19 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(3) Tschauner, O.; Kiefer, B.; Lee, Y.; Pravica, M.; Nicol, M.; Kim, E. Structural transition of PETN-I to ferroelastic orthorhombic phase PETN-III at elevated pressures. J. Chem. Phys. 2007, 127, 094502. (4) Dick, J. J. Effect of crystal orientation on shock initiation sensitivity of pentaerythritol tetranitrate explosive. Appl. Phys. Lett. 1984, 44, 859. (5) Dick, J. J.; Ritchie, J. P. Molecular mechanics modeling of shear and the crystal orientation dependence of the elastic precursor shock strength in pentaerythritol tetranitrate. J. Appl. Phys. 1994, 76, 2726. (6) Cawkwell, M. J.; Sewell, T. D.; Zheng, L. Q.; Thompson, D. L. Shock-induced shear bands in an energetic molecular crystal: Application of shock-front absorbing boundary conditions to molecular dynamics simulations. Phys. Rev. B 2008, 78, 014107. (7) Holian, B. L.; Lomdahl, P. S. Plasticity induced by shock waves in nonequilibrium molecular dynamics simulations. Science 1998, 280, 2085. (8) Erpenbeck, J. J. Molecular dynamics of detonation. I. Equation of state and Hugoniot curve for a simple reactive fluid. Phys. Rev. A 1992, 46, 6404. (9) Brennan, J. K.; Rice, B. M. Efficient determination of Hugoniot states using classical molecular simulation techniques. Mol. Phys. 2003, 101, 3309–3322. (10) Maillet, J.-B.; Mareschal, M.; Soulard, L.; Ravelo, R.; Lomdahl, P. S.; Germann, T. C.; Holian, B. L. Uniaxial hugoniostat: a method for atomistic simulations of shocked materials. Phys. Rev. E 2001, 63, 016121. (11) Ravelo, R.; Holian, B. L.; Germann, T. C.; Lomdahl, P. S. Constant stress hugoniostat. Phys. Rev. B 2004, 70, 014103. (12) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. A method for tractable dynamical studies of single and double shock compression. Phys. Rev. Lett. 2003, 90, 235503. 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(13) Coe, J. D.; Sewell, T. D.; Shaw, M. S. Nested Markov chain Monte Carlo sampling of a density functional theory potential: Equilibrium thermodynamics of dense fluid nitrogen. J. Chem. Phys. 2009, 131, 074105. (14) Desbiens, N.; Bourasseau, E.; Maillet, J.-B.; Soulard, L. Molecular based equation of state for shocked liquid nitromethane. J. Hazardous Mater. 2009, 166, 1120–1126. (15) Wallace, D. C. Thermodynamics of Crystals; Dover, 1998. (16) Greeff, C. W. Phase changes and the equation of state of Zr. Modelling and Simulation in Materials Science and Engineering 2005, 13, 1015–1027. (17) Greeff, C. W.; Graf, M. J. Lattice dynamics and the high-pressure equation of state of Au. Physical Review B - Condensed Matter and Materials Physics 2004, 69, 054107. (18) Cawkwell, M. J.; Luscher, D. J.; Addessio, F. L.; Ramos, K. J. Equations of state for the alpha and gamma polymorphs of cyclotrimethylene trinitramine. J. Appl. Phys. 2016, 119, 185106. (19) Zerilli, F. J.; Kuklja, M. M. Ab initio equation of state of the organic molecular crystal: beta-Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine. J. Chem. Phys. A 2010, 114, 5372–5376. (20) Vinet, P.; Ferrante, J.; Smith, J. R.; Rose, J. H. A universal equation of state for solids. J. Phys. C 1986, 19, L467. (21) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. Journal of Computational Chemistry 2011, 32, 1456– 1465. (22) Sorescu, D. C.; Rice, B. M. Theoretical Predictions of Energetic Molecular Crystals at Ambient and Hydrostatic Compression Conditions Using Dispersion Corrections to Conventional Density Functionals (DFT-D). J. Phys. Chem. C 2010, 114, 6734–6748. 20


Page 20 of 23

Page 21 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(23) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (24) Goedecker, S.; Teter, M.; Hutter, J. Seperable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703. (25) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Computer Phys. Commun. 2005, 167, 103. (26) Luscher, D. J.; Addessio, F. L.; Cawkwell, M. J.; Ramos, K. J. A dislocation densitybased continuum model of the anisotropic shock response of single crystal alphacyclotrimethylene trinitramine. J. Mech. Phys. Solids 2017, 98, 63–86. (27) Velizhanin, K. A.; Kilina, S.; Sewell, T. D.; Piryatinski, A. First-principles-based calculations of vibrational normal modes in polyatomic materials with translational symmetry: Application to PETN molecular crystal. J. Phys. Chem. B 2008, 112, 13252. (28) Olinger, B.; Halleck, P. M.; Cady, H. H. The isothermal linear and volume compression of PETN to 10 GPa (100 kbar) and the calculated shock compression. J. Chem. Phys. 1975, 62, 4480. (29) Sun, B.; Winey, J. M.; Hemmi, N.; Dreger, Z. A.; Zimmerman, K. A.; Gupta, Y. M.; Torchinsky, D. H.; Nelson, K. A. Second-order elastic constants of pentaerythritol tetranitrate and cyclotrimethylene trinitramine using impulsive stimulated thermal scattering. J. Appl. Phys. 2008, 104, 073517. (30) Vasil’ev, M. Y.; Balashov, D. B.; Mokrousov, L. N. Isothermal compressibility of explosives under pressures up to 22000-kg-cm2. Russ. J. Phys. Chem. 1960, 34, 1159. (31) Baytos, J. F. Specific heat and thermal conductivity of explosives, mixtures, and plasticbonded explosives determined experimentally (LA-8034-MS). 1979. 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(32) Cady, H. H. Coefficient of thermal expansion of pentaerythritol tetranitrate and hexahydro-1,3,5-trinitro-s-triazine (RDX). J. Chem. Engin. Data 1972, 17, 369. (33) Marsh, S. P., Ed. LASL Shock Hugoniot Data; University of California Press, 1980. (34) Nadykto, B. A. Equations of state of unreacted explosives: PETN, RDX, and TATB of state. EPJ Web of Conferences 2010, 10, 00007. (35) Yoo, C. S.; Cynn, H.; Howard, W. M.; Holmes, N. Equations of state of unreacted high explosives at high pressures. Eleventh International Detonation Symposium. 1998; p 951. (36) Nagayama, K.; Kubota, S. Equation of state of energetic materials for unreacted compression. Proceedings of the International Symposium on Interdisciplinary Shock Wave Research. 2004; p 377.

22


Page 22 of 23

Page 23 of 23

0

Compression

Vibrational density of states

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39


1000

2000


3000 -1

Frequency (cm )

Free Energy Based Equation of State for Pentaerythritol Tetranitrate

Recommend Documents