Numerical simulation for determining detonation parameters of explosive substances using EXPLO5 thermo-chemical prediction software

The article presents the theoretical and practical research results regarding the development of the IT and methodological infrastructure for assessing the parameters which define explosive substances. The software us intended for estimating, based on the chemical composition of substances, parameters which define such substances (detonation velocity, detonation energy, detonation pressure etc.), as well as for determining the probability of explosion. EXPLO5 is a specialized software which allows to solve thermo-dynamic equations between different reaction products in order to determine the balance composition, having incorporated a filtering algorithm for estimating Jones-Wilkins-Lee (JWL) state equations coefficients which are used for calculating the detonation energy of explosives. For the proper performance of pilot tests of numerical simulations upon certain representative types of explosives (civil use explosives and double-use explosives: civil and military) there have been taken into account the special characteristics of EXPLO5 provided by the EOS BKW (BKW standard with set of BKWN constants and BKW EOS modified with the set of constants and volumes BKWN-M) and EXP-6 EOS function which is based more on theoretical concepts.


Introduction
EXPLO5 is thermochemical computer code that predicts detonation (e.g. detonation velocity, pressure, energy, heat end temperature etc.) and combustion (e.g. specific impulse, force, pressure etc.) performance of energetic materials. The calculation of detonation parameters is based on the chemical equilibrium steady state model of detonation. The equilibrium composition of detonation and combustion products is calculated by applying modified White, Johnson, and Dantzig's free energy minimisation technique 1.
The program uses the Becker-Kistiakowsky-Wilson (BKW) and Exp-6 EOS for gaseous detonation products, the ideal gas and virial equations of state of gaseous combustion products, and the Murnaghan equation of states for condensed products. EXPLO5 is designed so that enables calculation of chemical equilibrium composition and thermodynamic parameters of state along the shock adiabat of detonation products, the C-J state and the detonation parameters at the C-J state, as well as the parameters of state along the expansion isentrope. The program has non-linear curve fitting program built in to fit relative volume-pressure data along the expansion isentrope according to the Jones-Wilkins-Lee (JWL) model, enabling the calculation of the detonation energy available for performing mechanical work. EXPLO5 is thermochemical computer code that can predict the performance of ideal high explosives and various explosive formulations, propellants, and pyrotechnic mixtures. It predicts the performance under constant volume combustion, constant pressure combustion, and ideal detonation 2.

Material and method
EXPLO5 has possibility to calculate detonation parameters of an explosive, as well as combustion parameters under isobaric and isochoric conditions 3,4,5,6.

Detonation
The detonation is a process of layer-by-layer, supersonic propagation of chemical reactions through an explosive. According to the generally accepted Zeldovich-von Neumann-Doering (ZND) model of detonation, chemical reactions occur at a definite rate in the chemical reaction zone, under the action of a shock wave (Fig.1). Under the influence of the dynamic action of the shock wave, a thin layer of the explosive is compressed from the initial specific volume V 0 (V 0 =1/ 0 ,  0 is initial density of explosive) to the volume V 1 in accordance with the shock (or Hugoniot) adiabat for a given explosive (Fig.2). The equation defines the relationship between the density (or volume) and the pressure during the shock compression of the explosive. As a consequence of shock compression, the increase of the pressure to the value p 1 occurs, resulting in a significant temperature increase in the compressed explosive layer, where consequently the initiation of the chemical reactions takes place. When the chemical reactions are at their end, the volume and the pressure of the reaction products reach the V 2 and p 2 values. This state corresponds to the point lying on the shock adiabat for detonation products. From that point the products expand isentropically (Taylor wave) into the surrounding medium. According to the steadystate model of detonation, the points (V 0 ,p 0 ), (V 1 ,p 1 ) and (V 2 ,p 2 ) lie on one straight line. This line is called the Rayleigh line. The slope of the Rayleigh line is determined by detonation velocity of a given explosive. According to the Chapman and Jouguet hypothesis, the Rayleigh line is a tangent to the shock adiabat of the detonation products at the point that corresponds to the end of the chemical reactions. That point is assigned as the Chapman-Jouguet point (C-J point).

Non-ideal detonation model
The detonation properties of "ideal explosives" can be accurately modeled by the ChapmanJouguet (C-J) detonation theory. The C-J theory assumes instantaneous chemical reactions and thermodynamic equilibrium of the detonation products. However, in the case of commercial intermolecular explosives, polymer bonded and metalized high explosives, the C-J theory is not able to reproduce accurately experimental behaviour. The reason for this is long reaction time and a wide reaction zone in such explosives (microsecond timescale, comparing to nanosecond scale in the case of ideal explosives), i.e. huge deviation from the C-J theory assumption on instantaneous thermodynamic equilibrium. The Zeldovich-von Neumann-Doering (ZND) detonation model takes into account existence of a finite width reaction zone, bounded by the von Neumann spike and the C-J point. The ZND theory considers a plane detonation wave, and it predicts the same detonation velocities as the C-J theory. In other words, it also fails to predict accurately detonation velocity of explosives which have highly curved detonation front and exhibit strong dependence of detonation properties on charge diameter (so-called "non-ideal" explosives). To model detonation parameters of "non-ideal" explosives we applied Wood and Kirkwood slightly-divergent detonation (WK) theory. The theory considers cylindrical explosive charge of infinite length and solves Euler hydrodynamic flow equations along the central streamline of the cylinder. It treats radial expansion as a first order perturbation to a perfect one-dimensional flow along the streamline and predicts detonation velocity as a function of rate of chemical reactions and rate of radial expansion. Short theoretical background of the theory is given bellow.
In the case of non-ideal detonation the shock front is rather curved and the flow of the reaction products significantly diverges. The result of this is that reactions are never complete in the detonation driving zone (DZZ) -between the shock front and the sonic line. The flow in DDZ is subsonic so that perturbation such are compression waves and rarefaction waves, which moves at the local sound speed, can contribute to support the shock front. On the other hand, the flow behind the sonic line is supersonic, so that perturbations can never catch up with the DDZ and contribute to the shock front (Fig.3).   Fig. 3. Two-dimensional representation of non-ideal detonation.
The reactions still continue in the Taylor rarefaction wave, between the sonic line and end of reaction zone, and they cannot contribute to the shock front but they can significantly contribute to the blast.

Isobaric combustion
Combustion of an energetic material is irreversible process in which it transforms into gaseous (mainly) and condensed products. The combustion is followed by the heat liberation. The combustion process of C C H H N N O O types of energetic materials can be presented by the following general equation: EXPLO5 enables calculation of composition of combustion products under different combustion conditions (isobaric or isochoric conditions), as well as thermodynamic parameters of combustion products. The calculation of thermodynamic parameters of combustion products under constant pressure conditions (isobaric combustion.) is based on the assumptions that the combustion of an energetic material proceeds without heat losses to the surrounding (i.e. adiabatically), and that in the combustion products the state of chemical equilibrium establishes. The calculation is performed as follows: (a) The system of equations mathematically describing the state of equilibrium in the combustion products is formed on the basis of the free energy minimization method. The system of equations is solved applying modified Newton's method; (b) State of gaseous combustion products is described by the ideal gas equation of state; (c) The thermodynamic functions of state of combustion products in their standard state are calculated from the enthalpy; (d) Condensed combustion products (if any) are considered as incompressible.

Isochoric combustion
The combustion of an energetic material is irreversible process in which it transforms into gaseous and solid products. The calculation of thermodynamic properties of combustion products under constant volume conditions (isochoric combustion) is based on the assumptions that the combustion of an energetic material proceeds without heat losses to the surrounding (i.e. adiabatically), and that in the combustion products the state of chemical equilibrium establishes. The calculation is performed as follows: (a) The system of equations mathematically describing the state of equilibrium in the combustion products is formed on the basis of the free energy minimization method. The system of equations is solved applying modified Newton's steepest descent method; (b) The state of gaseous combustion products is described by the virial equation of state, or by the ideal gas equation of gas (optionally); (c) Condensed combustion products (if any) are considered as incompressible; (d) The thermodynamic functions of combustion gases, as real gases, are derived applying thermodynamics laws and the virial equation of state; (e) The thermodynamic functions of combustion products in the standard state are derived from the enthalpy.

Results and discussion
The main results of the substances testing are provided in Fig. 4-9:    -The users interested in optimisation of parameters in equations of state of gas products can calculate shock Hugoniots of individual products (i.e. shock compression from initial state to specified maximum pressure) and compare calculated results (usually p-V and p-T data) with experimentally determined. In this way one can adjust parameters in EOS to fit the best experimental p-V data (Fig.4).
-Implemented the "golden-section search", technique to find the minimum of the merit function Y(D). The first step in the minimisation is to make an initial scan through a range of detonation velocities (maximum and minimum detonation velocities and step are specified by user) to check whether a minimum (or more minima) exists, and to bracket the minimum (Fig.5a). Once the minimum is bracketed, the minimiser search for the exact location of minimum and detonation velocity corresponding to that minimum Y(D), Fig.5b.
-The progress of calculation along the shock adiabat of HMX (=1.906 g/cm 3 ) using Modified BKW, can be followed from the Progress tab which shows the shock adiabat graph (p(V)-V and D(V)-V, and the flow chart (Fig.6).
-Program calculates composition of detonation products and parameters of state (pressure, thermodynamic functions, detonation velocity etc.) along the entire shock adiabat of detonation products, starting from the initial density of explosive to a final maximum pressure on the shock adiabat and increasing density for an arbitrary chosen density increase step; The calculation along the shock adiabat terminates when detonation pressure becomes higher than the maximum pressure chosen (default value of maximum pressure is 1.4p CJ , where p CJ is estimated pressure at the C-J point; p CJ =10.5 0 2 GPa) or when density becomes larger than 2 0 . Then the program fits p- data in the vicinity of minimum value of D to fourth degree polynomial, and then calculates minimum of so-obtained function D=f() applying the "golden section search" method. The value of density corresponding to D min is density at the C-J point ( CJ ). Finally, the program calculates parameters of state at the  CJ , i.e. at the C-J point (Fig.7 a,b).
-When the Isobaric combustion of Nitrocellulose (12.5%N) calculated using Ideal gas EOS, running mode is selected then the Progress of calculation shown in Fig.8. The progress of calculation can be followed in two ways: graphically, following temperature convergence (i.e. (H T -H 298 )-Q p in combustion chamber, or (S T -S ch ) in the nozzle throat and the nozzle exit, and by following temporary composition of combustion products.
-In the same, the Isochoric combustion of Nitrocellulose (12.5%N) calculated using Virial EOS running mode is selected then the Progress of calculation shown in Fig.9.

Conclusions
Evaluation / verification of ballistic / thermodynamic / physico-chemical parameters can be done either by testing explosives in specialized laboratories or by using computer applications. In this respect, EXPLO5's specialized thermochemical prediction software presents a large database and provides the ability to predict the performance of disruptive (ideal and non-ideal) explosives, propellants and pyrotechnic mixtures based on chemical formula, heat of formation and of density, being useful tools in the synthesis, formulation and numerical modelling of energetic materials. EXPLO5 can calculate the chemical composition of equilibrium and product-specific thermodynamic parameters at a given state V, p, T applying the techniques of minimizing free energy. These data, together with the Chapman-Jouguet detonation theory, allow the detonation parameters to be calculated, such as detonation velocity, explosion pressure, explosion energy, etc. Thus, from the equilibrium composition and thermodynamic parameters of the state over the isentropic expansion, EXPLO5 can determine the coefficients in the Jones-Wilkins-Lee state equation (JWL) through a built-in JWL installation subroutine and available energy for mechanical work. The paper highlights a series of thermodynamic parameters, using established thermodynamic functions (BKW EOS, Modified BKW EOS and EXP-6), which were obtained by computerized numerical simulation on several types of explosives (ANFO, HMX, RDX and Nitrocellulose), as well as CO2 (to determine the Hugoniot Shock).
This paper was developed within the NUCLEU -Programme, carried out with the support of ANCSI, Project no. PN-19 21 02 02.