Simulation of combustion-to-explosion transition processes in Octogen-based explosive

Introduction

During their lifetime cycle, products containing explosives may be exposed to various unanticipated emergency impacts (falling, fire action, bullet impact, etc.). The most frequently occurring emergency is fire. It is known that at explosive ignition, various outcomes can be possible: product containment failure and explosive burning out, convective combustion, low-velocity explosive transformation, or detonation [1].

The most severe consequences of accidents are associated with detonation. In this case, high-velocity fragments are produced, which may cause detonation in products located nearby. At the same time, in case of ignition of octogen-based explosives, a low-velocity explosive transformation scenario can be possible. In this case, larger fragments will be formed, but flying at a lower velocity than fragments produced by detonation. In view of this, it is vital to be able to predict product reaction to the impact of thermal field generated by fire.

The product reaction can be predicted using advanced finite-element analysis software applications, such as ANSYS Autodyn. The ANSYS Autodyn software features a fairly flexible combustion model Powder Burn, suitable for describing combustion processes in a wide range of pressures typical for transient conditions of explosive transformation (0.002...0.030 Mbar).

Powder Burn model structure

The model includes equations for the state of solids and gaseous products of decomposition, as well as for solid burnout rate [2]. A solid specimen is presumed to consist of particles of certain shape and corresponding characteristic size r₀.

The equation for the state of gaseous products of decomposition can be set in the form of the JWL equation. To describe the behaviour of a solid, the compaction equation can be used, for which it is necessary to set two tabular functions: the first one is the relationship between solid-phase pressure vs. density, and the second – relationship between sound velocity vs. density. The burnout rate is calculated by the formula

where F – burnout magnitude (varying from zero to unity);
t – time;
H (P_g) – tabular dependence of particle layerby-layer combustion rate H vs. gas-phase pressure P_g .

Here also, constant G characterises particle size. For example, for a spherical particle G = 3/r₀ (r₀ – radius), for a cubical one – G = 1 / r₀ (r₀ – rib length). Constants c characterise particle shape. For example, c = 2/3 for a spherical or cubical particle, c = 1/2 – for a cylindrical one.

Experimental data

Using the Powder Burn model, we shall try to reproduce the experiment results [3] with regard to the combustion-to-explosion transition in high-density octogen-containing LX-10 explosive samples confined in a solid assembly.

The experiment setup schematic is shown in Fig. 1. Five cylindrical samples of the LX-10 explosive (octogen – 95 %, Viton binder – 5 %) were placed in a steel tube with inside diameter of 44.93 mm, length of 203.2 mm, and wall thickness of 2.93 mm. One of the explosive samples contained an igniter.

After igniter actuation, an explosion occurred. The recorded velocities V of tube and fragments separation in the plane passing through the igniter and perpendicular to the assembly axis are shown in Fig. 2. The velocity increased for about 250 µs (see Fig. 2), which is a proof of non-detonative character of the explosion. A substantial increase in tube velocity, from 10 m/s to 400...700 m/s (before fragmentation), occurred during 30...41 µs. The values of fragments velocity was changing within a range of 800…2100 m/s.

Initial data for calculations

The purpose of the study was to describe the dynamics of steel tube separation based in the experiment [3], using the Powder Burn model. For that, the initial data are required, which were taken from papers [3], [4]. For example, the paper [4] provides parameter values for the JWL equation of state both for solid LX-10 explosive and for the gaseous products of LX-10 detonation (Table 1). The JWL equation appears as follows:

where P – pressure;
ρ – density;
E – internal energy;
A, B, R₁, R₂ – parameters;
ρ₀ – initial density of the LX-10 explosive.

For description of strength characteristics of the LX-10 explosive, a von Mises strength model with the shear modulus equal to 0.05 Mbar and yield strength Y₀ = 0.002 Mbar was selected.

The parameters for detonation products taken from Table 1 can be used right away, because in the Powder Burn model, the explosive decomposition products are also described by the JWL equation.

For solid behaviour description, the Powder Burn model employs the compaction equation. This equation was originally intended for modelling the behaviour of porous substances, such as loose sand or gunpowder of bulk density. It is known that substances in the bulk state become considerably compacted under load, and with load removed, their density drops, although not back to the initial state but to a certain value which is greater than the initial density. The compaction equation shall be adapted to address the problem under consideration, in which the initially high-density low-porosity explosive samples are used. In order to determine the compaction equation parameters, we shall use the JWL equation parameters for a solid explosive LX-10 (see Table 1).

For the compaction equation, it is necessary to set two tabular functions: function P(ρ), which defines the dependence of solid-phase pressure vs. density in the process of loading; function c (ρ^(A)), which describes unloading of earlier loaded substance down to a certain pressure and defines the dependence of sound velocity vs. density of the unloaded substance. To calculate the values of tabular functions included in the compaction equation, a curve for the JWL equation should be plotted for the case of blast-shock wave loading with the parameters taken from Table 1 (Fig. 3).

The values of tabular function P(ρ) can be select right away. The curve in Fig. 3 is obtained from the JWL equation for solid explosive, with the points corresponding to tabular function P(ρ). Let us take ten points belonging to the curve that corresponds to the JWL equation. The values of this tabular functions that defines loading of the explosive are given in Table 2.

To fill in tabular function c(ρ^(A)), we can proceed as shown in Fig. 3, where the curve corresponds to the JWL equation for a solid explosive, and the points correspond to tabular function P(ρ) (see Table 2). The line passing, for example, through points (ρ₄ ; P₄) and (ρ₅ ; P₅), will intersect the abscissa axis in point (ρ₄^(A); P = 0). The tangent of the slope angle of this straight line to the abscissa axis is equal to sound velocity square c₄ . Unloading of the explosive, initially loaded up to the state (ρ₅ ; P₅), will be occurring along this line. As a result, a pair of values (ρ₄^(A);c₄) is obtained. The rest of the pairs of values can be obtained in a similar way. Table 3 contains the values of tabular function с(ρ^(A)), defining explosive unloading.

For steel AerMet 100, the Grüneisen equation of state was used:

where Γ₀ = 1.84;
ρ₀ = 7.923 g/cm³;
c₀ = 0.4534 cm/µs;
S₁ = 1.5 [3].

We introduce the Steinberg – Guinan strength model as well:

where Y – dynamic yield strength;
– equivalent plastic strain;
G_сд – dynamic shear modulus;
T – temperature.

The parameters for steel AerMet 100 [3] are as follows. At temperature T₀ = 300 K and pressure P = 0, the initial shear modulus G₀ = 0.748 Mbar. The initial yield strength Y₀ = 0.01016 Mbar at temperature T₀ = 300 K and pressure P = 0. The maximum yield strength Y_max = 0.03 Mbar, β = 2, b = 0, n = 0.5, h = 0.

It is known that for quantitative description of the process of a relatively strong shock wave transformation into a detonation wave in high-density explosives, the Ignition and Growth model [4] is used, in which the burnout equation (1) can be written as follows:

Here, ρ – current density;
ρ₀ – initial explosive density;
I, G₁, G₂, a, b, c, d, e, g, x, y, z, F_{ig max}, F_G1max, F_G2min – parameters.

This model is represented in the ANSYS Autodyn and LS-DYNA software applications. For the LX-10 explosive, it is reasonable to use it in pressure ranges of 0.015…0.5 Mbar. To apply the Powder Burn model in this pressure range, it is necessary to adjust its parameters to comply with the Ignition and Growth model parameters, whose values are given in Table 4 [4].

In the Powder Burn model, combustion rate H(P_g) is set by the tabular method, therefore, for each value of pressure P_g magnitude it is possible to set such a value of combustion rate H, which will allow to describe on the average the behaviour of curve dF/dt, specified by the equation (2). For that, given a specific pressure value, we shall find the mean value of function dF/dt in the interval of burnout values F = 0–1. The mean value of function dF/dt is computed by the formula

Let us assume that in equation (1) G = 1, and c = 0.001. Then,

The value of parameter c = 0.001 corresponds, for example, to cubical particles burning on one of the six surfaces only. Coefficient G is taken equal to unity for the convenience of determining H(P_g). Its actual values for the LX-10 explosive lie within a range of 100…200 1/cm (explosive particle size r₀ ~ 1/G = 0.005...0.01 cm). Thus, values of table function   will contain information on parameter G . Now we can select the desired number of typical pressure values and calculate parameter  for each value selected (Table 5).

Then we select the values of function H(P_g) for pressures P_g < 0.01 Mbar. The layer-by-layer combustion rates, measured at P_g = 0.03 Mbar and P_g = 0.12 Mbar, are equal to V = 7∙10^-4 cm/µs and V = 0.03 cm/µs [5]. Multiplying the combustion rate values by 100, we have 0.07 and 3, respectively, which is close to the values of functions H (0.03 Mbar) = 0.077 µs^–1 and H (0.12 Mbar) = 2.792 µs^–1 (see Table 5).

In a formal way, it turns out that the equation (1) must have G ≈ 100. Such a value corresponds to explosive particle size r₀ ~ 1/G = 0.01 cm and to real sizes of octogen particles. Let us assume that G = 100 for pressures P_g < 0.01 Mbar as well. Then we multiply layerby-layer combustion rate values 2∙10^–5 cm/µs, 5.5∙10⁵ cm/µs, and 1∙10^–4 cm/µs [5] by 100 at the values of P_g equal to 0.0016 Mbar, 0.004 Mbar, and 0.007 Mbar. In so doing, we take the constant G = 1. The values of function H(P_g) for pressures P_g < 0.01 Mbar (Table 6) will be obtained.

The first two values of the tabular function are selected such that in the calculation, at the initial moments of time when the pressure is low, there is a weak reaction within the entire explosive. Then, as compression waves or shock waves pass through the explosive, the reaction rate will increase by orders of magnitude.

Computation results

The computational geometry is shown in Fig. 4. The grid size is 0.1 cm. Computations were conducted in the axisymmetric ALE formulation using the ANSYS Autodyn software (licence No. 774-2013-ША dated 01.10.2013). Ignition would start in several elements in which the initial pressure of 0.002 Mbar was applied for 30 µs. The impact of the igniter, as well as of the channels in the explosive samples in which wires for the igniter were laid, was not considered in the computations.

As a result of computation using function H(P_g) (see Tables 5, 6), and at G = 1, it was determined that tube separation velocity was increasing from 10 to 100 m/s during 16 µs and from 100 to 700 m/s during 15 µs, with transition to detonation occurring quite early. In the experiment, velocity was increasing from 10 to 100 m/s and from 100 m/s to 700 ms during 19 µs. If G = 0.7, then tube separation velocity increases from 10 to 100 m/s during 26 µs and from 100 to 700 m/s during 21 µs, with detonation occurring closer to the cover. If G = 0.5, then tube separation velocity increases from 10 to 100 m/s during 30 µs and from 100 to 700 m/s during 36 µs with no transition to detonation occurring.

Computations at G = 0.7 are the most consistent with the experiment. The results of this computation are shown in Fig. 5. A comparison between the computational and experimental data is given in Fig. 6.

The computational and experimental displacement of the tube wall is ~1 cm. The maximum computed tube velocity is ~2100 m/s, which corresponds to the maximum velocities of fragments recorded in the experiment. According to [3], destruction of the tube starts at 15 % radial strain. Under such deformation, tube velocity in the computation reached 500 m/s, which corresponds to experimental velocities of 400...700 m/s under which the fragmentation starts.

Conclusion

Using the ANSYS Autodyn finite-element analysis software and its integrated Powder Burn combustion model, we managed to satisfactorily describe the results of an experiment for simulating the combustion-to-exposition transition in LX-10 explosive samples enclosed in a solid assembly [3]. In the Powder Burn model, the burnout rate is selected by the tabular method, therefore a possibility appears to describe energy release for a wide range of pressures (P_g = 0.002...0.5 Mbar) where transition from combustion to detonation may occur.