5070
J. Phys. Chem. B 2008, 112, 5070-5078
Microscopic Approaches to Liquid Nitromethane Detonation Properties Anaı1s Hervoue1 t, Nicolas Desbiens, Emeric Bourasseau, and Jean-Bernard Maillet* CEA, Centre DAM - Ile de France, De´ partement de Physique The´ orique et Applique´ e, Bruye` res-le-Chatel, 91297 Arpajon Cedex, France ReceiVed: September 10, 2007; In Final Form: January 14, 2008
In this paper, thermodynamic and chemical properties of nitromethane are investigated using microscopic simulations. The Hugoniot curve of the inert explosive is computed using Monte Carlo simulations with a modified version of the adaptative Erpenbeck equation of state and a recently developed intermolecular potential. Molecular dynamic simulations of nitromethane decomposition have been performed using a reactive potential, allowing the calculation of kinetic rate constants and activation energies. Finally, the Crussard curve of detonation products as well as thermodynamic properties at the Chapman-Jouguet (CJ) point are computed using reactive ensemble Monte Carlo simulations. Results are in good agreement with both thermochemical calculations and experimental measurements.
I. Introduction A classical description for the propagation of a planar detonation wave in a liquid explosive is provided by the Zel’dovich-von Neumann-Do¨ring (ZND) model. A shock wave propagates in the fresh explosive, associated with a sudden increase of pressure and temperature. The shock brings the system initially at rest onto a particular point of its Hugoniot curve (defined by the Rankine-Hugoniot relations), the ZND state. At this point, corresponding to high pressure and temperature, chemistry starts to occur, generally with endothermic reactions as a first step followed by exothermic reactions, which correspond to the formation of detonation products. This brings the system from the ZND state down to a point on the Crussard curve, which is still defined by the Rankine-Hugoniot relations but with different chemical composition. Only one point of this curve could support a stationary behavior of the reactive wave, the Chapman-Jouguet (CJ) point, where the system is supposed to be at chemical equilibrium. From this point, the system expands isentropically, with a decrease of both temperature and pressure. The stationary behavior of the detonation wave could then be partitioned in three arbitrary steps: In the first step, the fresh explosive is shocked up to the ZND point. In the second step, chemical reactions occur in the system, with the composition evolving with pressure and temperature along the Rayleigh line. In the last step, the system can be described by a mixture of detonation products at chemical equilibrium. These three steps are obviously linked together; for example, chemical reactions are thermally activated processes, and kinetic rates strongly depend on the temperature reached at the ZND point. The rate of chemical reactions, although not relevant for the prediction of detonation velocity, are important if one wants to describe the length of the reaction zone and the curvature of the detonation wave when the explosive is confined by another material. Finally, the CJ properties can be calculated by thermochemistry, in the assumption that the explosive has completely reacted (λ ) 1). The calculation of the CJ point and the Crussard curve is important when the detonation wave interacts with another inert material, so the pressure work * Corresponding author.
achieved by the explosive can be evaluated. In this paper we will successively describe microscopic approaches that have been used to study separately these three steps. In a first step, we want to describe thermodynamic properties of the inert explosive on its Hugoniot curve. Direct simulations of the propagation of a shock wave in a material are possible using nonequilibrium molecular dynamics (NEMD)1 but remain time-consuming, as stationary properties take time to establish. As we are not directly interested in nonequilibrium properties, equilibrium simulations appear to be adequate in this case. Several techniques have been developed in order to compute the Hugoniot curve of a system in a molecular dynamic framework.2-4 These methods do not require the explicit simulation of the shock propagation, thus saving a considerable amount of CPU time spent for the calculation of one Hugoniot point. These methods are conceptually similar, based on the idea of applying a constrain that maintains the system on its Hugoniot curve (or on its Rayleigh line). The constrain is applied through the use of modified equations of motion or by employing an extended Lagrangian. Within the Monte Carlo framework, the adaptative Erpenbeck equation of state (AEEOS) method, based on the Erpenbeck-EOS method,5 has been shown to be a valuable tool to compute a Hugoniot point in only one simulation.6 In this work, we have used the AE-EOS method in order to compute Hugoniot properties of nitromethane before any chemical reaction occurs (inert system). These simulations provide direct properties of liquid nitromethane on its Hugoniot, and particularly at the ZND conditions. These data could be used as input parameters to study the chemical behavior of nitromethane at these thermodynamics conditions, i.e., equivalent to the beginning of the reaction zone. The main objective of these reactive molecular simulations is to extract chemical descriptors, such as activation energy, that are used in a coarse grained model to take into account the evolution of chemical composition of the system in a statistical way. However, the activation energy for the decomposition reaction of an energetic material is naturally dependent on the reaction path followed by the system. Several different reaction paths could be in competition; some of them would be favored in a particular thermodynamic range while they would play no role
10.1021/jp077250n CCC: $40.75 © 2008 American Chemical Society Published on Web 04/01/2008
Liquid Nitromethane Detonation Properties in other thermodynamic conditions. In particular, concerning nitromethane, several decomposition mechanisms have been found, depending on density and temperature conditions. Indeed, since the C-N bond is the weakest bond of the nitromethane molecule, it has been thought for a long time that the decomposition mechanism would involve the breaking of this bond as a first step. This is the case, for example, for the decomposition in the gas phase, when unimolecular decomposition pathways are favored. Concerning condensed liquid phase in standard conditions, this mechanism seems to remain predominant as it is associated with the lowest activation energy. Nevertheless, as pressure (or density) increases, different mechanisms may appear, involving two or more nitromethane molecules. In particular, the scission of a C-H bond could lead to the production of the nitromethane aci ions CH2NO2- together with the CH3NO2H+ molecule. This has been found in ab initio simulations in condensed-phase crystals at high temperature and pressure near the CJ conditions7 and ab initio simulations of multimolecular nitromethane collisions.8 Mass spectroscopic studies in the chemical reaction zone of nitromethane9 conclude that condensation reactions between two or more nitromethane molecules occur during the detonation. In the same paper, a good overview is given about questions arising from nitromethane reactivity and decomposition mechanisms. Finally, it seems that the decomposition mechanism at high pressure involves several molecules, and becomes favored as the density increases (negative volume of activation). On the other hand, the decomposition mechanism associated with the lowest activation energy in standard conditions or in the gas phase seems to be unimolecular. The transition between these two identified mechanisms with density is not well-known, because they coexist on a particular thermodynamic domain. An exhaustive study of dominant decomposition mechanism with density is not possible at this time, but our aim is to calculate the decomposition reaction rate for specific thermodynamic conditions. For this purpose, ab initio simulations cannot be used because we are interested in describing condensed-phase chemistry on disordered (liquid) systems. However, this could be achieved with classical reactive potentials (“bond order potentials”), that allow covalent bond breaking and forming without the explicit treatment of the electrons. For this particular study we have used the ReaxFF potential developed at CalTech,10 which enables simulation of hundreds of atoms for hundreds of picoseconds. The decomposition mechanism can be extracted from this type of simulation,11,12 as well as the activation energy from cook-off simulations.13 The third part is concerned with the calculation of the thermodynamic properties (equation of state) of the detonation product mixture. Rankine-Hugoniot relations still apply to the system, the difference with the inert Hugoniot being in the evolution of the chemical composition of the system with temperature and pressure. The related curve is called the Crussard curve. This curve can be calculated only from thermochemical data, and knowledge of decomposition mechanism and kinetic rate is not required (this explains why thermochemical codes can predict detonation properties of energetic material). Only one point on this curve, the CJ point, leads to a stationary behavior of the detonation wave. It corresponds to the tangential point between the Rayleigh line and the Crussard curve. The slope of the Rayleigh line is proportional to the square of the detonation velocity. The detonation behavior of a particular explosive is then given by its properties at the CJ point. Other parts of the Crussard curve are still important to compute because they rely on the
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5071 “mechanical” behavior of the explosive. In a real detonation, the mixture of detonation products follows a rarefaction isentrope from the CJ point. This isentrope gives the amount of pressure (versus time) that can be transmitted to a connex material. Unfortunately, a microscopic expression of the entropy is not available, and the isentrope cannot yet be computed using an equilibrium constrained technique. However, in the approximation of a weak shock, the isentrope is similar to the Crussard curve up to the second derivative. The calculation of the Crussard curve below the CJ point would then give good indications about the pressure release in the reacted mixture. The difficulty at this point arises from the fact that the reacted mixture, although at chemical equilibrium, exhibits an evolution of the relative concentration of detonation products with temperature and pressure. Thus it cannot be simulated using a fixed concentration of chemical species. In order to overcome this obstacle, we have made use of the reactive Monte Carlo method, which allows one to simulate the chemical equilibrium of a system. The reactive Monte Carlo method has been coupled to the AE-EOS method such that the Crussard curve can be directly evaluated.14 The paper is organized as follow: first, the different techniques used are described in the next section. Only a brief description of existing techniques is given with related references, and particular attention is given to the coupling of different techniques (which could lead, in some cases, to numerical instabilities). Then results are presented for nitromethane, concerning its unreacted Hugoniot curve, its chemical decomposition behavior, and finally the Crussard curve of its detonation products. Discussion and conclusions are provided at the end. II. Methods and Models A. The Inert Hugoniot. 1. The AE-EOS Method. In order to compute the Hugoniot curve of the inert explosive, we employed the AE-EOS method proposed by Brennan et al.6 As discussed in this reference, it is possible to implement this method in several fashions. We did not follow the idea of Brennan, who chose a succession of NPT simulations to converge to the Hugoniot pressure at a given temperature. We preferred to use a succession of NVT simulations to converge to the Hugoniot temperature at a given compression. This method, and its advantages, is detailed in reference 14. The main idea of the AE-EOS method is to minimize expression 1, which measures the gap between the simulated thermodynamic state and the real Hugoniot state. In this equation, E0, P0, and V0 are the energy, pressure, and volume, respectively, of the inert explosive before the detonation (the pole of the Hugoniot).
1 Hg ) E - E0 - (P + P0)(V0 - V) 2
(1)
So, starting from an initial configuration, the system is simulated in the canonic ensemble (NVT1) and the Hugoniot difference Hg(1) is evaluated following expression 1. Then, the temperature of the simulation is slightly modified from T1 to T2 ((10 K). The simulation NVT2 carries on, and Hg(2) is evaluated. At the end of this second step, the following derivative is also evaluated:
Hg(1) - Hg(2) dHg (2) ) dT T1 - T 2
(2)
From here, the new temperature of the simulation T3 is calculated through
5072 J. Phys. Chem. B, Vol. 112, No. 16, 2008
T3 ) T2 -
Hervoue¨t et al. TABLE 1: Geometry of Nitromethane
Hg(2) dHg (2) dT
(3)
The simulation NVT3 carries on, and Hg(3) is evaluated. This process is automatically iterated every 5 × 105 steps, evaluating each time
Hg(n - 1) - Hg(n) dHg (n) ) dT Tn-1 - Tn
(4)
and calculating a new temperature using
Tn+1 ) Tn -
Hg(n) dHg (n) dT
(5)
until the Hugoniot difference Hg(n) has converged to a required accuracy. Once the Hugoniot temperature THug is reached, an additional NVTHug simulation could be performed in order to compute the accurate pressure PHug. The advantage of this method is its rapidity. The convergence on the Hugoniot curve is usually obtained after 5 × 106 steps only in the NVT ensemble. Using this root-finding algorithm, we managed to calculate the Hugoniot curve of nitromethane. 2. Simulation Details and Molecular Models. All the simulations performed in this part of work have been done with the program GIBBS, owned by the Institut Franc¸ ais du Pe´trole, the Universite´ Paris Sud, and the CNRS, and developed in collaboration between those three owners and the CEA.15 All the simulations were performed with 150 nitromethane molecules. The respective probabilities of choosing a translation move, a rotation move, or a volume change were 0.48, 0.48, and 0.04. The convergence of AE-EOS simulations were generally obtained after 3 × 106 iterations, changing the temperature every 2 × 105 iterations and calculating the Hg value on 1 × 105 iterations before evaluating the new temperature. The potential has been developed especially for this study. It consists of intermolecular interactions modeled by EXP-6 type potentials and point charges. The CH3 group is coarse-grained, and the molecule is kept in a fixed geometry. As it will be seen, flexible models are not necessary in our case. Both force centers and charges are located on the N, O, and C atoms. The interaction energy between two molecules i and j is calculated as (a and b are N, O, and CH3)
Uij(rab) )
∑a ∑b
{ ( ab
Rab - 6
6e
Rab 1-{rab/σab)
( ))
- Rab
σab
6
+
rab
q a qbe 2 4π0rab
}
(6) The ab, σab, and Rab parameters are expressed using the following mixing rules:
ab ) xab Rab ) xRaRb σab )
σa + σb 2
(7)
The reaction field method has been used to evaluate the electrostatic part of the energy, with r ) 100. The molecular bond lengths and angles of nitromethane are given in Table 1. These values are taken from ref 16, except
dNO dNC dCM ONO
1.225 Å 1.49 Å 0.21584 Å 125.0°
TABLE 2: Potential Parameters N O CH3
σ (Å)
(K)
a
q (e)
3.28714 2.79701 4.53075
59.3182 95.4856 99.1565
15.0766 12.4018 13.1211
+0.54 -0.37 +0.20
dCM. M corresponds to the position of the center of force for the methyl group and is placed on the CN axis toward the hydrogen atoms. This value corresponds to the AUA displacement proposed by Ungerer et al. to model the CH3 group in alkanes.17 The values of the potential parameters are given in Table 2. The values of the repulsion-dispersion parameters have been obtained by an optimization procedure detailed in ref 18. The charges are the average charges deduced by a Bader analysis of electronic density maps of liquid nitromethane.19 B. The Reactive Part. In order to model the chemical behavior of nitromethane, we used the updated version of ReaxFF exposed in ref 20. The different contributions of the potential are not discussed here; the parameter fitting was performed on a large database including specific data on highenergy materials. Cook-off simulations of liquid nitromethane were performed for temperatures between 2000 and 3000 K and for a density of 1.97 g‚cm-3. The initial configuration is the one obtained after a NVT simulation in the liquid state at standard conditions (the liquid is stable at room temperature). The simulation box contains 32 molecules, and periodic boundary conditions are used. Simulation in the canonical (N, V, T) ensemble are performed using a Berendsen thermostat. Equations of motion are integrated with a time step of 0.1 × 10-15 s. The chemical analysis of the system is performed using a criterion based on the bond order parameter. In most cases, the threshold between bonded and nonbonded interaction is chosen to be 0.3. For a dense system, the bond order criterion is increased (a finer criterion such as the one used in ref 11, considering the relative atomic velocities, has not been used here). The concentration of different chemical species is then collected as decomposition occurs. Using the definition of the kinetic rate of reaction, different fitting forms corresponding to different reaction order are tested against the evolution of an explosive’s concentration with time. The rate constant k can then be computed for the given thermodynamic conditions (F, T). Repeating the same procedure for different initial temperatures allows the calculation of the activation energy in an Arrhe´nius framework, as well as the pre-exponential factor. We can input these microscopic data to compute the width of the (macroscopic) reaction zone for a one-dimensional detonation. Usually, the microscopic activation energy does not match the macroscopic length of the reaction zone. This can be explained by involving a mesoscopic effect such as hot spot formation. Indeed, the microscopic activation energy would lead to a correct prediction of the reaction zone only in the case of a perfect, homogeneous detonation. Unfortunately, this is never the case, and mesoscopic effects have to be taken into account.21 A common solution to avoid this type of problem is simply to consider a global effective kinetic rate constant, fitted to reproduce detonation properties. This is obviously not our goal here, where microscopic methodologies are used in a bottomup fashion in order to predict some properties related to
Liquid Nitromethane Detonation Properties
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5073
detonation when possible and/or feed coarse grain (hydrodynamic) models otherwise. C. The Crussard Curve. To obtain the Crussard curve of nitromethane, we used the combination of two recent methods: the AE-EOS method to obtain thermodynamic points satisfying Ranking-Hugoniot equations, and the reaction ensemble Monte Carlo (ReMC) method to reach the chemical equilibrium of the mixture of products from nitromethane detonation. 1. The ReMC Method. The goal of the ReMC method is to obtain macroscopic properties of a multicomponent system constrained by chemical equations. To assume this constraint, a particular statistical ensemble is defined: the reaction ensemble. The complete definition of this ensemble, and the rigorous way to obtain its density probability, were first described by Smith and Triska.22 Johnson has also given a very nice derivation of the reaction ensemble acceptance probability in reference 23. It is also interesting to consider the pioneer work of Shaw, who proposed earlier a method similar in nature to simulate chemical equilibrium of molecular mixtures.24,25 Here, we will only give a brief description of this ensemble, and we invite interested readers to see a previous article we wrote about the microscopic calculation of Crussard curves for more details.14 We consider a system of s different chemical species ai. In the reaction ensemble, the temperature and the volume (or pressure) are fixed. In comparison with the canonical ensemble, two more constraints are applied to satisfy the chemical equilibrium: • the number of atoms is fixed for each atom type, and it is possible to link the number of different molecules in the system, following the chemical equation that defines the chemical equilibrium: s
νiai ) 0 ∑ i)1
(8)
It can be shown that, in the case of these usual moves, only the energy changes between the two configurations. We obtain
Pacc ) min(1, exp(-β∆U))
where ∆U ) U(new) - U(old). To obtain the chemical equilibrium, we used the so-called reaction move, proposed at the same time by Smith et al.22 and Johnson et al.23 This move consists of first choosing a direction to perform the reaction, second, deleting a set of reactant molecules randomly chosen in the system, and finally inserting product molecules. Following the example given before, the reaction move can be
2NH3 f N2 + 3H2
2NH3 a N2 + 3H2
(9)
• the sum of chemical potentials over the different molecule species implied in the chemical reaction, weighted by stoichiometric coefficients, is equal to zero: s
νiµi ) 0 ∑ i)1
(10)
µN2 + 3µH2 ) 2µNH3
(11)
for example,
The key point is to establish a Metropolis algorithm to assume those constraints over the simulation. To obtain the mechanical equilibrium of the system, the usual Monte Carlo moves can be used: translation, rotation, and internal relaxation of molecules. The probability to accept a move is given by26
( (
Pacc ) min 1, exp
))
Fens(new) Fens(old)
(12)
where Fens(new) (and, respectively, Fens(old)) is the density probability of the new (old) configuration in the statistical ensemble.
(14)
when two randomly chosen NH3 molecules are deleted on one hand, and one N2 molecule and three H2 molecules are inserted on the other hand, or
N2 + 3H2 f 2NH3
(15)
when one N2 molecule and three H2 molecules randomly selected are deleted, and two NH3 molecules are inserted. During this move, not only is the energy modified, but the respective number of molecules is also implied in the reaction. Introducing the parameter ξ, which is positive if the reaction move is performed in the forward direction, or negative if performed in the reverse direction, it can be shown that the acceptance probability is given by (see ref 14 for details)
(
(
Pacc ) min 1,(P0βV) ‚exp -ξ ξjν
s νi∆fG0i (T) ∑i)1
RT Ni!
s
)
‚
exp(-β∆Uext) ∏ i)1 (N + ξν )! i
for example,
(13)
i
)
(16)
s where νj ) ∑i)1 νi, ∆fG0i (T) is the standard free enthalpy of formation of i at T, Ni is the number of molecules i in the system, and ∆Uext is the intermolecular energy difference between the old and new configurations. It is important to note that this expression is true only if rigid molecules are used (i.e., molecules with internal energy equal to zero). This is our case in this work when performing ReMC simulation of detonation product mixtures modeled by simple EXP-6 potentials. Practically, when a reaction move is chosen, the algorithm is the following: • The direction of the move is randomly chosen (ξ ) 0 or 1). The types of the R reactant molecules (to be deleted) and the P product molecules (to be inserted) are defined. • νi molecules for each reactant type of molecule are randomly chosen and deleted from the configuration. • νi molecules of each product type of molecule are randomly inserted in the configuration (i.e., the place of insertion and the orientation of the inserted molecule are randomly chosen). • The energy of the new configuration is calculated, and the acceptance probability of the move is obtained using eq 16. • Following the Metropolis algorithm, the new configuration is accepted or refused in the Markovian chain comparing the value of Pacc with a real randomly taken between 0 and 1. Combining translation, rotation, and reaction moves, it is possible to simulate a chemical equilibrium at a given temperature and a given density. It is also possible to use statistical
5074 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Hervoue¨t et al.
biases to solve the problem of insertion in dense phase. As an example, Brennan proposed the use of the cavity bias sampling to improve the insertion of molecules during a reaction move.27 We preferred to use the pre-insertion bias,28 already implemented in our Monte Carlo code, but the two algorithms are quite similar. Those biases applied to the reaction move consists first of all of inserting the first product molecules at the locations released by the deleted reactant molecules. If other product molecules have to be inserted (i.e., if νj > 0), the insertion is performed at a preselected location. The details of the algorithm and the expression of the acceptance probability when using the pre-insertion bias are presented in ref 14. It is interesting to note that the only input data needed in the reaction move is the ∆rG0(T) of the simulated reaction. This standard free enthalpy of reaction is obtained from an experimental database. It is also possible to compute a chemical equilibrium combining several chemical reactions. To do this, we used the same algorithm as presented here, but we added a preliminary step consisting of randomly choosing the reaction to perform before each reaction move. Finally, note that it is also possible to simulate a chemical equilibrium with a constant pressure. To do this, another move has to be used during the simulation: the usual volume change Monte Carlo movement. As in the NPT ensemble, it consists of modifying the volume of the simulation box, performing an homothetic transformation of the system. The acceptance probability of this move is given by29
[ ( (
( )))]
Vnew Pacc ) min 1, exp -β ∆U + P∆V + N ln Vold
3. The AE-EOS ReMC Method. When using the ReMC method with the composite ensemble, it is difficult to use the AE-EOS with a succession of NVT simulations as presented in section (2.1). Indeed, to fix the volume of the system (i.e., the volume of the fluid phase plus the volume of the virtual solid phase) implies that the volume of the fluid phase would be changed when a reaction move involving a carbon atom is performed. This presents some avoidable algorithmic complications. As a consequence, we used the AE-EOS method with a succession of NPT simulations. Nevertheless, two different algorithms have been used. Some simulations have been done with a succession of NPT simulations to converge to the Hugoniot temperature at a given pressure, and others have been done with a succession of NPT simulations to converge to the Hugoniot pressure at a given temperature. Results show that the two methods are similar. It is important to understand that the ReMC AE-EOS method used with the composite ensemble tries to minimize the following expression, similar to eq 1, but where the properties of the solid phase are introduced:
1 Hg ) Etot - E0 - (P + P0)(V0 - Vtot) 2 with
Etot ) Efluid + ECsol Vtot ) Vfluid + VCsol
(17)
2. The Composite Ensemble. The calculation of the Crussard curve of Nitromethane shows a particular difficulty in comparison with other explosives (see ref 14 for examples). Indeed, experiments show that a solid phase of carbon appears in the detonation product mixture. Unfortunately, we are not actually able to explicitly simulate a chemical equilibrium involving a solid phase. Nevertheless, following the work described in our previous article,14 we have developed a new method to implicitly take into account the solid phase of carbon in our calculations. This method is strongly inspired by the work of Shaw published in 2002 and is called the composite Monte Carlo method.30 In our method, we consider that the solid phase of carbon is taken into account in a virtual box containing NCsol carbon atoms, the properties of which (volume and energy) are given by an equation of state regarding the equilibrium pressure of the system. To assume the chemical equilibrium between the solid and the fluid phases, a particular reaction move must be added: a reaction implying a carbon atom as, for example, 2CO a Csol + CO2. When this particular Monte Carlo move is chosen, 2 CO molecules are deleted from the liquid mixture simulation box, and a O2 molecule is randomly inserted. At the same time, NCsol is increased by one, and the properties of the solid phase are modified as a consequence. Obviously, the reverse move is also performed from time to time, in order to respect the microreversibility. The acceptance probability of this move can be computed using eq 16. Nevertheless, as the intermolecular energy of the solid phase of carbon is neglected in the simulation, this must be included in the ∆fG0Csol(T). Thus, this data corresponds, in fact, to the chemical potential of the virtual solid phase of carbon in the given conditions, calculated from the equation of state.
(18)
(19)
4. Simulations Details. All the simulations performed in this part of our work have also been done with the program GIBBS. To reach the chemical equilibrium of a system using the ReMC method, it is necessary to define a priori the set of chemical equations driving the system. To do this, the first step is to choose which type of molecule will be present a the equilibrium. In our case, it is known that the detonation products of nitromethane are the following: CO2, H2O, CO, N2, H2, O2, NO, NH3, and CH4 in a fluid phase, and C in a solid phase. Nevertheless, it is known that O2 and NO exist in really small quantity, usually smaller than 0.1%. As a consequence, we decided to neglect those two types of molecules, and finally the simulated system is composed of seven molecular types in the fluid mixture: CO2, H2O, CO, N2, H2, NH3, and CH4, and solid carbon in a virtual phase. The second step consists of finding a set of linearly independent chemical equations involving those molecules, and satisfying the stoichiometry and the chemical specificity of the system. To do this, we applied the formula-vector matrix stoichiometric algorithm devised by Smith and Missen.31 This method leads to the following set of reactions:
2NH3 a N2 + H2 CO + 2NH3 a N2 + CH4 + H2O CO2 + H2 a CO + H2O 2CO a Csol + CO2
(20)
To described the molecules in the fluid mixture, we used the exponential-6 potential proposed by Fried et al.32 and optimized to reproduce the thermodynamic behavior of pure compounds under high temperature and high pressure. A cutoff equal to
Liquid Nitromethane Detonation Properties
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5075
Figure 1. Hugoniot curves of nitromethane. Open symbols are experiments, blue circles are from ref 40, red circles are from ref 41, and green circles are from ref 42. Black dots are the results of our simulations. The dashed line is a guide to the eyes.
Figure 2. Pressure vs temperature for the shock Hugoniot of nitromethane. Orange line: results from Lysne and Hardesty;43 red circles: results from Jones;38 blue circles: results from Soulard;37 green circles: results from Liu et al.;39 and black dots: results from our simulations.
half of the box length has been used to reduce the computing time. Limit conditions and long-range corrections have been applied. To described the virtual solid phase of NCsol carbon atoms, we have used the Ree-van Thiel equation of state.33 This allows us to compute the energy, volume, and chemical potential of the solid phase of carbon in regard to NCsol, P, and T. All the simulations performed in this part of our work were initiated with 250 NH3 molecules and 250 CO2 molecules. This corresponds to a system of 250 nitromethane molecules. The respective probabilities of choosing a translation move, a volume change, and a reaction move were 0.8, 0.02, and 0.18, respectively. The convergence of AE-EOS ReMC simulations were generally obtained after 107 iterations, changing the pressure every 0.5 × 106 iterations and calculating the Hg value of 0.25 × 106 iterations before evaluating the new pressure.
TABLE 3: Main Characteristics of the Force Fields Used in MD Works
III. Results A. The Inert Hugoniot. To compute the Hugoniot curve, one needs to know the initial state. In our case, the initial state is chosen as P0 ) 1 atm and T0 ) 298 K. To determine the specific volume of nitromethane in the initial state, we performed an NPT simulation of 150 molecules, which gave us a density of 1144 kg‚m-3 and a vaporization enthalpy of 36.59 kJ‚mol-1. Concerning the density, this corresponds to a difference of 0.6% with the experimental values found in literature: 1137 kg‚m-3 from ref 34, 1134 kg‚m-3 from ref 35, and 1140 kg‚m-3 from ref 36. The value of V0 is thus 53.32 × 10-6 m3‚mol-1. Concerning the vaporization enthalpy, the difference is about 4.4% with the experimental value (38.31 kJ‚mol-1 from ref 34). Once we have calculated the initial state, we can simulate the whole Hugoniot line. The results are reported in Figure 1. The agreement between our simulation results and the experimental data is very satisfactory over the whole range of pressure. Other Hugoniot curves of nitromethane, based on molecular dynamics (MD) simulations, have been reported in the literature.37-39 The different force fields used in MD simulations are rather more complex than our own. For example, they all contain intramolecular interactions (bending, stretching, etc.), as can be seen in Table 3. Although recent, the force field developed by Liu et al.39 seems to be the poorer one since the simulated density at ambient conditions is more than 3% lower than the experimental one, and the simulated isotherm at 300
this work ref 37 ref 38 ref 39
RPa
coulb
X X X X
X X X
stretchc
bendc
torsc
X X X
X X X
X X
a RP stands for repulsion-dispersion part. b Coul stands for coulombic interactions. c Stretch, bend and tors correspond to intramolecular stretching, bending and torsion interactions, respectively.
K is far from the experiments above 2 GPa. On the contrary, the force fields developed by Soulard37 and Jones38 and our own force field are much more accurate. Concerning the simulated P-V Hugoniot curves (not shown here), the results obtained by Soulard and Jones are in good agreement with experiments, while the ones of Liu are markedly different. On an other hand, simulated shock temperatures are a good indicator of the quality of a force field. The prediction of the shock temperature is an important issue in the studies of energetic materials. The results obtained by Soulard, Jones, and Liu et al.37-39 are compared to our results and to the Lysne and Hardesty results43 (experimental results and physical equation of state) in Figure 2. For a given pressure, temperatures obtained by MD simulations are lower than the experimental data and the discrepancy increases with the pressure. Our results are in good agreement with Lysne and Hardesty’s results. The slope of the curve is very close from the experimental one, and there’s only a difference of 100 K at 10 GPa. Soulard37 suggested to improve the description of the intramolecular part of the potential, and Jones38 suggested to take into account the quantum effects (instantaneous polarization, transfer, correlation, exchange, etc.). Our results clearly show that an adequate functional form (EXP-6) for the repulsion-dispersion interaction, point charges, and a rigid molecule can predict accurate quantities along the Hugoniot. Finally, the us(up) curve has been studied (but not shown here). Our results and the results of Soulard and Jones compare very well to Lysne and Hardesty’s data. Second shock properties of nitromethane have also been examined and compared to experimental results provided by Crouzet.40 The results are reported in Table 4. It appears that the agreement between experimental and simulated properties is very satisfactory, indicating the transferability of the potential
5076 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Hervoue¨t et al.
TABLE 4: Experimental and Simulated Properties of Reshocked Nitromethane
a
data
experimentsa
simulations
P2ndshock (GPa) us (m‚s-1) T (K)
14.6 875 980 ( 50
16.13 ( 0.04 925 ( 7 1057
Reference 40.
over a broad range of pressure. It is important to note that, to our knowledge, no other potential is able to accurately reproduce shock temperatures and even more second shock temperatures. As a consequence, shock temperature appears to be relevant to evaluate the quality of a potential. B. The Reactive Part. In the gas phase, it is known that the mechanism of decomposition of nitromethane involves the breaking of the C-N bond. Indeed, this bond is the weakest bond of the molecule and, consequently, the easiest to break. Moreover, the associated activation energy of 234.3 kJ‚mol-1 is almost equal to the bond energy. Analysis in a condensedphase system is more complicated. A good review of experimental results is provided in ref 9. As pressure or density increases, different decomposition mechanisms are activated. One in particular involves the presence of the aci ion. Experiments described in this paper concern the mass spectroscopy analysis of reaction products in the reaction zone of a detonation wave propagating in nitromethane. These experiments evidence species resulting from condensation reactions implying several nitromethane molecules, corresponding to at least a bimolecular mechanism. On the theoretical side, Manaa and co-workers performed DFT simulations of early thermal decomposition of dense and hot nitromethane.7 The observed decomposition process consists in the rapid formation of the two ions CH3NO2H+ and CH2NO2-, coming from the ejection of an hydrogen atom from one molecule. Interestingly, this process involves three molecules. In condensed phase, and particularly at high density, a caging effect, i.e., the constraints coming from the surrounding molecules to the intramolecular degrees of freedom, may enhance vibrations in the C-H mode. Eventually these vibrations could lead to the breaking of the covalent bond. In our simulations we observed similar mechanisms, i.e., a high concentration of CH3NO and CH2NO as intermediate products of reactions. These molecules are the leading products of the deprotonation process of one nitromethane molecule. The proton is then captured by a surrounding oxygen atom, eventually pertaining to the same molecule. This leads to the formation of CH3NOOH or CH2NOOH. These molecules then quickly release their OH group, which eventually captures another hydrogen atom in order to form one water molecule. Indeed, a water molecule is not only the most stable decomposition product but is also formed very rapidly. The mechanism can involve one or two molecules. In Figure 3 the concentration of nitromethane molecule is plotted as a function of time for different temperatures (T ) 2200, 2400, 2600, 2800, and 3000), together with a fit corresponding to first-order chemical reactions. Different fitting functions have been tested (order zero up to second order). The initial temperatures of the system determines the time required for the decomposition to occur. These relatively high temperatures allow a direct observation of nitromethane decomposition in MD simulations. For each temperature, the effective rate constant of nitromethane chemical decomposition can be extracted from the fit. “Effective” has a statistical meaning in the sense that it corresponds to an “observable” and not to a specific route of
Figure 3. Evolution of the concentration of nitromethane molecules (normalized by the initial concentration) with time, for temperatures equal to 2200 K (blue curve), 2400 K (cyan), 2600 K (green), 2800 K (orange), and 3000 K (red).
TABLE 5: Properties of the Detonation Product Mixture Including the Solid Phase of Carbon along the Crussard Curve Obtained Using the AE-EOS ReMC Method and Compared with Thermochemical Results thermochemical results P (GPa) 10 12 13.75 T (K) 3500 3700
Monte Carlo results
THug (K)
VHug (cm3‚g-1)
THug (K)
VHug (cm3‚g-1)
3279 3443 3581 PHug (GPa) 12.81 15.32
0.7123 0.6711 0.6431 VHug (cm3‚g-1) 0.6566 0.6207
3325 3488 3625 PHug (GPa) 12.15 14.72
0.7025 0.6583 0.6286 VHug (cm3‚g-1) 0.6555 0.6146
decomposition. This effective rate constant is an average of the rate constants associated with different (activated) decomposition mechanisms. In the Arrhenius framework, the activation energy is the slope of the evolution of the rate constant with temperature. In this case the activation energy associated with the chemical decomposition of nitromethane is equal to 136.8 kJ‚mol-1. This result appears to be much lower than the activation energy corresponding to the breaking of the C-N bond at low density or pressure. This confirms the fact that the activation energy for a particular mechanism varies with density. It is particularly the case when bimolecular processes are involved, as they are expected to be greatly favored as density increases. C. The Crussard Curve. Using the method described in section 2.3, we have performed five AE-EOS ReMC simulations for the following five pressures and temperatures: 10, 12, and 13.75 GPa, and 3500 and 3700 K. Of course, in those simulations we used the same pole conditions used in part 3.1 to calculate the inert Hugoniot of nitromethane. The converged Hugoniot temperatures THug and corresponding volumes VHug for the three first simulations, and the converged Hugoniot pressures PHug and corresponding volumes VHug for the two other simulations are presented in Table 5. The Crussard curve of nitromethane detonation products is presented in Figure 4 and compared to the inert Hugoniot of nitromethane calculated in part 3.1, as well as to the crussard curve of nitromethane calculated neglecting the solid phase of carbon. Table 5 and Figure 4 show that the properties calculated along the Crussard curve using our Monte Carlo method are in good agreement with thermochemical results. We can also see that the effect of the solid phase of carbon is important, especially when high pressures are involved.
Liquid Nitromethane Detonation Properties
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5077 important to notice that, considering detonation velocity, the Monte Carlo result is closer to experiment than the thermochemical one. IV. Conclusion
Figure 4. Inert Hugoniot of nitromethane (circles) and Crussard curves with (squares) and without (triangles) the solid phase of carbon calculated using the AE-EOS ReMC method, compared with thermochemical results. Open symbols are the CJ points calculated through Monte Carlo simulations. The star is the thermochemical CJ point.
TABLE 6: CJ States of the Nitromethane Calculated With or Without the Solid Phase of Carbona PCJ VCJ (GPa) (cm3‚g-1) Monte Carlo results with carbon 11.84 Monte Carlo results without carbon 12.13 thermochemical results 12.31 experimental results
0.66 0.668 0.653
TCJ (K)
DCJ (m‚s-1)
3430 3274 3513
6438 6641 6493 6300b
a Monte Carlo results are compared to thermochemical results and, considering the detonation velocity, to experimental measurements. b Reference 40.
Once the Crussard curve has been calculated, the CJ point can be determined. To do this, we used the fitting procedure presented in ref 14. The CJ points calculated from the results obtained with and without the solid-phase of carbon and from the thermochemical results are presented in Table 4. From the calculated PCJ and VCJ, we have also determined the detonation velocity DCJ, using the following equation derived from the Hugoniot relations:
DCJ ) V0
x
(PCJ - P0) (V0 - VCJ)
(21)
This property is interesting because it can directly be measured experimentally. Table 4 shows a comparison between the calculated and measured value of DCJ. Table 6 shows that the CJ state calculated using our Monte Carlo method and with the solid-phase of carbon are in good agreement with thermochemical results. The discrepancies between the two methods represent 3.8% concerning pressures, 1.0% concerning volumes, 2.3% concerning temperatures, and 0.8% concerning detonation velocities. We can also note that the calculated value of DCJ is in good agreement with the experimental value. The difference between the Monte Carlo simulation result with the solid phase of carbon and the experimental one is about 2.2%, which is very satisfactory considering that simple potential models have been used. It is already anticipated that these differences would vanish as realistic potential models are used. On the contrary, when neglecting the solid phase of carbon, the error concerning the detonation velocity is more significant, about 5.5%. It is also
In this paper, we have used several atomistic methods to compute the different stages encountered during a detonation. In particular, the Hugoniot curve of the inert explosive has been computed with a good accuracy using a recently developed intermolecular potential for nitromethane. Moreover, second shock properties including temperature measurements have been computed and compared directly to experimental data. The good agreement between our calculations and experimental data is a good indication of the transferability of the potential over a large range of pressure. Second, reactive properties of nitromethane have been computed using a reactive force field. Cook-off simulations at different densities have shown that the breaking of the C-N bond is not the decomposition mechanism in condensed phase reactions. Instead, we observed the breaking of the C-H bond leading to the appearance of CH2NOOH and CH3NOOH intermediates. Water is then the first decomposition product to be formed. The calculated activation energy confirms that the decomposition mechanism exhibits a negative volume of activation. There is no experimental data for the kinetic rate of decomposition at these extreme conditions of pressure and temperature, and a validation at other conditions does not appear to be appropriate (as the mechanism of decomposition depends on thermodynamic conditions). Nevertheless, a global validation remains possible; for that, one should use the kinetic parameters extracted from our reactive simulations in a reactive hydrodynamic simulation of detonation. Results concerning the curvature of the wave or critical diameter could then be compared to their experimental counterparts. Finally, the Crussard curve of the detonation products mixture at chemical equilibrium has been calculated using a simple potential model and including a carbon solid phase. Results for the Crussard curve are in good agreement with thermochemical calculations. The CJ point can be determined from the Crussard curve, and an additional simulation at CJ conditions has shown good agreement with experimental data. From this study, we see that microscopic methods are well suited for the prediction of detonation properties of liquids, including inert, reactive, and reacted thermodynamic properties of explosives. Results obtained so far compared well with experimental data and thermochemical calculations, and serve as a basis for future developments. First, the potential models used for the simulations of detonation products are the same as those used in thermochemical calculation, i.e., simple models. Although these potentials may be accurate enough to model the behavior of nitrogen under pressure, we also anticipate that more complex potentials may be required to model the behavior of polar molecules (such as water and hydrogen fluoride) over a large range of temperatures and pressures. These more complex potentials could easily be included in our simulations. Second, the methodology proposed in this paper is well suited for the calculation of homogeneous properties. When heterogeneous structure is considered, as, for example, in solid explosives, the mesoscopic structure of the solid should be taken into account. Indeed, for solid explosives, mesoscopic defects lead to the appearance of hot spots that control the detonation process. Simulations of these processes could be envisaged through a mesoscopic modelization of a reactive system.44
5078 J. Phys. Chem. B, Vol. 112, No. 16, 2008 Acknowledgment. All Monte Carlo simulations have been performed with the Gibbs code from IFP, CNRS, and the Universite´ Paris-Sud.15 References and Notes (1) Holian, B. L.; Lomdahl, P. S. Science 1998, 280, 2085. (2) Soulard, L. In Shock Compression of Condensed Matters1999, Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Snowbird, UT, June 27-July 2, 1999; Furnish, M. D., Chhabildas, L. C., Hixson, R. S., Eds.; AIP Press: New York, 2000. (3) Maillet, J.-B.; Mareschal, M.; Soulard, L.; Ravelo, R.; Lomdahl, P. S.; Germann, T. C.; Holian, B. L. Phys. ReV. E. 2001, 63, 016121. (4) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. Phys. ReV. Lett. 2003, 90, 235503. (5) Erpenbeck, J. J. Phys. ReV. A 1992, 46, 6406. (6) Brennan, J. K.; Rice, B. M. Mol. Phys. 2003, 101, 3309. (7) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Galli, G.; Gygi, F. J. Chem. Phys. 2004, 120, 10146. (8) Decker, S. A.; Woo, T. K.; Wei, D.; Zhang, F. In Proceedings of the 12th International Detonation Symposium; Office of Naval Research: Arlington, VA, 2002. (9) Blais, N. C.; Engelke, R.; Sheffield, S. A. J. Phys. Chem. 1997, 101, 8285. (10) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Chem. Phys. 2001, 105, 9396. (11) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Phys. ReV. Lett. 2003, 91 (9), 098301. (12) van Duin, A. C. T.; Zeiri, Y.; Dubnikova, F.; Kosloff, R.; Goddard, W. A. J. Am. Chem. Soc. 2005, 127, 11053. (13) Maillet, J.-B.; Soulard, L.; Germann, T. C.; Strachan, A.; Van Duin, A. Computational Modeling and Simulation of Materials III, Proceedings of the 3rd International Conference on Computational Modeling and Simulation of Materials (CIMTEC), Acireale, Italy, May 30-June 4, 2004; Techna Press: Faenza, Italy, 2004; Part A, p 347. (14) Bourasseau, E.; Dubois, V.; Desbiens, N.; Maillet, J.-B. J. Chem. Phys. 2007, 127, 084513. (15) Ungerer, P.; Boutin, A.; Tavitian, B. Applications of Molecular Simulation in the Oil and Gas Industry. IFP Publications: Rueil Malmaison, France, 2005. (16) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. J. Comput. Chem. 2001, 22, 1340. (17) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Fuchs, A. J. Chem. Phys. 2000, 112, 5499. (18) Desbiens, N.; Bourasseau, E.; Maillet, J.-B. Mol. Simul. 2007, 33, 1061. (19) Bourasseau, E.; Maillet, J.-B. In Proceedings of the APS Conference on Shock Compression of Condensed Matter, Baltimore, MD, July 31August 5, 2005; AIP Press: New York, 2005. (20) Nielson, K. D.; Van Duin, A. C. T.; Oxgaard, J.; Deng, W.-Q.; Goddard, W. A. J. Phys. Chem. A 2005, 109, 493.
Hervoue¨t et al. (21) Stoltz, G.; Maillet, J.-B.; Soulard, L. To be submitted for publication. (22) Smith W. R.; Triska, B. J. Chem. Phys. 1994, 100, 3019. (23) Johnson, J. K.; Panagiotopoulos, A. Z.; Gubbins, K. E. Mol. Phys. 1994, 81 (3), 717. (24) Shaw, M. S. J. Chem. Phys. 1991, 94, 7550. (25) Shaw, M. S. In Proceedings of the 12th Symposium of Detonation, San Diego, CA, August 11-16, 2002; Office of Naval Research: Arlington, VA, 2002. (26) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: New York, 2002. (27) Brennan, J. K. Mol. Phys. 2005, 103, 2647. (28) Bourasseau, E.; Ungerer, P.; Boutin, A. J. Phys. Chem. B 2002, 106, 5483. (29) McQuarrie, D. Statistical Mechanics; Harper and Row: New York, 1976. (30) Shaw, M. S. In Proceedings of the AIP Conference on Shock Compression and Condensed Matter, Atlanta, GA, June 24-29, 2001; AIP Press: New York, 2001. (31) Smith, W. R.; Missen, R. W. Chem. Eng. Educ. 1979, 13, 26-32. (32) Fried, L. E.; Howard, W. M.; Souers, P. C. In Proceedings of the 12th International Detonation Symposium, San Diego, CA, August 1116, 2002; Office of Naval Research: Arlington, VA, 2002. (33) van Thiel, M.; Ree, F. H. Phys. ReV. B 1993, 48 (6), 3591. (34) Lide, D. R. Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1997. (35) Leal-Crouzet, B.; Baudin, G.; Presles, H. N. Combust. Flame 2000, 122, 463. (36) Riddick, J. A.; Bunger, W. B.; Sakano, T. K. Techniques of Chemistry, Organic SolVents, Physical Properties and Methods of Purification; Wiley: New York, 1986. (37) Soulard, L. In Shock Compression of Condensed Matters2001, Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Atlanta, GA, June 2429, 2001; AIP Press: New York, 2002. (38) Jones, H. D. In Shock Compression of Condensed Matters2003, Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Portland, OR, July 20-25, 2003; AIP Press: New York, 2004. (39) Liu, H.; Zhao, J.; Ji, G.; Gong, Z.; Wei, D. Physica B 2006, 382, 334. (40) Crouzet, B. Private communication. (41) Marsh, S. P. LASL Shock Hugoniot Data; University of California Press: Berkely/Los Angeles, 1980. (42) Klebert, P. Etude expe´ rimentale et the´ orique de la transition choc de´ tonation dans les ex- plosifs homoge` nes. Ph.D. Thesis, Paris VI-Pierre et Marie Curie University, 1998. (43) Lysne, P. C.; Hardesty, R. D. J. Chem. Phys. 1973, 59, 6512. (44) Maillet, J.-B.; Soulard, L.; Stoltz, G. Euro. Phys. Lett. 2007, 78, 68001.
