New developments of the CARTE thermochemical code: I-parameter optimization

We present the calibration of the CARTE thermochemical code that allows to compute the properties of a wide variety of CHON explosives. We have developed an optimization procedure to obtain an accurate multicomponents EOS (fluid phase and condensed phase of carbon). We show here that the results of CARTE code are in good agreement with the specific data of molecular systems and we extensively compare our calculations with measured detonation properties for several explosives.


Introduction
Thermochemical methods have been used with considerable success to calculate the properties of the detonation products [1][2][3][4].The advantages are a very fast calculation of the properties (few seconds) and a direct evaluation of the Gibbs free energy.The thermochemical approach then becomes of great interest for choosing new molecules to be synthesized or for selecting an explosive for a given practical application.These methods can be very accurate provided that the different parameters are carefully chosen.The parameters are generally optimized to reproduce a very large variety of experimental (and ab initio) data.In this work, we present a new calibration of the CARTE thermochemical code [2].

Models
The CARTE code computes the equilibrium composition of chemical species by minimization of the Gibbs free energy of a system of detonation products.The system can be formed of one or possibly two fluid phases and of a condensed carbon phase.The fluid phase is composed of The potential energy is given by an EXP-6 potential.In the case of mixtures, an average potential is obtained through a mixing law.We apply an equation of state for the EXP-6 fluid to obtain the free energy of the system.In this version of the CARTE code, we use the following equation : -Buckingham-Keesom method [5] to treat dipolar molecules -Van der Waals one-fluid approximation proposed by Ree [6] to treat the mixture of chemical species -Statistical physics method MCRSR [7], a variational perturbation theory which is based on the inverse 12th-power potential as a reference system to calculate the free energy.
We have checked the accuracy of the treatment of fluid phase by comparison with molecular simulations results.The performance are satisfactory but suggest some new developments which are presented in an other work [8].
a e-mail: vincent-jp.dubois@cea.frThis is an Open Access article distributed under the terms of the Creative Commons Attribution-Noncommercial License 3.0, which permits unrestricted use, distribution, and reproduction in any noncommercial medium, provided the original work is properly cited.

EPJ Web of Conferences
The decomposition of a carbon-rich explosive produces an important amount of carbon, made of nanoparticles with different shapes, sizes, environments, . . .[9].In order to characterize the intrinsic properties of an explosive C x H y O z N w (of molar mass M), we use the oxygen balance (OB), which is the excess oxygen mass after the formation of CO 2 and H 2 O in g/kg of explosive: If OB is negative, the amount of oxygen is not sufficient to consume all the carbon into CO 2 and carbon nanoparticles are produced.To model the condensed carbon, the distribution of clusters is replaced by an effective cluster in equilibrium with fluid phase and treated by a modified 3 phases (graphite/diamond/liquid) Ree-van Thiel EOS [10]: -the EOS of each bulk phase is modified by the 0 K Birch potential's parameters (V 0 ,E 0 ) to take into account the size effect of clusters, -the EOS of liquid carbon is obtained by a scaling model from effective graphite/diamond solids, -we use a Boltzmann law to mix the carbon phases.
To take into account the evolution of the carbon EOS, we postulate that the size of the effective cluster depends on the quantity of carbon in excess, and thus the oxygen balance.The energy of formation (E 0 ) and the molar volume (V 0 ) of the 0 K isotherm also evolve with oxygen balance of the explosives.

Parameters
The parameters have to be optimized in order for the thermochemical code to reproduce some reference data.We minimize error functions representing the difference between reference data and calculated properties by a quadratic normalized function.We developed an original method which performs the minimization in two steps.An exploration step (Monte Carlo scheme) allows to visualize the surface parameters and to extract particular points which act as initial points of the minimization step (Nelder-Mead method).A similar procedure has been recently applied to fluid force fields for molecular simulations [11].In the CARTE code, the parameters can be separated in three groups: -the parameters of EXP-6 potentials describing the interactions between like molecules in fluid phase (pure fluids) -the parameters describing the interactions between unlike molecules (cross potentials).Indeed, the cross potential between two unlike molecules is usually given by the Lorentz-Berthelot rules.Nevertheless, the accuracy of these rules decreases with increasing pressure and temperature.We thus adjusted these rules by three parameters (m i j , k i j , l i j ): -the parameters related to the carbon EOS.
Our approach is based on a sequential procedure in three steps.We follow a logical order for the fluid phase parameters for which we have specific reference data: first the pure fluids and secondly the cross potentials.The lack of specific reference data for the carbon EOS implies the use of global data (CJ detonation properties of carbon-rich explosives) and puts back the optimization of the carbon EOS parameters at the last step, assuming that our fluid phase EOS is correct.

00031-p.2
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter

Pure fluids
We built our database with all the available experimental data in the thermodynamic range (0.05-100 GPa/400-6000 K).We have thus applied our minimization routine with these large databases.For prior species (CO 2 , H 2 O, N 2 ), we use nearly 150 data per molecule with direct or derivative properties along isothermal, isochoric or isobaric curves.The potentials obtained are able to reproduce all the data in the thermodynamic domain (0.05-100 GPa/400-6000 K) with a standard deviation less than 2%.

Cross potentials for mixtures
These parameters are obtained by matching reference data.Nevertheless, this calibration is not easy due to the lack of specific results on binary mixtures in the domain of detonation products.The parameters are generally obtained by matching the detonation velocity of some explosives [2], but this global optimization can hide some balancing errors between the different contributions.Ab initio data appear as an alternative way to supply us with specific binary data.We used ab initio results [12] on binary mixtures involving (CO 2 , H 2 O and N 2 ) and some other experimental specific results to calibrate the cross potentials between the prior species.The other cross parameters are taken equal to unity.
The cross potentials H 2 O-N 2 and CO 2 -N 2 remain close to the Lorentz-Berthelot rules.On the contrary, the mixture H 2 O-CO 2 at high pressure and temperature strongly deviates from idealy (decrease of volume) with the formation of new species CO 3 H 2 [12].To model the properties of mixture, we have introduced a two-species model where the water is represented by two species in chemical equilibrium : non-associated water (standard water H 2 O) and associated water (H 2 O i ).The agreement between DFT results and thermochemistry data is shown in Fig. 1.The decrease of volume observed at 3000 K and 20 Gpa in the DFT results is fairly reproduced by our empirical model.This behaviour is obtained by an increase of H 2 O i species with the pressure and temperature.

Carbon EOS
The EOS bulk of carbon is modified by changing the energy of formation (E 0 ) and the molar volume (V 0 ) of the 0 K isotherm which evolve with oxygen balance of the explosives.The optimization 00031-p.3procedure of (E 0 ,V 0 ) values is decomposed in three steps:

EPJ Web of Conferences
-We picked out a very large database of CJ points (150 values of D CJ and 50 values of P CJ ) and separated these data by class of oxygen balance (OB) between −800 and +70.-For each OB class, we fitted (E 0 ,V 0 ) by matching the CJ data.
-At last, we adjusted analytical functions to describe the evolution of (E 0 ,V 0 ) with OB.
With this new calibration, the CARTE code is able to reproduce the detonation properties of a very large panel of high explosives.As an illustration of the efficiency of our code, we reproduced with a very good precision the evolution of experimental CJ properties (velocity and pressure) with the initial density for carbon-poor explosive (like PETN, see Fig. 2), as well as for carbon-middle explosives (like HMX/RDX, see Fig. 3), as well as for carbon-rich explosives (like HNS/TNT, see Fig. 4).In particular, we are able to reproduce the smooth transitions in condensed carbon phase as shock conditions vary.

Conclusion
We have improved the thermochemical code CARTE dedicated to the computation of thermodynamic properties of CHON systems with high accuracy at high pressures and temperatures.CARTE code is proven to reproduce both specific experimental or ab initio data (shock or static compression of 00031-p.The CARTE code is able to calculate the performances of a wide variety of explosives, and can help in choosing new explosives to be synthesized or in implementing an explosive for a given practical application.Whatever the performances of the code, several new developments are in progress: -Bayesian theory to evaluate the uncertainties on the parameters, -Improvement of fluid phase EOS, presented in paper [8], -Development of a new equation of state of Carbon to better describe the ultra dispersed particles of carbon in detonation products.

2 Fig. 1 .
Fig. 1.Top: isotherm at 3000 K for the mixture H 2 O-CO 2 obtained from CARTE code (red line) and compared to DFT results[12] (black points with error bars).Bottom: proportion of associated water (H 2 O i ) along the 3000 K isotherm.

Fig. 2 .Fig. 3 .
Fig. 2. Evolution of detonation velocity (top) and pressure (bottom) as a function of the initial density for PETN (carbon-poor).The experimental data are represented by red crosses and solid black lines represent results calculated with CARTE code.

4
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter Evolution of detonation velocity (top) and pressure (bottom) as a function of the initial density for HNS/TNT (carbon-rich).The experimental data are represented by red crosses and solid black lines represent results calculated with CARTE code.pure fluids for example) and global experimental results (such as CJ state properties of explosives).