Molecular dynamics simulation of the burning front propagation in PETN

One of the models of detonation development in condensed explosives under shock loading is the concept of "hot spots." According to this model, the reaction initially starts at various defects and inhomogeneities, where energy is localized during shock wave propagation. In such a region the reaction may start and the heat flux sufficient for the ignition of the adjacent layers of matter may be formed. If the reaction propagates fast enough, the merging of the burning fronts from several hot spots may lead to detonation. So there is an interest in determining the burning propagation rate from the hot spot in various conditions. In this work we investigate the propagation of plane burning front from initially heated layer in PETN single crystal using molecular dynamics method with the reactive force field (ReaxFF). The burning rate depends on the direction in crystal. The kinetics of chemical transformations is considered. The dependence of the burning front propagation rate along [100] direction on the external pressure in the pressure range from normal to 30 GPa is calculated, it is shown that it grows linearly in the considered range from 50 m/s to 320 m/s. The results are compared with the data from experiments and quantum chemical calculations.


Introduction
The hot spot model of condensed explosive initiation assumes that initial heating and initial chemical decomposition reactions occur in defect regions, like voids or plastic-explosive boundaries, of the explosive structure. Estimations show that the heating due to shear deformations in such regions may reach more than 1000 K in shock waves of about 1 GPa intensity [1], forming the "hot spot". Such temperatures may start the sequence of chemical transformations, and the heat yield from these transformations can be high enough to give rise to the decomposition of the adjacent layers of explosive material. Thus, the reaction may propagate from the regions initially heated by a shock wave. If the hot spot concentration and the reaction propagation rate are high enough, growing hot spots may lead to the formation of a detonation front. The hot spot concentration for a particular explosive may be determined experimentally by studying the microstructure of the material, but the evaluation of the reaction propagation rate is more complicated. While there are novel experimental techniques that promise the precise numerical data on reaction kinetics on picosecond time scales [2], the decomposition mechanisms in various conditions are still in question. In the present work, we apply the molecular dynamics (MD) method with the ReaxFF interatomic potential [3] to study the reaction propagation in PETN and the kinetics of chemical events in a heated region.

Simulation model
First part of the work is the calculation of the reaction propagation rate in PETN single crystal. We consider a unidimensional propagation of the plane reaction front. To minimize size effects, the simulation cell was made large along reaction propagation direction [100] (256 unit cells). Along transversal directions the cell was 4 × 4 unit cells. One unit cell contains 58 atoms, so the total number of atoms in these calculations was 237568. First, simulation cell was posed in Nose-Hoover thermo-and barostat to obtain normal temperature and target pressure. Then we chose a thin layer (about 5 unit cell wide) in the middle part of the cell and set the temperature there about 2000 K. It is considered as a model of the formed hot spot. Its evolution is then simulated in NVE ensemble.
For the second part of the work (kinetics of decomposition) we considered 4 × 4 × 4 unit cells (3712 atoms) and 8 × 8 × 8 unit cells (29696 atoms) systems. Similarly, normal temperature and pressure in the cell were first obtained with Nose-Hoover technique, then higher temperature from 500 to 2100 K was set in the whole cell. After setting the high temperature MD simulations were performed in NVE ensemble. So there was no energy leakage and temperature could grow as a result of chemical reactions.
ReaxFF is a reactive interatomic potential. It is based on the concept of chemical bond order between atoms, that depends on the instantaneous chemical environment of each atom. Initially, bond orders are evaluated from the interatomic distances, then they are corrected to give consistent valence values for each atom. On the basis of the information on bond orders, energy and forces are calculated in the way similar to classic potentials for bonded systems: potential energy is decomposed to several terms corresponding to bond lengths, angles, dihedrals and so on. As total bond order for each atom may change during simulation, charges on atoms do not remain constant and are equilibrated on each simulation step, too. There is screened Coulomb term in potential energy to take charge interaction into account, and there is also the LJ term to describe van der Waals interactions. ReaxFF is fitted to the energies of large number of atomic configurations characteristic for a given class of materials. Here we use the ReaxFF parametrization for PETN in the form used in the work [4]. Calculations were done in LAMMPS molecular dynamics package [5] with USER-REAXC module [6] as a ReaxFF implementation.
During simulation, we periodically printed out atomic coordinates and bond orders between each atom and its neighbors. We then determined bonded clusters and their positions in the simulation cell at the considered timesteps.

Reaction propagation rate
The character of the reaction propagation from heated region depends on the energy barriers and heat yields of the reactions, thermal conductivity of the material and its response to the pressure applied. To determine the reaction propagation rate (RPR), we compared several criteria based on the constant level of concentration of chemical species that take part in the reactions as well as on constant temperature value. All the criteria gave very similar results. There may be different mechanisms of the decomposition in different conditions, so we primarily used PETN concentration instead of concentration of individual products to calculate RPR. We divided simulation cell to 128 layers (each layer 2 unit cells wide) and considered that when PETN concentration in a layer equals half of the initial one, it corresponds to the reaction front position.
On figure 1 the PETN concentration in the middle part of the simulation cell is shown. External pressure in this calculation is 10 GPa. Initially heated layer corresponds to x values from 1200 to 1250Å. On the z axis there is the number of PETN molecules in each layer (64 at simulation start). Levels of constant PETN concentration are shown on x-t plane. We place the reaction front in the coordinate where PETN concentration equals 32, i.e. half of the initial one, at a particular time step. Values of RPR were checked for convergence with respect to We also performed one calculation of the RPR along [001] direction at 1 GPa, and the RPR proved out to be about 2 times lower than along [100] direction. It may be related to geometric parameters of PETN unit cell, that is shorter in [001] direction than in [100] one. A layer normal to [001] contains more molecules and thus requires more energy for reactions to start than the one normal to [100] that has the same thickness. There may also be differences in thermal conductance along different axes. Finally, as PETN molecules have complex geometry, their orientation with respect to the reaction front will affect rates of reactions. Additional simulations are needed to see if the relation between RPR along different axes will remain under high pressures.
The molecular dynamics results on the reaction propagation rate is of the same order with diamond anvil cell experiments in the high pressure region. However, the experimental values were sometimes several times higher than calculated ones. It may be partially explained by the fact that there was a compressed disperse PETN in experiments that contained a lot of pores, boundaries and other defects. Another feature of the MD model is that it does not explicitly take into account electronic degrees of freedom, so there is no electron thermal conductivity that might also facilitate reaction propagation at high pressures.

Kinetics of the decomposition in a heated region
PETN decomposition was studied previously by DFT [9,10] and DFTB [11] approaches. In the report of Wu et al. [11] a scheme of reactions was proposed leading to the formation of water   Concentration of water per PETN molecule for different initial temperatures versus time. The dependencies on figure 3 correspond to the process that looks like relatively long induction time with slow product accumulation, followed by fast self-accelerating stage of complete transformation to final products. As this second stage is faster, we describe the process with only one parameter τ that is essentially the mentioned induction time. As product formation is not a first-order reaction with constant rate for each initial temperature, we do not calculate real activation energy of a particular chemical process using figure 4, but rather the effective activation energy of the transformation as a whole. It describes change of induction time instead of reaction rate constant, and reflects changes caused by all factors including temperature growth during reaction and possible catalytic effect of initial products that stay in the system. The slope of the linear approximation in figure 4 is 5.8, so the effective activation energy for water formation in PETN decomposition reaction chain appears to be 48 kJ/mole. It is not directly connected to any of intramolecular bond energies and is significantly lower than the energy of weakest bond O-NO 2 in isolated PETN molecule, that equals 150-180 kJ/mole [10,12]. That weakest bond energy value was also reproduced in our ReaxFF calculation for one PETN molecule. Difference between any of bond energies and the calculated effective activation energy is natural as the latter describes none of elementary acts but the process going under varying temperature and chemical composition).

Conclusions
We determined the reaction propagation rate in PETN single crystal using molecular dynamics with ReaxFF interatomic potential. The RPR along [100] direction grows linearly with external pressure from 50 to 320 m/s in pressure range from 3 to 30 GPa. These values are in good agreement with known experimental data. We also considered the kinetics of thermal decomposition in uniformly heated bulk PETN in NVE ensemble. The process goes in two stages. The first one is slow product formation during relatively long induction time τ , that is succeeded by the second stage of fast self-accelerating reaction. Dependence of τ on initial temperature in the system is Arrhenius-like with the effective activation energy of 48 kJ/mole, with τ values in range of 10 −10 -10 −11 s for initial temperatures from 800 to 2000 K.