Application of the PengâRobinson Equation of ... - ACS Publications

Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI

Article

Application of the Peng-Robinson equation of state to energetic materials RDX and TNT: pure components, liquid mixtures, and solid mixtures Philip Myint, Matthew McClelland, and Albert L. Nichols Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.5b04808 • Publication Date (Web): 01 Feb 2016 Downloaded from http://pubs.acs.org on February 4, 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.

Industrial & Engineering Chemistry Research 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 47

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

Industrial & Engineering Chemistry Research

Application of the Peng-Robinson equation of state to energetic materials RDX and TNT: pure components, liquid mixtures, and solid mixtures Philip C. Myint,∗,† Matthew A. McClelland,‡ and Albert L. Nichols III‡ †Design Physics Division, Lawrence Livermore National Laboratory, Livermore, CA, USA ‡Materials Science Division, Lawrence Livermore National Laboratory, Livermore, CA, USA E-mail: [email protected] Abstract Energetic materials are substances that can undergo rapid, exothermic reactions when subjected to an external stimulus, such as heating. In this work, we show that the well-known Peng-Robinson equation of state can be applied to energetic materials, whether they are pure components, liquid mixtures, or solid mixtures. We are specifically interested in two energetic materials: hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX) and 2,4,6-trinitrotoluene (TNT). We model RDX and TNT in both their liquid and solid phases, as well liquid and solid mixtures of the two compounds. Our work examines temperatures and pressures as high as about 500 K and 2500 bar, respectively. The Peng-Robinson equation of state provides a good representation of experimental volumetric (e.g., density and bulk modulus), thermal (heat capacity), and phase behavior (melting temperature and solubility) data. It can be applied to other energetic

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

materials, ranging in complexity from pure components to multicomponent mixtures, by adapting the procedures described in this study.

Introduction Energetic materials may be broadly defined as substances that can undergo rapid, exothermic reactions when subjected to an external stimulus, such as heating or mechanical impact. 1–3 This class of materials includes explosives, pyrotechnics, and propellants. Numerical simulations are increasingly being used as a means to study their behavior. One active area of research concerns cookoff models, in which energetic materials placed inside a sealed container are heated on timescales ranging from microseconds to days until they explode. 4–8 This area has applications to the behavior of explosives in naval shipboard fires and to the safe handling and storage of explosives in general. The heating raises the temperature of the energetic materials, causing them to eventually reach temperatures where chemical decomposition becomes significant. The formation of product gases increases the pressure, which further accelerates the reactions and may lead to an explosion that ruptures the container walls. The pressure typically reaches around 2000 bar before the explosion occurs, but this number can vary significantly since it depends on many factors, such as the rate of heating, the quality of the container’s seal, and the identity of the energetic material. Cookoff models, like virtually every other continuum model, rely on thermodynamic information to close the governing equations. The thermodynamic properties of the product gases released by energetic materials can be predicted well by thermochemical codes such as Cheetah. 9 In contrast, thermodynamic models of unreacted materials in their liquid and solid forms are still relatively undeveloped. Many authors have performed molecular simulations to study both liquid and solid energetic materials (see the reviews by Fried et al. 1 and by Politzer and Boyd, 2 as well as the references cited later in this study). The simulations can


Page 2 of 47

Page 3 of 47

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


be either directly or indirectly coupled with continuum models through various techniques. However, molecular simulations of energetic materials have been limited to single-component systems only, not mixtures. Macroscopic thermodynamic models therefore represent the only currently feasible way to model both pure components and mixtures of energetic materials. For the sake of discussion, we divide this broad class of models into two subclasses: 1) excess Gibbs energy models 10–14 and 2) equations of state based on departure functions, such as the various cubic, cubic-plus-association (CPA), and statistical associating fluid theory (SAFT) equations of state. 15–19 The central quantities of interest in the former subclass are the activity coefficients, while they are the fugacity coefficients in the latter subclass. The focus of this study is the Peng-Robinson equation of state (EOS), which is a cubic EOS. 20 Excess Gibbs energy models have enjoyed enormous success in predicting the phase equilibria of many systems, but their utility for liquid and solid mixtures of energetic materials may be more limited, particularly at higher pressures. The underlying reason is because experimental data in the field of energetic materials are much more scarce than they are for traditional systems, such as carbon dioxide (CO2 )-hydrocarbon mixtures. Experimental measurements, especially those at higher temperatures and pressures, are difficult to obtain because of the highly reactive nature of the materials. This has implications for excess Gibbs energy models, whose purpose is to determine the activity coefficients that relate the chemical potentials of the components in a mixture to their pure-component values. The chemical potentials of the pure components can be calculated (with respect to an arbitrary reference state) by integrating the constant-pressure molar heat capacity over a temperature interval. 13,14 The limits of integration involve temperatures where a phase change occurs, such as the melting point. This presents a serious difficulty for energetic materials because, as we demonstrate later, their melting points at the elevated pressures of interest in this study are not known to a high degree of precision. Furthermore, heat capacity vs. temperature data for some compounds (particularly liquids) are not available even 3 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

at atmospheric pressure, let alone at higher pressures. The uncertainties in these quantities are reflected in the pure-component chemical potentials and are consequently propagated to the activity coefficients and the mixture chemical potentials. These difficulties may be why excess Gibbs energy models for the phase behavior of energetic material mixtures are currently limited to a few studies at atmospheric pressure. 12,14 No such models exist for higher pressures, such as those relevant to cookoff near the explosion time. Equations of state based on departure functions can typically be expressed in a form that is applicable to both pure components and mixtures. The only difference between the two cases is whether the parameters of the EOS refer to a pure-component system or to a mixture. Since both types of systems are represented in a unified manner, these equations of state offer a more theoretically elegant alternative to the excess Gibbs energy approach described above. The mixture parameters are calculated from the pure-component parameters through a set of mixing rules that may involve one or more interaction coefficients. On a practical level, the equations of state are suited for energetic materials because the parameters can be fitted to experimental results that are more readily available than heat capacities. For example, the pure-component parameters can be calibrated with vapor pressure and volumetric (e.g., density) data, while the interaction coefficients of the mixtures can be adjusted to match solubility data. If vapor pressure information is not available, as is the case for one of the compounds in this study, we can use melting point and enthalpy of fusion measurements to infer the vapor pressure. An important detail we emphasize is that only atmospheric-pressure volumetric, melting point, or enthalpy of fusion data are necessary to calibrate the equations of state. This means that unlike the case with excess Gibbs energy models, the melting point as a function of pressure does not need to be known a priori to make predictions at higher pressures. Instead, the effect of pressure on the phase behavior (like the melting point) is something that we determine a posteriori with the equations of state. In this work, we show that the Peng-Robinson equation of state 20 can be applied to 4 ACS Paragon Plus Environment

Page 4 of 47

Page 5 of 47

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


energetic materials, whether they are pure components, liquid mixtures, or solid mixtures. The Peng-Robinson EOS is perhaps the most widely used cubic EOS in the world. Its success is such that some authors 21 have marveled at its “unreasonable effectiveness.” We are specifically interested in two energetic materials: hexahydro-1,3,5-trinitro-1,3,5-triazine and 2,4,6-trinitrotoluene. These two compounds are commonly abbreviated as RDX and TNT, respectively. Past studies have used the Peng-Robinson EOS to represent the gaseous products released by RDX and TNT during cookoff, 8 or to predict the solubility of RDX in CO2 . 22 But until now, there has not been an attempt to use this equation of state to model pure RDX and TNT in their condensed phases, as well as liquid and solid mixtures of the two compounds. This study is motivated by our research group’s recent efforts 6,7 to simulate the cookoff of Composition B, which is a mixture whose overall composition can be approximated as 64% RDX and 36% TNT by mass. Our work analyzes temperatures and pressures up to about 500 K and 2500 bar, respectively. These are conditions typically encountered in cookoff simulations before explosion occurs. The Supporting Information summarizes key expressions for the Peng-Robinson EOS that we reference repeatedly throughout this paper.

Pure components Solid RDX, solid TNT, and liquid TNT The next three subsections describe how we have applied the Peng-Robinson equation of state to model the thermodynamic properties of three different compounds: solid RDX, solid TNT, and liquid TNT. Our approach relies on vapor pressure and volumetric (density or thermal expansion) measurements, which are available for these three energetic materials. There is much less experimental data, including vapor pressure data, on liquid RDX. This is because pure liquid RDX is very reactive at the high-temperature conditions where it exists. The last subsection in this major section on pure components discusses how we have 5 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

Page 6 of 47

calibrated the Peng-Robinson EOS for liquid RDX with the limited data that is available. Determining the pure-component parameters ai and bi from vapor pressure and volumetric data The Supporting Information states that the parameters a and b for a mixture are calculated from the pure-component parameters ai and bi through the van der Waals mixing rules in (S2)–(S4), but it does not explain how ai and bi are obtained. For vapor-liquid equilibria, ai and bi at the critical point (Pc,i , Tc,i ) of component i must satisfy

ai (Tc,i ) = 0.45724 bi (Tc,i ) = 0.0778

2 R2 Tc,i , Pc,i

(1)

RTc,i , Pc,i

(2)

where R is the gas constant. These constraints follow 23,24 from (S5) and the requirements ∂P ∂v

!

T =Tc,i

=

∂2P ∂v 2

!

= 0,

(3)

T =Tc,i

in which P is the pressure, ρ is the molar density, and v = 1/ρ is the molar volume. At any other temperature T , it has been found that equations of the form

ai (T ) = 0.45724 bi (T ) = 0.0778

2 R2 Tc,i ψa (Tr,i ), Pc,i

(4)

RTc,i ψb (Tr,i ), Pc,i

(5)

can closely reproduce vapor-liquid equilibria measurements for a number of different systems. Here, Tr,i = T /Tc,i is the reduced temperature of i, and the functions ψa , ψb are chosen so that ψa (Tr,i = 1) = ψb (Tr,i = 1) = 1. Traditionally, ψb (Tr,i ) is chosen to be unity for all temperatures so that bi is independent of T (i.e., is a constant), while ψa (Tr,i ) is calculated from


Page 7 of 47

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


a general correlation involving the acentric factor ωi , which is related to the vapor pressure of i evaluated at Tr,i = 0.7. While this approach demonstrates the universality of the EOS, it can potentially cause systematic errors (5% or more) in the liquid-phase density predictions. One way to overcome the error in the density predictions is to apply volume translation. 25,26 Another possibility that Xu and Sandler 23,27 have explored is to allow ψb to be temperaturedependent and to calculate ψa and ψb from fluid-specific correlations fitted to experimental vapor pressure and density data. Martynov et al. 28 have followed a similar approach to show that the Peng-Robinson EOS, which has already been extensively applied to CO2 in the liquid and vapor phases, can also model solid CO2 . They report good agreement between their predictions of derived properties, such as the heat capacity, and experimental data. Except in the case of liquid RDX, we have not attempted to model ai and bi according to (4) and (5). One reason is because we are interested in modeling solids as well as liquids with the Peng-Robinson EOS, and the same set of ai and bi is not applicable to both the solid and liquid phases of a compound. The EOS can model only one condensed phase for a given set of ai and bi because if there are three real roots for the compressibility factor Z in (S9), only two of them are valid; one represents the vapor phase and the other represents the condensed phase. This means that solid RDX must be modeled by different parameters than liquid RDX, for example. Since the conditions in (3) are not relevant to solids, their parameters do not need to satisfy (1) and (2). For liquid TNT, we do not calculate ai and bi from critical properties because no such experimental data exists in the literature. In fact, due to the extremely reactive nature of energetic materials at higher temperatures, no direct experimental measurements of critical properties have been reported. A couple of studies have used theoretical relations to predict the critical properties. 29,30 Maksimov 29 presents critical property predictions for two energetic materials: RDX and HMX. Using group contribution methods, Toghiani et al. 30 are able to closely match Maksimov’s reported critical temperature for both RDX and HMX (840 K and 927 K, respectively). Their predicted 7 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

Page 8 of 47

critical pressure of 58 bar is, however, significantly larger than Maksimov’s value of 37 bar. This implies that their critical volume will differ from Maksimov’s result of 0.442 L/mole by roughly the same factor. Toghiani et al. also provide group contribution predictions of the acentric factor, but their results cannot be verified with experimental data because vapor pressure measurements 31,32 are available only up to temperatures as high as 450 K for RDX or 480 K for HMX, well below the value of Tr,i = 0.7 for these compounds. We determine the functional forms of ai and bi for solid RDX, solid TNT, and liquid TNT using vapor pressure and density (or other volumetric) data. Our formulation is similar to that of Martynov et al. 28 mentioned earlier, except that bi in our study does not depend on temperature. At every point (P, T ) along the vapor pressure curve of compound i, the condensed phase (superscript l, s) is in equilibrium with the vapor phase (superscript v). This equilibrium condition can be expressed in terms of the fugacity coefficients as

v φl,s i (T, P ) = φi (T, P ).

(6)

The fugacity coefficients are evaluated using (S14). The parameters ai and bi are related to their dimensionless analogues A and B in (S14) through (S7) and (S8). Both φl,s i (T, P ) and φvi (T, P ) are calculated from the same values of A and B, but φl,s i (T, P ) is a function of Z l,s = P/ρl,s RT , while φvi (T, P ) is a function of Z v = P/ρv RT . Martynov et al. 28 show that since Z l,s and Z v both satisfy (S9), this cubic polynomial equation can be written as

(Z − Z l,s )(Z 2 + qZ + r) = 0, where q = B − 1 + Z l,s , r = (A − 3B 2 − 2B) + (B − 1)Z l,s + Z l,s · Z l,s .


(7)

Page 9 of 47

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 compressibility factor Z v of the vapor phase is the larger root of the quadratic in (7): q Zv = − + 2

s

q2 − r. 4

(8)

The value of the constant bi is chosen to match density data of compound i, as we discuss in more detail below. Hence, the density data also determines B. This leaves only two unknowns — A and Z l,s — since Z v can be calculated from A, B, and Z l,s using (8). The two equations that we solve to find A and Z l,s are (S9) and (6). Vapor pressures are presented in two reviews 31,32 that have compiled and fit data from the literature to the Clausius-Clapeyron equation in the form

log10

P Bvap = Avap − , 0 P T

(9)

where P 0 = 1 is a reference pressure that determines the units of P . Equation (9) is valid because the vapor pressures of energetic materials tend to be very low (10−10 –10−2 bar, depending on T ) so that their vapors can be approximated as ideal gases. The values of Avap and Bvap that we use to determine ai and bi are presented in Table S4 of the Supporting Information. Our values for solid RDX are obtained from Östmark et al., 31 who have fitted vapor pressure measurements from five different references over the temperature range 310– 411 K. Östmark et al. also present results for solid and liquid TNT. For liquid TNT, they incorporate data only from Edwards, 33 while for solid TNT, they use data from Edwards in addition to several other references. Their two vapor pressure curves intersect at the point (P = 8.38 × 10−6 bar, T = 351.80 K), which represents the triple point of TNT predicted by the curves. The normal melting point of TNT will be at a temperature that is only slightly higher than the triple point temperature. Indeed, we have found that if we select the Avap and Bvap values provided by Östmark et al. for solid and liquid TNT, the normal melting point obtained from the Peng-Robinson EOS is 351.84 K. This is a little lower than 9 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

Page 10 of 47

the melting temperature cited in the literature, 34–38 which lies somewhere in the range 353– 354 K, depending on the reference. In contrast, the Avap and Bvap values from Edwards lead to a triple point of (P = 1.09 × 10−5 bar, T = 354.65 K), so the normal melting point we expect to obtain from his data may be a little high. In order to have a normal melting point in the desired range, we adopt the Avap and Bvap from Edwards, but slightly change the Avap of solid TNT from Edwards’ value of 14.34 to 14.35. This change in Avap leads to a triple point of (P = 9.97 × 10−6 bar, T = 353.62 K). It also raises the vapor pressure of solid TNT by a small amount (the average deviation over the range 293–353 K is 2.3%), which actually improves the agreement with the solid TNT vapor pressures from Östmark et al. Densities along the vapor pressure curves are not available for energetic materials, but due to their non-volatile (low vapor pressure) nature and the relatively weak dependence of ρl,s on pressure, we can approximate these densities with data at atmospheric pressure. For solid RDX and liquid TNT, we fix the value of bi to match the linear relations of density vs. temperature provided by Sun et al. 39 and Cady and Rogers, 40 respectively. We are unable to find density vs. temperature data for solid TNT, and instead estimate its density based on information regarding the thermal coefficient of expansion α, which is defined as 1 α= v

∂v ∂T

!

.

(10)

P

Vrcelj et al. 41 have found α at room temperature (293 K) for the more stable (i.e., monoclinic) form of TNT to be 2.06 × 10−4 K−1 , which agrees with the α measured by Heberlein, 42 but is higher than the values reported by Dobratz and Crawford 36 (1.8 × 10−4 K−1 ) and Sorescu et al. 43 (1.36 × 10−4 K−1 at 300 K). Over the temperature range 323–337 K, Rosen et al. 44 report α to be 3.03 × 10−4 K−1 . Based on these results, we assume that α varies linearly from the value obtained by Vrcelj et al. at 293 K to that reported by Rosen et al. at 330 K, which is the midpoint of the temperature range in their study. In addition, we take the molar


Page 11 of 47

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


volume 36 at 293 K to be 1.376 × 10−4 m3 /mole, corresponding to a density of 1650 kg/m3 . We solve the differential equation (10) with these assumptions on v and α to estimate the molar volume of solid TNT (in units of m3 /mole) at a temperature T to be h

i

v(T, P = 1 bar) = (0.1376 × 10−4 ) exp (1.31 × 10−6 )(T − 293)2 + (2.06 × 10−4 )(T − 293) .

Applying the steps detailed above, we are able to represent ai (T ) as

ai (T ) = a1,i T + a0,i ,

(11)

ai (T ) = (a1,i T + a0,i )2 ,

(12)

for solid RDX and solid TNT, and

for liquid TNT. Table 1 presents the coefficients used in these equations, as well as the values for the constants bi . Once the critical temperature and pressure of TNT have been accurately determined, one can adjust bi of liquid TNT and employ cubic splines to extend (12) up to the critical point so that it is compatible with (1). Table 1 also shows the average absolute deviation (AAD) and maximum deviation between the calculated density and experimental data over the indicated temperature range. We can achieve fairly accurate density results by selecting appropriate values for bi ; neither of the two density corrections described near the beginning of this section (i.e., adding volume translation functions 25,26 or making the bi themselves temperature-dependent 23,27,28 ) are necessary. We have found that allowing for temperature-dependent bi can further improve the density results, but the other results shown in this study remain comparable; see the Supporting Information for more details on this alternative formulation. Since temperature-dependent bi offer only a marginal improvement at the expense of a more complicated model, we elect to use a constant value for all bi .



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 47

Table 1: Value of bi and coefficients for calculating ai (units of L2 ·bar/mole2 ) in (11) and (12). The table also shows the average absolute deviation (AAD) and maximum deviation between the calculated density and experimental data 39–41,44 over the indicated temperature range. Compound Solid RDX Solid TNT Liquid TNT

bi Deviation in ρ T range for ρ (L/mole) (AAD %/Max %) (K) -0.2754556 246.4571708 0.1187689 0.49/1.9 293–478 -0.3763838 254.3240616 0.1309095 0.14/0.3 293–353 -0.0120155 15.8402463 0.1459449 0.71/1.1 354–478 a1,i

a0,i

Melting behavior of pure TNT The melting temperature Tmelt of TNT as a function of pressure is shown in Figure 1. Along the melting curve, the fugacity coefficients of the liquid (l) and solid (s) satisfy φlTNT (T, P ) = φsTNT (T, P ). The normal melting temperature Tmelt,n that we obtain from the Peng-Robinson EOS is 353.65 K, which is comparable to the experimental values of 353.35 K, 353.5 K, and 354.05 K reported by Campbell and Kushnarov, 34 Dattelbaum et al., 38,45 and Parker and Thorpe, 35 respectively. Experimental data of the melting behavior at higher pressures has recently become available with the work of Dattelbaum et al. They have used optical/spectroscopic methods (Raman, far-IR, and x-ray diffraction) to determine the melting point of purified TNT crystals heated in a pressurized diamond anvil cell. They estimate the precision in their temperature measurements to be within 5 K, while the uncertainty in the pressure readings is ± 0.1 GPa (± 1000 bar). Their study examines the behavior up to pressures as high as 12 GPa (120,000 bar), including one data point at an elevated, but low enough pressure (1500 bar) where we expect the Peng-Robinson EOS to still be valid. Our results lie within the experimental uncertainties, and in fact, we argue that Tmelt at 1500 bar must be higher than the experimental value of 373.15 K indicated in Figure 1. Our reasoning involves consideration of the dashed line in the figure, whose slope is calculated from three quantities evaluated at atmospheric pressure: the normal melting temperature Tmelt,n , the volume change upon melting


Page 13 of 47

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


420

354.2

410 ! $

400

354

" #

"#

353.8 390 380

353.6

370

353.4

360

!

353.2

$

350 0

500

1000

1500

2000

0

2500

(a)

5

" #

"#

10

15

20

(b)

Figure 1: The melting temperature Tmelt of TNT vs. pressure. A magnified view of (a) near atmospheric pressure is shown in (b). Our results lie within the uncertainties of Dattelbaum et al., 38 who estimate the error in their temperature and pressure measurements to be ± 5 K and ± 0.1 GPa (± 1000 bar), respectively. The slope of the dashed line is calculated from three quantities evaluated at atmospheric pressure: the normal melting temperature 0 Tmelt,n , the volume change upon melting ∆vfus , and the molar enthalpy of fusion ∆h0fus . By analyzing the slope of this line, we conclude that Tmelt at 1500 bar must be higher than the indicated experimental value of 373.15 K (see the text for more details).

0 ∆vfus = v l (Tmelt,n , P = 1 bar) − v s (Tmelt,n , P = 1 bar),

and the molar enthalpy of fusion

∆h0fus = hl (Tmelt,n , P = 1 bar) − hs (Tmelt,n , P = 1 bar). We denote the last two quantities, if they are evaluated at any other point along the melting curve, as ∆vfus and ∆hfus , respectively. We calculate the molar enthalpy departure function for both the liquid and solid forms of pure component i according to

hi −

hig i

= −RT

2

∂ ln φi (T, P ) ∂T

!

.

(13)

P

This equation can be derived from (S13) and the fundamental relation µi (T, P ) = hi − T si ,



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 47

where si = −(∂µi /∂T )P is the molar entropy of i. 19,26 From (13), we have

∆hfus =

2 −RTmelt

"

∂ ln φl (Tmelt , P ) ∂T

!

P

∂ ln φs (Tmelt , P ) − ∂T

! #

.

(14)

P

The volume change upon melting ∆vfus can be obtained from the compressibility factors Z l and Z s , or equivalently, it can be calculated from a relation that is analogous to (14):

∆vfus = RTmelt

"

∂ ln φl (Tmelt , P ) ∂P

!

T

∂ ln φs (Tmelt , P ) − ∂P

! #

.

(15)

T

The equality between the chemical potential of the liquid and that of the solid along the melting curve requires that its derivative be given by Tmelt ∆vfus dT = . dP ∆hfus

(16)

For the pressures shown in Figure 1, ∆vfus and ∆hfus undergo only moderately small changes 0 (less than 10%) compared to their values ∆vfus and ∆h0fus at atmospheric pressure. Thus, the

dashed line virtually overlaps the melting curve at lower pressures, and the deviation becomes noticeable — though still small (difference is less than 2%) — only at higher pressures. The fact that the dashed line in Figure 1 approximates the behavior of the melting curve (the solid line) implies that the melting temperature of TNT at 1500 bar must be higher than 373.15 K. The slope of the dashed line would have to be significantly smaller than that depicted in the figure for one to reasonably expect the melting curve to pass through the point (P = 1500 bar, T = 373.15 K). Since our normal melting temperature Tmelt,n closely matches experimental values, we examine the possibility that we have considerably 0 underestimated (overestimated) ∆h0fus (∆vfus ). By evaluating (14) at atmospheric pressure,

we obtain ∆h0fus = 20.5 kJ/mole, which is close to the experimental value of 21.2 kJ/mole reported by Acree 46 and the ∆h0fus of 20.1 ± 0.4 kJ/mole predicted by Neyertz et al. 47 in their molecular dynamics simulations. Cady and Rogers 40 reference several other studies 14 ACS Paragon Plus Environment

Page 15 of 47

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


that have measured ∆h0fus ; the results they cite range from 20.2 cal/g (19.2 kJ/mole) to 25.2 cal/g (23.9 kJ/mole). We have adjusted the value of bi using experimental liquid TNT density data (the deviation is 0.9% at Tmelt,n ), so we expect most of the uncertainty concern0 ing ∆vfus to be associated with the density of the solid. Density vs. temperature information

is not readily available for solid TNT, and instead we have estimated its density from thermal coefficient of expansion data. 41,44 We calculate the density of solid TNT at Tmelt,n to be 1617.1 kg/m3 , which is similar to the density of 1632.6 kg/m3 predicted Neyertz et al. at the melting temperature. Thus, our results for Tmelt,n , ∆h0fus , v l (Tmelt,n , P = 1 bar), and v s (Tmelt,n , P = 1 bar) have all either been fitted to experimental data or closely match previously published results. Unless there are significant errors in these past studies, it is clear that the melting temperature at 1500 bar must be higher than 373.15 K. Volumetric and thermal properties We have described how ai and bi are calibrated with atmospheric-pressure volumetric data in the form of density vs. temperature correlations (solid RDX and liquid TNT) or coefficient of expansion results (solid TNT). We can investigate the volumetric behavior at higher pressures by examining the isothermal coefficient of compressibility β, which is defined as 1 β=− v

∂v ∂P

!

.

T

The reciprocal of β is the isothermal bulk modulus KT = 1/β. Our predictions for the variation of KT with temperature and pressure are presented in Figure 2. The bulk moduli of energetic materials have been studied both experimentally and computationally using one of two different approaches. One approach is to measure or simulate the components of the elastic tensor of a crystal, and then estimate KT by applying an averaging scheme, such as the Voigt or Ruess schemes, to the tensor components. Another approach is to measure (e.g.,



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 47

through diffraction techniques) or simulate the change in the crystal lattice parameters as the compound is pressurized, calculate the corresponding change in the unit cell volume, and fit the volume vs. pressure results to a semi-empirical relation, such as the Murnaghan or Vinet equations, which yields an estimation of KT at atmospheric pressure. (Technically, this second approach produces KT at zero pressure, but the difference between KT at P = 1 bar and that at P = 0 bar is negligible.) 14

30 12

25 10 8

!" # $%%& " ' ( '' $%%% ) ' $%% ( " ' $% *

20 15

6 4

10

2

5 0

0

0

300 325 350 375 400 425 450 475

(a)

500

1000

1500

2000

2500

(b)

Figure 2: Peng-Robinson EOS predictions for the variation of the isothermal bulk modulus KT with (a) temperature at atmospheric pressure and (b) with pressure at 293.15 K. Representative values from the literature are shown for comparison. 48–53 Most of the past studies report the value of the bulk modulus at atmospheric pressure and room temperature, which we take to be 293.15 K, unless otherwise stated (i.e., some authors specify T to be 298 K instead). Our predictions for RDX differ from the results of Yoo et al. 48 and Sorescu et al. 49 by 0.7% and 4.2%, respectively. Sewell and Bennett 50 have calculated the Young’s modulus and Poisson ratio at 304 K and 333 K through Monte Carlo simulations. Their values for KT that we infer from these quantities are 8.1% and 6.3% different from our results at 304 K and 333 K, respectively. The KT obtained by fitting the results of Bowden et al. 53 to the Vinet equation and averaging over the three trials they have carried out yields a TNT bulk modulus of 7.9 GPa. This is 9.9% larger than our predicted value. Figure 2 shows that our predictions are comparable to the results from past studies. The elastic tensor components of RDX at ambient conditions have been measured by various authors through different techniques (ultrasonic resonant measurements, resonant ultrasound spectroscopy, impulsive stimulated thermal scattering, and Brillouin spectroscopy). 51,54 Some 16 ACS Paragon Plus Environment

Page 17 of 47

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


studies also present results for high pressures relevant to shock formation in energetic materials (e.g., see the work of Sewell and Bennett 50 ), but these pressures are typically greater than 1 GPa, beyond the expected range of validity for the Peng-Robinson EOS. Haycraft et al. 55 compare results for KT from a number of experimental and computational studies. Although there is a sizeable disagreement among a few studies over the value of a particular elastic constant (about a 40% deviation 54–56 ), the reported results for KT all lie within the range 11.2–12.3 GPa, if one takes the true value of KT to be the average of that produced by the Voigt and Ruess schemes. Sewell and Bennett present Young’s modulus E and Poisson ratio ν predictions for RDX at 304 and 333 K. These quantities are obtained from Monte Carlo simulations of the elastic constants. If we average their E and ν results along the different directions, we can infer KT through the equation KT = E/[3(1 − 2ν)]. Their KT differ from our predictions at 304 K and 333 K by 8.1% and 6.3%, respectively. The second approach described above yields similar, though somewhat higher numbers for KT . Molecular simulation results 49,57 fitted to the Murnaghan equation have predicted KT of RDX to be 12.3–12.9 GPa, comparable to the experimental results of 13.0 and 13.9 GPa obtained by Olinger et al. 58 and Yoo et al., 48 respectively. This last value closely agrees with our prediction of 14.0 GPa. For TNT, Neyertz et al. 47 report from their molecular simulations performed at 298 K that KT = 10.4 GPa. Stevens et al. 52 have correlated their volume vs. pressure experimental data to the Murnaghan equation, yielding a KT of 8.52 GPa. Bowden et al. 53 point out, however, that the particular choice of the semiempirical relation can lead to noticeably different estimations of the bulk modulus. For three different experimental trials, their results for KT at 298 K vary between 9.6 ± 0.9 GPa and 10.8 ± 0.6 GPa if fitted to the Murnaghan equation, while they vary between 7.2 ± 1.6 GPa and 8.8 ± 0.7 GPa if fitted to the Vinet equation. The average KT for the three trials are 10.1 GPa and 7.9 GPa upon applying the Murnaghan and Vinet equations, respectively. The latter is about 10% different from our prediction of 7.12 GPa at 298 K. 17 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

Page 18 of 47

Figure 3 illustrates that even if there are modest, but non-negligible disagreements in the bulk modulus at atmospheric pressure, they ultimately translate to small deviations in the density at the pressures for which the Peng-Robinson EOS is applicable. Our predicted value of KT differs from those of Sorescu et al. and Stevens et al. by 4.2% and 16.4%, respectively. Yet in both cases, the deviation in the density is less than 1% even at P = 2500 bar. 1850

1

1800

0.8

1750

0.6

1700

!

0.4

"##$

1650

0.2

1600 0

500

1000

1500

2000

0

2500

0

500

1000

1500

2000

2500

(b)

(a)

Figure 3: Comparison of the Peng-Robinson EOS predictions for the densities of solid RDX and TNT with the Murnaghan equation-fitted results of Sorescu et al. 49 and Stevens et al. 52 Our constant-pressure molar heat capacity cP predictions agree well with experimental data (Figure 4). It follows from (13) that cP = (∂h/∂T )P of a pure component is given by

cP =

cig P

− 2RT

∂ ln φ(T, P ) ∂T

!

P

− RT

2

∂ 2 ln φ(T, P ) ∂T 2

!

.

(17)

P

Two studies have performed quantum simulations to compile ideal gas molar heat capacity cig P predictions for many energetic materials. 59,60 We use the results from Osmont et al. because their work involves fewer simplifications and more sophisticated basis sets. We have fit their data in the temperature range 300–800 K to cubic polynomials so that

3 2 cig P = c3 T + c2 T + c1 T + c0 .


(18)

Page 19 of 47

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 coefficients are listed in Table S5. Nearly all the experimental studies have measured heat capacities directly through differential scanning calorimetry. 61–67 One study has inferred heat capacities from thermal conductivity and diffusivity measurements. 68 Baytos has examined many different energetic materials, and for each material he has obtained a linear fit of his data to temperature. 63 The average absolute deviations (AADs) between our predictions and his linear correlations are 6.7%, 7.0%, and 9.8% for solid RDX, solid TNT, and liquid TNT, respectively. Compared to the measurements of Ellison et al. (Figure 4 presents their results for flake TNT), the AADs are 1.1% and 1.4% for solid TNT and liquid TNT, respectively. 66 42

50 45

38

! # & (

" $

%

40

'

34 35 30

30

25

26 20

300

325

(a)

350

375

400

425

450

300

475

325

350

375

400

425

450

475

(b)

Figure 4: Constant-pressure molar heat capacity cP of (a) RDX and (b) TNT vs. temperature at atmospheric pressure. Representative experimental values are included for comparison. 61–68 The results are non-dimensionalized by R. The average absolute deviations (AADs) between our predictions and Baytos’ linear fits to temperature 63 are 6.7%, 7.0%, and 9.8% for solid RDX, solid TNT, and liquid TNT, respectively. The AADs compared to the results of Ellison et al. 66 are 1.1% and 1.4% for solid TNT and liquid TNT, respectively.

Liquid RDX As we have discussed, the critical properties of energetic materials have not been measured (and perhaps cannot be measured) experimentally. The few theoretical results that are available in the literature are therefore subject to a great deal of uncertainty. One possible exception may be the critical temperature Tc,RDX of RDX, which has been determined 19 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

Page 20 of 47

theoretically by a few authors 29,30 to be about 840 K. This is also the accepted value in studies that use the critical temperature to predict other thermophysical properties of RDX. 69 Using this value for Tc,RDX , we have initially attempted to represent aRDX and bRDX of liquid RDX according to the functional forms in (4) and (5), respectively. As done in many past vapor-liquid equilibria studies on the Peng-Robinson EOS, we set ψb (Tr,RDX ) = 1 for all temperatures so that bRDX is a constant, and have tried to calculate ψa (Tr,RDX ) from a correlation that involves the acentric factor and some adjustable parameters. Unlike the case for the other compounds we have considered, the variation of both vapor pressure and density with respect to temperature are not known for liquid RDX. Instead, we can calibrate aRDX and bRDX to fit three quantities that are known experimentally: Tmelt,n , 0 ∆vfus , and ∆h0fus . Just like for TNT, we calculate Tmelt,n of RDX by finding the temperature

at which the fugacity coefficients of the liquid and solid phases are equal to each other at 0 P = 1 bar. The molar enthalpy of fusion ∆h0fus and volume change upon melting ∆vfus are

calculated by evaluating (14) and (15), respectively, at P = 1 bar. In addition, the heat capacity cp of liquid RDX must be greater than that of solid RDX. We have found that we cannot satisfy all of the required conditions by representing aRDX and bRDX as described in the previous paragraph. We have tried many different functional forms for ψa (Tr,RDX ), and irrespective of our choice, if we adjust the parameters to match experimental values for 0 Tmelt,n , ∆vfus , and ∆h0fus , our predicted liquid RDX heat capacity is 20–30% less than that

of the solid at Tmelt,n . Although experimental data on the liquid-phase heat capacity is not available, it is very likely to be larger than cp of solid RDX. We can obtain more physically realistic liquid RDX heat capacity predictions by combining melting point data with our solid RDX results to infer the vapor pressure of the liquid. Vapor-liquid equilbria can be approximately described by the Clausius-Clapeyron equation P2 ∆hvap ln = P1 R

1 1 − , T1 T2


(19)

Page 21 of 47

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


Table 2: Comparison of our values for the RDX normal melting point Tmelt,n and molar enthalpy of fusion ∆h0fus with experimental data. As explained in the text, our results are obtained by using the experimental data to infer the liquid RDX vapor pressure and represent it according to the Clausius-Clapeyron equation (19). Nauflett et al. do not report 0 ∆h0fus results, but from their high-pressure optical setup, they have found that ∆vfus /∆h0fus is approximately 4 × 10−4 cm3 /cal = 9.56 × 10−8 m3 /kJ. We have adjusted the parameter bRDX 0 so that our molar volume change upon melting ∆vfus satisfies this ratio. Study This work Kishore 70 Zeman 71 Hall 72 Nauflett et al. 73 Nauflett et al. 73

Method Tmelt,n (K) ∆h0fus (kJ/mole) Peng-Robinson EOS 478.15 33.90 Differential scanning calorimetry 477.4 30.71 ± 0.26 Differential scanning calorimetry 477.5 ± 0.2 32.90 ± 0.73 Differential scanning calorimetry 478.5 35.65 ± 0.29 Differential scanning calorimetry 480.2 – High-pressure optical setup 484.2 –

where ∆hvap is the molar enthalpy of vaporization. Equation (19) is valid for conditions where ∆hvap can be treated as a constant and the vapor phase can be modeled as an ideal gas. We verify both assumptions later. Table 2 shows experimental values for Tmelt,n and ∆h0fus . Since our solid RDX model predicts the molar enthalpy of sublimation at Tmelt,n to be roughly 124 kJ/mole, we set ∆hvap at this temperature to be 90 kJ/mole so that ∆h0fus lies in the range of experimental values. The vapor pressure curves of the liquid and solid intersect at the triple point temperature, which is very close to Tmelt,n . We take the triple point temperature to be 478.15 K, and the vapor pressure of solid RDX at this temperature is 9.84 × 10−4 bar. We designate these two quantities as T1 and P1 , respectively, and fit the parameters aRDX and bRDX to the vapor pressure curve in (19) and to volumetric data. 0 Liquid RDX volumetric data regarding ∆vfus can be deduced from the measurements of

Nauflett et al., 73 who have studied the melting behavior with differential scanning calorimetry and with a high-pressure optical setup that is similar to the one used by Dattelbaum et al. 38 for TNT. Although the optical setup likely produces less accurate results than does differential scanning calorimetry, as seen by the fact that the Tmelt,n obtained from this method is a little higher than the other values in Table 2, it is still valuable because it allows for melt-



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 22 of 47

ing temperature measurements at elevated pressures. The derivative of the melting curve is given by (16), and for relatively low pressures this curve can be approximated as a straight 0 line with slope dT /dP ≈ Tmelt,n ∆vfus /∆h0fus . From examining the slope of their melting 0 curve, Nauflett et al. have determined ∆vfus /∆h0fus to be approximately 4 × 10−4 cm3 /cal 0 = 9.56 × 10−8 m3 /kJ. We choose bRDX = 0.1216931 L/mole so that our ∆vfus /∆h0fus ratio

agrees with this result. In addition, aRDX in units of L2 ·bar/mole2 is given by

aRDX (T ) = exp(a1,RDX T + a0,RDX ),

(20)

where a1,RDX = −1.1373415 × 10−3 K−1 , and a0,RDX = 5.3091496. In addition to matching Tmelt,n (we note that our value is higher than the triple point 0 temperature by only 0.004 K), ∆vfus , and ∆h0fus data, the calculated heat capacity cp of

liquid RDX is also greater than that of the solid. We predict cp to be 39.6R and 41.1R at Tmelt,n for solid RDX and liquid RDX, respectively. We have also verified that cp of both phases increases monotonically with temperature. Our results are based on (19), which is predicated on the assumptions that ∆hvap is a constant and the vapor phase is an ideal gas. In order for these requirements hold, we have fit aRDX to (19) in only the narrow range 478.15– 480.15 K. Over these temperatures, we find that ∆hvap decreases slightly from 90.0 kJ/mole to 89.9 kJ/mole, while Z = P v/RT of the vapor phase is approximately unity (0.99992). It is therefore justified to represent the vapor pressure of liquid RDX with (19). One can interpolate (20) to the critical point once it has been accurately determined. We have seen that there is much uncertainty regarding the critical pressure. Although there is agreement on the critical temperature Tc,RDX , one may also question the methods by which it has been obtained. For example, Maksimov 29 calculates Tc,RDX from the following function of the normal boiling point Tb,RDX : Tc,RDX = 1.027Tb,RDX + 159. He estimates Tb,RDX to be about 664 K by extrapolating (19) to P2 = 1 bar. It is questionable whether the stated relation


Page 23 of 47

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


between Tc,RDX and Tb,RDX is valid, and also whether ∆hvap is a constant between 478–664 K. Similar questions can be raised for the group contribution methods of Toghiani et al. 30

Mixtures Phase diagrams and the binary interaction coefficients kij Experimental studies have found that RDX and TNT form a binary eutectic system, as presented in the phase diagrams of Figure 5. Campbell and Kushnarov 34 have determined the eutectic point at atmospheric pressure to be at an RDX mole fraction of 0.0511 and lying in the temperature range 351.25–351.65 K, while Parker and Thorpe 35 report it to be at 0.0347 mole fraction RDX and 352.75 K. By definition, the lowest temperature at which liquid may exist is the eutectic temperature. In the RDX/TNT system, three phases can coexist in equilibrium at this temperature: liquid RDX and TNT, RDX-rich solid, and TNTrich solid. The liquid phase is represented by the eutectic point. At any other temperature, an RDX/TNT mixture exists either as a single-phase or two-phase mixture, depending on the overall composition. Figure 5 also shows experimental results for the liquidus curve. The right branch of this curve represents the composition of the liquid phase that it is in equilibrium with RDX-rich solid, while the left branch (seen more clearly in Figure 5b) represents the composition of the liquid phase that it is in equilibrium with TNT-rich solid. To the best of our knowledge, the only study that has theoretically modeled the phase behavior of the RDX/TNT binary system is the recent work by Jaansalu. 14 He has followed the excess Gibbs energy approach that we have briefly described in the Introduction. This approach entails calculating the chemical potential of the liquid and solid phases of the pure components (with respect to an arbitrary reference state) by integrating the constantpressure molar heat capacity cP over a temperature interval. Jaansalu represents the excess Gibbs energy with a Margules model, 11,14 from which he obtains the activity coefficients that 23 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

Page 24 of 47

355

475

354

450 )

425 400

*

353 !" $ !%&'

#

352

!%'

375

351

(

350

350

!" $ !%&'

#

325

!%'

349

*

300

( (

348

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

(a)

0

1

0.01

0.02

!"

0.03

0.04

0.05

0.06

(b)

Figure 5: Phase diagram of the RDX/TNT binary eutectic system at atmospheric pressure. A magnified view of (a) near the eutectic point is shown in (b). The dashed tie lines are imaginary horizontal lines that connect the three phases that are in equilibrium at the eutectic temperature (liquid RDX and TNT, RDX-rich solid, and TNT-rich solid). The 2 two-phase regions (liquid-solid and solid-solid) are labeled in (a). The liquid-solid region is enclosed by the right (left) branch of the liquidus curve and the curve representing the composition of RDX-rich (TNT-rich) solid in equilibrium with the liquid. The solid-solid region is enclosed by the two solid-phase curves. For clarity, in (a) we include Jaansalu’s theoretical predictions for only the liquidus curve. We have reproduced his results, with permission, 74 by digitizing the phase diagrams presented in his recent paper. 14 The average absolute deviation and root-mean-square deviation (RMSD) between our results and Campbell and Kushnarov’s experimental results 34 for the right branch of the liquidus curve are 3.69% and 0.0119, respectively. relate the chemical potentials in the mixture to the pure-component chemical potentials. For all p phases of a multiphase, multicomponent mixture to be in equilibrium at a given P and T , the equality of chemical potentials across the phases requires that fugacity fij of component i in phase j satisfy 10,11,24

fi1 (T, P, z1 ) = fi2 (T, P, z2 ) = · · · = fij (T, P, zj ) = · · · = fip (T, P, zp ),

∀i.

(21)

Equation (21) is the fundamental condition for thermodynamic equilibrium that we solve to obtain all phase boundaries, including the eutectic point. We obtain the fugacity fi from the fugacity coefficient φi = fi /zi P , which for the Peng-Robinson EOS is given by (S12).


Page 25 of 47

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 mixture parameter A = aP/R2 T 2 in (S12) is calculated from the pure-component parameters ai according to (S2) and (S3). These mixing rules involve the binary interaction l coefficients kij = kij (T ). In this study, we use two binary interaction coefficients: kRDX-TNT for s liquid RDX-liquid TNT and kRDX-TNT for solid RDX-solid TNT. The three-phase equilibrium

at the eutectic temperature is described by (21), which denotes a set of four equations at this temperature — two for each component. We denote the eutectic temperature at P = 1 bar as Teutectic,n . Since the composition of the liquid phase is known experimentally at Teutectic,n , we solve (21) to find the four remaining unknowns. These unknowns are the mole fraction l of RDX in the two solid (i.e., RDX-rich and TNT-rich) phases and the values of kRDX-TNT s and kRDX-TNT at Teutectic,n . Mixture properties at Teutectic,n are presented in Table 3. The value

of Teutectic,n itself and the liquid composition are obtained from Campbell and Kushnarov. We prefer using the data from their study over Parker and Thorpe’s because the former has examined a much wider range of conditions, as evident in Figure 5a. Table 3: Our results for the eutectic temperature Teutectic,n at P = 1 bar. The compositions of the three phases listed here are expressed in terms of the RDX mole fraction. The value of Teutectic,n itself and the liquid composition are obtained from Campbell and Kushnarov, 34 while the other four quantities are calculated by solving the conditions for equilibrium in (21). Teutectic,n (K) 351.45

l s kRDX-TNT kRDX-TNT Liquid RDX-rich solid TNT-rich solid -0.0168302 0.0931277 0.05107 0.99481 0.00919

The average absolute deviation and root-mean-square deviation (RMSD) between our results and Campbell and Kushnarov’s experimental results for the right branch of the s liquidus curve are 3.69% and 0.0119, respectively. We treat kRDX-TNT as a constant, but l allow kRDX-TNT (T ) to depend on temperature so that we obtain good agreement with their

solubility measurements. This dimensionless binary interaction coefficient is given by

l (T ) = k3 (T − Teutectic,n )3 + k2 (T − Teutectic,n )2 + k1 (T − Teutectic,n ) + k0 , kRDX-TNT


(22)


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 26 of 47

where k0 is listed in Table 3, k1 = 5.185 × 10−4 , k2 = −5.925 × 10−6 , and k3 = −5.721 × 10−9 . Jaansalu’s phase behavior results and ours are qualitatively the same, but different quantitatively. Part of the reason is because he has adjusted the parameters in his excess Gibbs energy model to match the eutectic point reported by Parker and Thorpe, rather than Campbell and Kushnarov. As a result, his eutectic point is characterized by a slightly higher temperature (352.75 K vs. 351.45 K) and a lower RDX mole fraction (0.0347 vs. 0.0511). The lower RDX solubility at the eutectic point is reflected in lower RDX solubilities in all of his phase boundary curves. At the resolution of Figure 5a, the differences are not easily discernible for the solid-phase boundaries, but they are noticeable for the liquidus curve, especially at temperatures between 400–450 K. It is clear that more data is needed, especially for the solid-phase curves, to verify the limited experimental results that are presently available and to validate the theoretical predictions. We have asserted in the Introduction that the Peng-Robinson EOS is capable of making phase behavior predictions at higher pressures without requiring any additional data. For a given pressure P , we first solve the equilibrium conditions in (21) to find the corresponding eutectic temperature and the RDX mole fraction of the three phases that coexist in equilibrium at this temperature and pressure. Equation (21) also provides the equations from which the phase boundaries at all other temperatures are determined. Note that the binary interaction coefficients depend on temperature only, so they are already known from (22) and Table 3. Figures 6 and 7 illustrate the phase behavior at pressures up to 2500 bar. The eutectic point appears at a higher temperature and RDX solubility as the pressure is increased. Due to the increase in the eutectic temperature, the two solid phases that are in equilibrium with the liquid phase at this temperature become less pure. That is, as the pressure is increased the RDX-rich (TNT-rich) solid has a higher solubility of TNT (RDX) at the corresponding eutectic temperature. The left and right branches of the liquidus curve terminate at higher temperatures as well. The leftmost and rightmost points of this curve represent the 26 ACS Paragon Plus Environment

Page 27 of 47

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


475

475

450

450

425

425

400

400

375

375

350

350

325

325

300

300

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

(a)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

(b)

475

475

450

450

425

425

400

400

375

375

350

350

325

325 300

300 0

(c)

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

(d)

Figure 6: Peng-Robinson EOS predictions for the phase diagram of the RDX/TNT binary system at pressures of (a) 500 bar, (b) 1000 bar, (c) 2000 bar, and (d) 2500 bar.



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

374

386

372

384

370

382

368

380

366

378

364

376

362

374

0

0.04

0.08

0.12

0.16

0.2

(a)

0

0.04

0.08

0.12

0.16

0.2

0

0.04

0.08

0.12

0.16

0.2

(b) 410

420

408

418

406

416

404

414

402

412

400

410

398

408

0

(c)

Page 28 of 47

0.04

0.08

0.12

0.16

0.2

(d)

Figure 7: Peng-Robinson EOS predictions for the phase diagram of the RDX/TNT binary system near the eutectic point at pressures of (a) 500 bar, (b) 1000 bar, (c) 2000 bar, and (d) 2500 bar.


Page 29 of 47

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


melting points of pure TNT and RDX, respectively. One can verify, for example, that the TNT melting points indicated in Figure 7 are the same as those in Figure 1. The liquidus curve represents the freezing points of the RDX/TNT system. This is similar to how the liquid side of a vapor-liquid envelope represents the locus of bubble points. The average slope of the liquidus curve between 1 mole percent RDX and pure TNT (0 mole percent RDX) is -0.42 K/mole percent RDX at atmospheric pressure, and -0.45 K/mole percent RDX at 2500 bar. The increase in magnitude of the slope implies that the effect of RDX on the freezing point depression of TNT is more pronounced at the higher pressure. By doing the same comparison for the slope of the right branch of the liquidus curve, we find the effect of TNT on the freezing point depression of RDX is less pronounced at 2500 bar than at 1 bar.

Volumetric and thermal properties In this section, we examine the volumetric and thermal properties of Composition B-3 (Comp B-3), which has an overall composition of 60% RDX and 40% TNT, where the percentages indicate mass fractions. Comp B-3 is a relatively well-studied binary RDX/TNT mixture that is similar in its overall composition to Composition B (63% RDX, 36% TNT, and 1% wax by mass), the energetic material that is the central focus of some of our research group’s recent work. 6,7 We can find the phase(s) that make up Comp B-3 at any temperature and pressure by applying the lever rule to a vertical line that intersects the abscissa of the appropriate phase diagram at 0.605 (the overall RDX mole fraction of Comp B-3). This procedure reveals the composition of the phase(s) and their relative amounts, as shown in Figure 8. This figure portrays how the phase mole fraction xj , which represents the mole fraction of phase j in the overall Comp B-3 mixture, varies with temperature at atmospheric pressure. The composition of these phases, which we present in Figure 8b, are obtained from the phase boundaries in Figure 5. At room temperature, Comp B-3 is composed of the two solid phases, both of which are nearly pure. The solids become slightly less pure in 29 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

Page 30 of 47

composition as they are heated to the eutectic temperature Teutectic,n of 351.45 K. The TNTrich solid disappears and a new phase — the liquid — forms at this temperature. The liquid phase has a higher solubility of RDX than does the TNT-rich solid, which is why the phase mole fraction of the RDX-rich solid decreases in a discontinuous manner at Teutectic,n . As the temperature is raised further, the RDX solubility of both phases increases, resulting in a concomitant increase (decrease) in the phase mole fraction of the liquid (RDX-rich solid). Comp B-3 becomes a single-phase liquid at a temperature of about 451.3 K, which we can label as its normal melting point. The composition of this phase is identical to the overall composition of Comp B-3. 1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 300

325

350

375

400

425

450

475

300

325

350

375

400

425

450

475

(b)

(a)

Figure 8: (a) Mole fraction xj of the phase(s) j that constitute Comp B-3 at atmospheric pressure, and (b) RDX mole fraction of the phases in (a). Comp B-3 exists as a solid-solid mixture from 293 K to the eutectic temperature Teutectic,n of 351.45 K, a liquid-solid mixture from Teutectic,n to 451.3 K, and a single-phase liquid above 451.3 K (normal melting point). The molar volume v of a multiphase, multicomponent mixture of p phases and c components at a pressure P and temperature T is

v = RT

p X

j=1

 c X 1 xj zj  i

i=1

∂ ln φji + P ∂P

!

T,zj



,

(23)

where zij and φji are the mole fraction and fugacity coefficient of component i in phase j, 30 ACS Paragon Plus Environment

Page 31 of 47

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


respectively. 19 The quantity in brackets, when multiplied by RT , is the partial molar volume of component i in phase j. The molar density of the mixture is ρ = 1/v. Figure 9a compares the predicted density of Comp B-3 with that of pure RDX and TNT. The sharp drop in the density of TNT occurs at its normal melting point (353.65 K), which is slightly higher than Teutectic,n (351.45 K). There is a discontinuous decrease in the density of Comp B-3 at the eutectic temperature due to the simultaneous disappearance of the TNT-rich solid and formation of the liquid. The figure also includes calculations for the density of Comp B-3 if we assume that its constituent phase(s) form ideal mixtures. The difference between the Comp B-3 solid line and the dotted line therefore gives an indication of the excess volume of the phases. One defining feature of ideal mixtures is that the partial molar volumes are equivalent to the molar volumes of the pure components (i.e., that the volume change of mixing is zero). 19,24,26 As a result, the two solid phases, both of which are nearly pure in composition, have negligible excess volumes. The excess volume of Comp B-3 becomes noticeable (though still small) at temperatures above the melting point of 451.3 K. Figure 9a presents experimental data at room temperature from two sources. Baytos 63 and Dobratz 61 report the density of Comp B-3 to be 1725 kg/m3 and 1720 kg/m3 , respectively. The former is less than our predicted density of 1740.2 kg/m3 by 0.88%, while the latter is less than ours by 1.17%. We have not found density measurements at higher temperatures or pressures. Equation (23) can be combined with Figure 5 to find the volume fraction of the phase(s) that constitute Comp B-3. Figure 9b is the volumetric analogue of Figure 8a. We have included the phase mole fraction of the RDX-rich solid for comparison in this figure. The volume fraction of this phase is less than its mole fraction mainly because RDX is more dense than TNT. The volume fraction of the RDX-rich solid is particularly important when Comp B-3 is partially molten (exists as a liquid-solid mixture) because it can strongly influence the viscosity. Recently developed constitutive equations suggest that an increase of the solid volume fraction by only 10% can increase the Comp B-3 viscosity by severalfold. 6,75 31 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

1800

1

1750

0.9

1700

0.8

1650

0.7

1600

0.6 0.5

1550 1500 1450

#

0.4

!"! $ !"%

0.3

1400

0.2

1350

0.1 0

1300

300

300 325 350 375 400 425 450 475

(a)

Page 32 of 47

325

350

375

400

425

450

475

(b)

Figure 9: Peng-Robinson EOS predictions for the volumetric behavior of Comp B-3 and related substances at atmospheric pressure. The density of pure RDX, pure TNT, and Comp B-3 are compared in (a). The dotted line indicates our prediction for the density of Comp B-3 if its constituent phase(s) were to form ideal mixtures. The excess volume of Comp B-3, which is related to the difference between the Comp B-3 solid line and the dotted line, becomes noticeable at temperatures higher than the normal melting point of 451.3 K. The ambient-condition, experimental results of Baytos 63 and Dobratz 61 are less than our predicted density by 0.88% and 1.17%, respectively. The volume fraction of the phase(s) present in Comp B-3 are shown in (b). The mole fraction of the RDX-rich solid from Figure 8a is included for comparison. The volume fraction of this phase when Comp B-3 exists as a liquid-solid mixture (i.e., when the temperature is between 351.45 K and 451.3 K) is important because it can strongly affect the Comp B-3 viscosity. 6,75


Page 33 of 47

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


Our Comp B-3 heat capacity results, along with experimental data from Baytos, 63 are shown in Figure 10. We calculate the constant-pressure molar heat capacity cP from the following equation, which is essentially a multiphase, multicomponent extension of (17):

cP =

p X

j=1

 c X z j −2RT xj i

i=1

∂ ln φji ∂T

!

− RT 2

P,zj

∂ 2 ln φji ∂T 2

!

P,zj



 + cig P,i ,

(24)

ig where cig P,i is the ideal gas molar heat capacity of component i. The cP,i are given by (18)

and Table S5. Equation (24) is a thermal analogue of (23). Both can be derived in a similar manner from fundamental thermodynamic relations. 19 The quantity in brackets is the partial molar heat capacity of component i in phase j. Unlike the case for the excess volume, the excess heat capacity is noticeable well below the melting point. Baytos has fit his data for Comp B-3 to linear functions of temperature, just like he has done for pure RDX and TNT. These are also the same results published later in two widely-used compendia on energetic materials. 36,76 The AAD between Baytos’ linear correlation at temperatures below Teutectic,n and our predictions is 8.8%. The AAD is 10.3% at temperatures above Teutectic,n . 60 55 50 45 40 35 30 25 300

325

350

375

400

425

450

475

Figure 10: Comparison of Peng-Robinson EOS predictions for the Comp B-3 constantpressure molar heat capacity cP at atmospheric pressure with experimental data from Baytos. 63 He has fit his data to linear functions of temperature in two regimes: one below the eutectic temperature Teutectic,n and the other above Teutectic,n . The average absolute deviations between his linear fits in these two regimes and our results are 8.8% and 10.3%, respectively.



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 34 of 47

Conclusions We have shown that the Peng-Robinson equation of state can model pure RDX and TNT in their condensed phases, as well liquid and solid mixtures of the two compounds. Our work has analyzed temperatures and pressures as high as about 500 K and 2500 bar, respectively. These are conditions relevant to cookoff of mixtures containing RDX and TNT, such as Composition B, before explosion occurs. The Peng-Robinson equation of state in its most elementary form is given by the Helmholtz energy departure function in (S1), from which we derive all other expressions through fundamental thermodynamic relations. We have demonstrated that the equation of state provides a good representation of experimental volumetric (density and bulk modulus), thermal (heat capacity), and phase behavior (melting temperature and solubility) data. Predictions of derived properties — such as the TNT enthalpy of fusion, the isothermal bulk modulus (Figure 2), and the constant-pressure heat capacity (Figures 4 and 10) — lie within about 10% of the available experimental data, and many are within the error of the measurements. Heat capacity predictions can be improved further by allowing for temperature-dependent bi , as discussed in the Supporting Information. Although our results agree well with available data, more experimental measurements are needed to further validate the equation of state and to narrow the uncertainties in the existing data. We have seen that there is a fair amount of scatter in the experimental results of some properties. The uncertainties partly arise from differences in the samples used among the various studies. For example, single crystals can exhibit different properties than pressed materials. Impurities (e.g., contaminants) may have also been present in some samples. In addition, some of the scatter is undoubtedly due to the reactive nature of energetic materials. The decomposition of these substances is especially prominent at higher temperatures and pressures. Nonetheless, these difficulties must somehow be overcome to obtain the many vital pieces of information that are still currently missing. For example, the critical properties


Page 35 of 47

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


of energetic materials are almost absent in the literature, having been limited to a few theoretical studies 29,30,69 only and no experimental measurements. Vapor pressures 31,32 have been measured only to temperatures as high as about 480 K, well below the expected values of the critical temperatures (800 K and higher). Vapor pressure data, as well as density and heat capacity data, are still unavailable for liquid RDX. This is also likely true for several other compounds that exist only at high temperatures. More melting temperature measurements at pressures relevant to the Peng-Robinson equation of state are needed. Our melting point results for TNT (Figure 1) lie within the experimental uncertainties of Dattelbaum et al., 38 but we can only compare with two data points from their study, which is oriented towards pressures that are one or two orders of magnitude greater than those in this study. Moreover, we again note that experimental results for RDX/TNT mixtures are still quite limited. For example, no data exists for the solid-phase curves in the atmospheric-pressure phase diagram of Figure 5. It would also be important to have experimental validation of the phase behavior predictions at higher pressures (Figures 6 and 7). Density measurements of RDX/TNT mixtures like Composition B-3 are currently limited to ambient conditions (Figure 9). The experimental studies discussed in the previous paragraph are necessary to attain a more complete understanding of energetic materials, but the reality of the situation is that it may take decades — if not longer — for them to be completed. Given the data that is currently available, we conclude that the Peng-Robinson equation of state represents a good choice for modeling pure RDX and TNT, in addition to mixtures of these two compounds. One possibility to make the equation of state even more promising is to extend its applicability to higher pressure regimes. One way to do this may be to fit the pure-component parameters ai and bi to unreacted Hugoniot data. 77 The Hugoniot curve describes the relation between the pressure and density discontinuities that exist across shock waves generated in energetic materials. The pressure discontinuities typically lie somewhere between a few kilobars (say, 0.5 GPa) and 20 GPa. By adjusting ai and bi to fit measurements taken along 35 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

Page 36 of 47

both the Hugoniot and vapor pressure curves, it may be possible to extend the pressure range over which the Peng-Robinson equation of state is applicable. Finally, we end by pointing out that the equation of state can be used to model a multicomponent mixture of c components if the procedure described for the binary RDX/TNT system is applied to each of the c(c − 1)/2 binary pairs in the multicomponent system. The Peng-Robinson equation of state therefore provides the theoretical framework to model the thermodynamic behavior of energetic materials ranging in complexity from pure components to multiphase, multicomponent mixtures.

Acknowledgement The authors thank Drs. Craig Tarver and H. Keo Springer of Lawrence Livermore National Laboratory (LLNL) for helpful suggestions on an earlier draft of the manuscript. The authors also thank Drs. M. Riad Manaa, I-Feng (William) Kuo, and Sorin Bastea, all affiliated with LLNL, for their analysis and critique of the quantum calculations used in two studies 59,60 to determine ideal gas heat capacities. Professor Kevin Jaansalu of the Royal Military College of Canada and Dr. Dana Dattelbaum of Los Alamos National Laboratory graciously shared some of the data presented in their papers. 14,38 This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The research has been partially funded by the Joint DoD/DOE Munitions Technology Development Program (JMP).

Supporting Information Available Summary of important Peng-Robinson equation of state expressions referenced throughout the main text, as well as tables for Avap and Bvap (Table S4) and ideal gas heat capaci-


Page 37 of 47

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


ties cig P (Table S5). The Supporting Information also includes a comparison with results from an alternative Peng-Robinson model we have developed where bi of solid RDX, solid TNT, and liquid TNT vary with temperature.

This material is available free of charge via the

Internet at http://pubs.acs.org/.

References (1) Fried, L. E.; Manaa, M. R.; Pagoria, P. F.; Simpson, R. L. Design and synthesis of energetic materials. Annu. Rev. Mater. Res. 2001, 31, 291 – 321, DOI: 10.1146/annurev.matsci.31.1.291. (2) Politzer, P.; Boyd, S. Molecular dynamics simulations of energetic solids. Struct. Chem. 2002, 13, 105 – 113, DOI: 10.1023/A:1015748330357. (3) Fabbiani, F. P. A.; Pulham, C. R. High-pressure studies of pharmaceutical compounds and energetic materials. Chem. Soc. Rev. 2006, 35, 932 – 942, DOI: 10.1039/b517780b. (4) Yoh, J. J.; McClelland, M. A.; Maienschein, J. L.; Wardell, J. F. Towards a predictive thermal explosion model for energetic materials. J. Comput. Aid. Mater. Des. 2003, 10, 175 – 189, DOI: 10.1007/s10820-005-1750-z. (5) Tarver,

C.

M.;

Tran,

T.

D.

Thermal

decomposition

models

for

HMX-

based plastic bonded explosives. Combust. Flame 2004, 137, 50 – 62, DOI: 10.1016/j.combustflame.2004.01.002. (6) McClelland, M. A.; Glascoe, E. A.; Nichols III, A. L.; Schofield, S. P.; Springer, H. K. ALE3D simulation of incompressible flow, heat transfer, and chemical decomposition of Comp B in slow cookoff experiments. 15th International Detonation Symposium. 2014; 11 pp. https://e-reports-ext.llnl.gov/pdf/776766.pdf. 37 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

Page 38 of 47

(7) Nichols III, A. L.; Schofield, S. P. Modeling the response of fluid/melt explosives to slow cook-off. 15th International Detonation Symposium. 2014; 9 pp. https://e-reportsext.llnl.gov/pdf/776576.pdf. (8) Asante, D. O.; Kim, S.; Chae, J.; Kim, H.; Oh, M. CFD cook-off simulation and thermal decomposition of confined high energetic material. Propell. Explos. Pyrot. 2015, 40, 699 – 705, DOI: 10.1002/prep.201400296. (9) Cheetah thermochemical code, Lawrence Livermore National Laboratory, https:// www-pls.llnl.gov/?url=science_and_technology-chemistry-cheetah. (10) Elliott, J. R.; Lira, C. T. Introductory Chemical Engineering Thermodynamics, 1st ed.; Prentice Hall, Upper Saddle River, 1999. (11) Prausnitz, J. M.; Lichtenthaler, R. N.; Gomez de Azevedo, E. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall, Upper Saddle River, 1999. (12) McKenney Jr., R. L.; Krawietz, T. R. Binary phase diagram series: HMX/RDX. J. Energ. Mater. 2003, 21, 141 – 166, DOI: 10.1080/07370650390256656. (13) Sandler, S. I. Mixture phase equilibria involving solids; John Wiley & Sons, Hoboken, 2006; Chapter 12 of Chemical, Biochemical, and Engineering Thermodynamics, pp 658–702. (14) Jaansalu, K. M. Phase diagram modelling of cast energetic materials. IMEMTS 2015, Rome, Italy. 2015; 11 pp. http://www.imemg.org/wp-content/uploads/2015/06/8A217208-Phase-Diagram-Modelling-of-Cast-Energetic-Materials.pdf. (15) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. Fluid Phase Equilibr. 1989, 52, 31 – 38, DOI: 10.1016/0378-3812(89)80308-5. 38 ACS Paragon Plus Environment

Page 39 of 47

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


(16) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An equation of state for associating fluids. Ind. Eng. Chem. Res. 1996, 35, 4310–4318, DOI: 10.1021/ie9600203. (17) Li, Z.; Firoozabadi, A. Cubic-plus-association equation of state for water-containing mixtures: Is “cross association” necessary? AIChE J. 2009, 55, 1803–1813, DOI: 10.1002/aic.11784. (18) Kontogeorgis, G. M.; Folas, G. K. Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories, 1st ed.; John Wiley & Sons, Chichester, 2010; DOI: 10.1002/9780470747537. (19) Myint, P. C.; Hao, Y.; Firoozabadi, A. The CPA equation of state and an activity coefficient model for accurate molar enthalpy calculations of mixtures with carbon dioxide and water/brine; Technical Report LLNL-TR-672077, Lawrence Livermore National Laboratory, 2015; 43 pp., DOI: 10.2172/1194030, https://e-reports-ext.llnl.gov/pdf/ 790891.pdf. (20) Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundamen. 1976, 15, 59–64, DOI: 10.1021/i160057a011. (21) Wu, J.; Prausnitz, J. M. Phase equilibria for systems containing hydrocarbons, water, and salt: An extended Peng-Robinson equation of state. Ind. Eng. Chem. Res. 1998, 37, 1634 – 1643, DOI: 10.1021/ie9706370. (22) Boddu, V. M.; Maloney, S. W.; Toghiani, R. K.; Toghiani, H. In Energetic Materials: Thermophysical Properties, Predictions, and Experimental Measurements; Boddu, V. M., Redner, P., Eds.; CRC Press, Taylor & Francis Group, Boca Raton, 2010; Chapter 11: Solubility of RDX, HMX, and ǫ-CL20 in supercritical carbon dioxide, pp 199–220. 39 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

Page 40 of 47

(23) Xu, Z.; Sandler, S. I. Temperature-dependent parameters and the Peng-Robinson equation of state. Ind. Eng. Chem. Res. 1987, 26, 601–606, DOI: 10.1021/ie00063a030. (24) Sandler, S. I. Chemical, Biochemical, and Engineering Thermodynamics, 4th ed.; John Wiley & Sons, Hoboken, 2006. (25) Tsai, J.-C.; Chen, Y.-P. Application of a volume-translated Peng-Robinson equation of state on vapor-liquid equilibrium calculations. Fluid Phase Equilibr. 1998, 145, 193 – 215, DOI: 10.1016/S0378-3812(97)00342-7. (26) Firoozabadi, A. Thermodynamics and Applications in Hydrocarbons Energy Production; McGraw-Hill, New York, 2015. (27) Xu, Z.; Sandler, S. I. Application to mixtures of the Peng-Robinson equation of state with fluid-specific parameters. Ind. Eng. Chem. Res. 1987, 26, 1234–1238, DOI: 10.1021/ie00066a030. (28) Martynov, S.; Brown, S.; Mahgerefteh, H. An extended Peng-Robinson equation of state for carbon dioxide solid-vapor equilibrium. Greenhouse Gas. Sci. Technol. 2013, 3, 136–147, DOI: 10.1002/ghg.1322. (29) Maksimov, Y. Y. Boiling point and enthalpy of evaporation of liquid hexogen and octogen. Russ. J. Phys. Chem. 1992, 66, 280–281. (30) Toghiani, R. K.; Toghiani, H.; Maloney, S. W.; Boddu, V. M. In Energetic Materials: Thermophysical Properties, Predictions, and Experimental Measurements; Boddu, V. M., Redner, P., Eds.; CRC Press, Taylor & Francis Group, Boca Raton, 2010; Chapter 10: Prediction of physicochemical properties of energetic materials, pp 171–198.


Page 41 of 47

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


(31) Östmark, H.; Wallin, S.; Ang, H. G. Vapor pressure of explosives: A critical review. Propell. Explos. Pyrot. 2012, 37, 12–23, DOI: 10.1002/prep.201100083. (32) Ewing, R. G.; Waltman, M. J.; Atkins, D. A.; Grate, J. W.; Hotchkiss, P. J. The vapor pressures of explosives. TrAC, Trends in Analytical Chemistry 2013, 42, 35–48, DOI: 10.1016/j.trac.2012.09.010. (33) Edwards, G. The vapour pressure of 2:4:6-trinitrotoluene. Trans. Faraday Soc. 1950, 46, 423–427, DOI: 10.1039/TF9504600423. (34) Campbell, A. N.; Kushnarov, H. A. A study of systems involving trinitrotoluene, RDX, NENO, and dinitrobenzene, with particular reference to the composition of the binary eutectics. Can. J. Res. 1947, 25, 216–227, DOI: 10.1139/cjr47b-025. (35) Parker, R. P.; Thorpe, B. W. The phase diagram of the RDX/TNT system; Technical Note No. 140, Department of Supply, Australian Defence Scientific Service, Defence Standards Laboratories, 1970; 10 pp. (36) Dobratz, B. M.; Crawford, P. C. LLNL Explosives Handbook: Properties of Chemical Explosives and Explosive Simulants; Lawrence Livermore National Laboratory, 1985. (37) Miller, G. R.; Garroway, A. N. A review of the crystal structures of common explosives. Part 1: RDX, HMX, TNT, PETN, and Tetryl; Technical Report NRL/MR/6120– 01-8585, U.S. Naval Research Laboratory, 2001; http://www.dtic.mil/dtic/tr/ fulltext/u2/a396646.pdf. (38) Dattelbaum, D. M.; Chellappa, R. S.; Bowden, P. R.; Coe, J. D.; Margevicius, M. A. Chemical stability of molten 2,4,6-trinitrotoluene at high pressure. Appl. Phys. Lett. 2014, 104, 021911, DOI: 10.1063/1.4860395.



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 42 of 47

(39) Sun, J.; Shu, X.; Liu, Y.; Zhang, H.; Liu, X.; Jiang, Y.; Kang, B.; Xue, C.; Song, G. Investigation on the thermal expansion and theoretical density of 1,3,5trinitro-1,3,5-triazacyclohexane. Propell. Explos. Pyrot. 2011, 36, 341–346, DOI: 10.1002/prep.201000026. (40) Cady, H. H.; Rogers, W. H. Enthalpy, density, and coefficient of thermal expansion of TNT ; Technical Report LA-2696, Los Alamos National Laboratory, 1962; 16 pp., http://catalog.hathitrust.org/Record/012213646. (41) Vrcelj, R. M.; Sherwood, J. N.; Kennedy, A. R.; Gallagher, H. G.; Gelbrich, T. Polymorphism in 2-4-6 trinitrotoluene. Crystal Growth & Design 2003, 3, 1027–1032, DOI: 10.1021/cg0340704. (42) Heberlein, D. C. Thermal expansion of α-trinitrotoluene from 77 to 273 K. J. Chem. Phys. 1974, 61, 2346, DOI: 10.1063/1.1682313. (43) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Molecular packing and molecular dynamics study of the transferability of a generalized nitramine intermolecular potential to non-nitramine crystals. J. Phys. Chem. A. 1999, 103, 989–998, DOI: 10.1021/jp983847e. (44) Rosen, J. M.; Sickman, D. V.; Morris, W. W. The melting behavior of TNT ; Technical Report 6146, U.S. Naval Ordinance Laboratory, 1959; http://www.dtic.mil/dtic/ tr/fulltext/u2/309257.pdf. (45) Dattelbaum, D. M. Personal communication on August 27, 2015 to obtain a spreadsheet of the data used to create the phase diagram (Figure 3) published in Dattelbaum et al., Appl. Phys. Lett., 2014, 104, 021911. (46) Acree Jr., W. E. Thermodynamic properties of organic compounds: enthalpy of fusion 42 ACS Paragon Plus Environment

Page 43 of 47

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


and melting point temperature compilation. Thermochim. Acta 1991, 189, 37 – 56, DOI: 10.1016/0040-6031(91)87098-H. (47) Neyertz, S.;

Mathieu, D.;

Khanniche, S.;

Brown, D. An empirically opti-

mized classical force-field for molecular simulations of 2,4,6-trinitrotoluene (TNT) and 2,4-dinitrotoluene (DNT). J. Phys. Chem. A 2012, 116, 8374–8381, DOI: 10.1021/jp305362n. (48) Yoo, C.-S.; Cynn, H.; Howard, W. M.; Holmes, N. Equations of state of unreacted high explosives at high pressures. 11th International Detonation Symposium, Snowmass Village, CO, USA. 1998; 8 pp. https://e-reports-ext.llnl.gov/pdf/234078.pdf. (49) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Theoretical studies of the hydrostatic compression of RDX, HMX, HNIW, and PETN crystals. J. Phys. Chem. B 1999, 103, 6783–6790, DOI: 10.1021/jp991202o. (50) Sewell, T. D.; Bennett, C. M. Monte Carlo calculations of the elastic moduli and pressure-volume-temperature equation of state for hexahydro-1,3,5-trinitro-1,3,5triazine. J. Appl. Phys. 2000, 88, 88 – 95, DOI: 10.1063/1.373628. (51) Schwarz, R. B.; Hooks, D. E.; Dick, J. J.; Archuleta, J. I.; Martinez, A. R. Resonant ultrasound spectroscopy measurement of the elastic constants of cyclotrimethylene trinitramine. J. Appl. Phys. 2005, 98, 056106, DOI: 10.1063/1.2037865. (52) Stevens, L. L.; Velisavljevic, N.; Hooks, D. E.; Dattelbaum, D. M. The high-pressure phase behavior and compressibility of 2,4,6-trinitrotoluene. Appl. Phys. Lett. 2008, 93, 081912, DOI: 10.1063/1.2973162. (53) Bowden, P. R.; Chellappa, R. S.; Dattelbaum, D. M.; Manner, V. W.; Mack, N. H.; Liu, Z. The high-pressure phase stability of 2,4,6-trinitrotoluene (TNT). 18th 43 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

APS-SCCM and 24th AIRAPT, Seattle, WA, USA. 2014; DOI: 10.1088/17426596/500/5/052006, 6 pp. Published in J. Phys. Conf. Ser. 2014, 500, Part 5. (54) Bolme, C. A.; Ramos, K. J. The elastic tensor of single crystal RDX determined by Brillouin spectroscopy. J. Appl. Phys. 2014, 116, 183503, DOI: 10.1063/1.4901461. (55) Haycraft, J. J.; Stevens, L. L.; Eckhardt, J. J. The elastic constants and related properties of the energetic material cyclotrimethylene trinitramine (RDX) determined by Brillouin scattering. J. Chem. Phys. 2006, 124, 024712, DOI: 10.1063/1.2141958. (56) Johnson, J. A.; Manke, K. J.; Veysset, D. G.; Maznev, A. A.; Ramos, K. J.; Hooks, D. E.; Nelson, K. A. Photoacoustic determination of the speed of sound in single crystal cyclotrimethylene trinitramine at acoustic frequencies from 0.5 to 15 GHz. J. Appl. Phys. 2011, 110, 113513, DOI: 10.1063/1.3667291. (57) Agrawal, P. M.; Rice, B. M.; Zheng, L.; Thompson, D. L. Molecular dynamics simulations of hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using a combined SorescuRice-Thompson AMBER force field. J. Phys. Chem. B 2006, 110, 26185 – 26188, DOI: 10.1021/jp065241t. (58) Olinger, B.; Roof, B.; Cady, H. Symposium International Sur Le Comportement Des Milieux Denses Sous Hautes Pressions Dynamiques. Commissariat a l’Energie Atomique Centre d’Etudes de Vajours: Paris, France. 1978; p 3. (59) Burcat, A. Thermodynamic properties of ideal gas nitro and nitrate compounds. J. Phys. Chem. Ref. Data 1999, 28, 63–130, DOI: 10.1063/1.556033. (60) Osmont, A.;

Catoire, L.;

Gökalp, I.;

Yang, V. Ab initio quantum chemi-

cal predictions of enthalpies of formation, heat capacities, and entropies of gas-phase energetic compounds. Combust. Flame 2007, 151, 262 – 273, DOI: 10.1016/j.combustflame.2007.05.001. 44 ACS Paragon Plus Environment

Page 44 of 47

Page 45 of 47

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


(61) Dobratz, B. M. Properties of chemical explosives and explosive simulants; Technical Report UCRL-51319, Lawrence Livermore National Laboratory, 1972; 334 pp., DOI: 10.2172/4285272, http://www.osti.gov/scitech/biblio/4285272. (62) Velicky, R.; Lenchitz, C.; Beach, W. Enthalpy change, heat of fusion, and specific heat of basic explosives; Technical Report 2504, Picatinny Arsenal, Dover, New Jersey, 1959. (63) Baytos, J. F. Specific heat and thermal conductivity of explosives, mixtures, and plasticbonded explosives determined experimentally; Technical Report LA-8034-MS, Los Alamos National Laboratory, 1979; http://www.osti.gov/scitech/biblio/5913065. (64) Wilcox, J. D. Differential scanning calorimeter methods in the determination of thermal properties of explosives. M.Sc. thesis, Air Force Institute of Technology, Air University, Wright-Patterson Air Force Base, Ohio, 1976; GAW/ME/67B-3. (65) Shoemaker, R. L.; Stark, J. A.; Koshigoe, L. G.; Taylor, R. E. In Thermal Conductivity 18, 1st ed.; Ashworth, T., Smith, D. R., Eds.; Springer US, New York, 1985; Chapter 22: Thermophysical Properties of Propellants, pp 199–211, DOI: 10.1007/978-1-46844916-7. (66) Ellison, D. S.; Alcorn, R. A.; Neal, E. Effects of thermal cycling on trinitrotoluene and tritonal explosive compositions. J. Hazard Mater. 1980, 4, 57 – 75, DOI: 10.1016/03043894(80)80023-9. (67) Salo, J.; Cichinski, C.; Carignan, Y.; Turngren, E. Thermodynamic properties of energetic materials. III. Heat capacity of trinitrotoluene (TNT). 13th Symposium on Explosives and Pyrotechnics, Hilton Head Island, SC, USA. 1986; pp I–21 – I–24. (68) Miller, M. S. Thermophysical properties of RDX ; Technical Report ARL-TR-1319, Army Research Laboratory, Aberdeen Proving Ground, Maryland, 1997. 45 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

(69) Yang, L. C. Prediction of thermo-physical properties of energetic materials by similarity to inert substances. AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Nashville, TN, USA. 2010; pp 1 – 23, AIAA 2010-7006. http://arc.aiaa.org/doi/abs/10.2514/6.2010-7006. (70) Kishore, K. Comparative studies on the decomposition behaviour of secondary explosives RDX and HMX. Defence Sci. J. 1978, 28, 59 – 66, DOI: 10.14429/dsj.28.6535. (71) Zeman, S. Some predictions in the field of the physical thermal stability of nitramines. Thermochim. Acta 1997, 302, 11 – 16, DOI: 10.1016/S0040-6031(96)03101-2. (72) Hall, P. G. Thermal decomposition and phase transitions in solid nitramines. Trans. Faraday Soc. 1971, 67, 556 – 562, DOI: 10.1039/TF9716700556. (73) Nauflett, G. W.; Carlson, D.; Austin, T. D.; Brasch, J. W. Effect of pressure and melting behavior on the combustion of propellants containing RDX and RDX admixtures. 16th JANNAF Combustion Meeting, Naval Postgraduate School, Monterey, CA, USA. 1979; pp 95 – 111. (74) Jaansalu, K. M. Copyright permission obtained on January 8, 2016 through email. (75) Zerkle, D. K.; Núñez, M. P.; Zucker, J. M. Molten Composition B viscosity at elevated temperature. J. Energ. Mater. 2016, In press. (76) Baytos, J. F.; Craig, B. G.; Campbell, A. W.; Deal, W. E.; Dick, J. J.; Dinegar, R. H.; Engelke, R. P.; Larson, T. E.; Marshall, E.; Ramsay, J. B.; Rogers, R. N.; Soran, D. E.; Urizar, M. J.; Wackerle, J. D. LASL explosive property data, 1st ed.; University of California Press, London, 1980; Gibbs, T. R. and Popolato, A., Eds.; 479 pp. (77) Marsh, S. P. LASL shock Hugoniot data; University of California Press, Berkeley and Los Angeles, 1980; 674 pp. 46 ACS Paragon Plus Environment

Page 46 of 47

Page 47 of 47

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


Graphical TOC Entry 355

475

354

450 )

425 400

*

353 !" $ !%&'

#

352

!%'

375

351

(

350

350

!" $ !%&'

#

325

!%'

349

*

300

( (

348 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

0.01

0.02

!"

0.03


0.04

0.05

0.06

Application of the PengâRobinson Equation of ... - ACS Publications

Recommend Documents

Application of the PengâRobinson Equation of ... - ACS Publications