Parameter Determination of an

Ignition and Growth Model from

Emulsion Explosive Tests
C Yi1, U Nyberg2 and D Johansson3

ABSTRACT
Emulsion explosives are a common industrial explosive and have a non-ideal detonation
behaviour. The detonation performance for a given product changes with the charge diameter,
ground conditions, connement and density. The burning process of an emulsion explosive has
been modelled with an ignition and growth model in the LS-DYNA code. The parameters in the
burning rate function are calibrated by the detonation velocities and the detonation front curvature
radii from emulsion explosive tests. A Perl program was developed for an efficient estimation of
the parameters, and the results show that the calibrated parameters can predict both detonation
velocities and the detonation front curvatures.

INTRODUCTION
Emulsion explosives are widely used by the mining industry Much effort has been made to investigate the non-
due to their low manufacturing cost, resistance to water, ideality of explosives (Kennedy, 1995; Takahashi et al, 1999;
ease of handling and excellent safety features (low risk of Hamashima, Kato and Itoh, 2004; Braithwaite, Cunningham
accidental ignition) when being transported and stored. It and Sharpe, 2006; Watt et al, 2012). In order to simulate the
is important to predict the performance of these explosives non-ideality, the reactive flow and chemical kinetics of the
in a confined environment in order to optimise the blasting burning process should be included in the model. Several
results. It is well recognised that ammonium nitrate (AN)- non-ideal detonation models and/or theories have been
based explosives such as emulsion explosives exhibit strong, proposed to describe the detonation process of explosives
non-ideal detonation behaviour since their performance is (Bdzil, 1981; Bdzil, Aslam and Short, 2002; Bdzil and
influenced by charge diameter, density and confinement Stewart, 2007; Kirby and Leiper, 1985; Esen, 2008). Areview
(Kirby and Leiper, 1985; Braithwaite et al, 1990; Persson, of existing non-ideal detonation models was given by Esen
Holmberg and Lee, 1993; Forbes et al, 1992; Forbes and (2008). The popular Ignition and Growth (I&G) reactive
Lemar, 1998; Esen, 2004). The oxidiser and fuel in AN-based model (Lee and Tarver, 1980) is often used to study the
explosives are located on different molecules, leading to non- non-ideality of explosives (Souers, Garza and Vitello, 2002;
ideal detonation behaviour over a large range of charge radii Lu and Kennedy, 2003; Price and Ghee, 2009) and has
(Schoch et al, 2013). This means that in this range, the velocity been implemented in several hydrodynamic codes such as
of detonation (VoD) and the detonation front curvature LS-DYNA.
strongly depend on the confiner material, density and charge The Swedish Rock Engineering Research Foundation
diameter and cannot be accurately predicted by simple (SveBeFo) carried out a series of VoD and detonation front
detonation models (eg Chapman-Jouguet (CJ) theory). curvature measurements to investigate the properties
The conventional theories of ideal detonation, for example of emulsion explosives (Nyberg, Deng and Chen, 1994;
the CJ theory and the Zeldovich-Neumann-Doering (ZND) Nyberg and Deng, 1994). Numerical simulations for these
theory, cannot be used to describe the non-ideality of an experiments were carried out with a self-developed reaction
explosive. They describe a plane detonation wave of an rate equation, which was implemented in the LS-DYNA
instantaneous reaction (CJ theory) and a finite reaction zone code (Deng, Nie and Chen, 1995; Deng et al, 1999). In this
(ZND theory). The programmed burn model implemented paper, the reaction rate parameters for the I&G reactive
in some hydrodynamic codes such as LS-DYNA (Livermore model for emulsion explosives in light confinement were
Software Technology Corporation, 2014) is based on the CJ calibrated from the detonation velocities and the curvature
theory. With a constant input VoD value and CJ pressure, the radii of detonation fronts at different charge diameters
model simulates an idealised detonation under conditions that were measured earlier (Nyberg, Deng and Chen,
where the Jones-Wilkins-Lee (JWL) parameters have been 1994). Because there are several unknown parameters in
calibrated (eg by cylinder tests). the reaction rate function of the I&G reactive model, these

1. Associate Senior Lecturer, Swebrec at Lule University of Technology, Lule SE-971 87, Sweden. Email: [email protected]
2. Senior Research Engineer, Swebrec at Lule University of Technology, Lule SE-971 87, Sweden. Email: [email protected]
3. Senior Lecturer, Swebrec at Lule University of Technology, Lule SE-971 87, Sweden. Email: [email protected]

11TH INTERNATIONAL SYMPOSIUM ON ROCK FRAGMENTATION BY BLASTING / SYDNEY, NSW, 2426 AUGUST 2015 585
C YI, U NYBERG AND D JOHANSSON

unknown parameters need to be adjusted iteratively to fit the as it is compressed by the leading shock wave creating
experimental results. The conventional method (eg Deng, heated volumes (hot spots) as the voids in the material
Nie and Chen, 1995; Deng et al, 1999) to get a set of suitable collapse. The fraction of explosive ignited is approximately
parameters is often time-consuming. In this paper, a code equal to the original void volume. The second reaction term
written in the Perl programming language (Christiansen, models the rapid formation of the major reaction product
Wall and Orwant, 2012) was developed to run the LS-DYNA gases (CO2, N2, H2O, CO, etc) in highly vibrationally excited
code, adjust parameters and get the results automatically, states and their subsequent expansion and equilibration. The
which makes the parameter calibration more efficient. third term is used to describe the relatively slow diffusion-
controlled formation of the solid carbon particles in the form
IGNITION AND GROWTH REACTIVE MODEL IN of diamond, graphite or amorphous carbon (Tarver and
McGuire, 2002).
LS-DYNA
The parameters of two JWL EoS for the reacted and
The I&G reactive model in LS-DYNA originated with the
unreacted explosive can be determined by experiments.
work of Lee and Tarver (1980), and was refined in a number
There are several unknown parameters in Equation 3. In
of subsequent articles by Tarver and collaborators (Tarver,
order to calibrate these parameters, iterative simulations need
Hallquist and Erickson, 1985; Tarver, Kury and Breithaupt,
to be conducted using LS-DYNA by adjusting the parameters
1997). Although the model was initially built to study the
and repeating until a best fit to the experimental results is
sensitivity of ideal high explosives, it essentially simulates
obtained (Aldis, Quirk and Breithaupt, 1991).
the non-idealities of explosives. The I&G reactive model
contains two equations of state (EoS): one for the unreacted
explosive and the other for its reaction products. Both CALIBRATION METHODOLOGY
unreacted and product EoS are of the JWL form: A high resolution is required to capture the reaction of
explosives in the reaction zone, which is usually of the
order of millimetres. The elements of the model therefore
P = Ae-R1V + Be-R2V + ~Cv T V (1)
need to be small. An axi-symmetrical model was built and
discretised with rectangular elements. The element size was
where:
0.02 cm, which is sufficient for resolving the reaction details
P is pressure
of commercial explosives according to Deng, Nie and Chen
V is relative volume (1995). After the LS-DYNA calculation, the coordinates of the
T is temperature centroids of the elements with the maximum pressure along
is the Grneisen coefficient the element rows (in the detonation front direction) at a certain
Cv is the average heat capacity time were output. The fitted curve for these coordinates can
A, B, R1 and R2 are constants be treated as the detonation front at the selected time (Bergh,
The mixing rule is: personal communication, 2013). The loci of the detonation
fronts at two different times were output (see Figure 1) and
V = ]1 - Fg Ve + FVp (2) the distance (S in Figure 1) between the two sets of vertices
were calculated. Since the time difference is known, the VoD
can then be estimated.
where:
Ve is the unreacted explosive relative volume The detonation front can be fitted to several forms, such
as circles, second order polynomials, ellipses and the
Vp is the reacted explosive relative volume
natural logarithm of a Bessel function (Kennedy, 1998). In
F is the fraction reacted, where complete reaction means
this paper, the detonation front is fitted with an ellipse of the
that F = 1 and no reaction means that F = 0
form:
This parameter is calculated by the I&G reactive rate
equation. The reaction rate law has three terms and each
R2 = aZ2 + bZ + c (4)
is active for different values of the reacted fraction. This
three-term rate law describes the three stages of reaction
where:
generally observed in the shock initiation and detonation of
heterogeneous solid explosives: R and Z are the radial and axial coordinates respectively
a, b and c are the fitting constants
dF I ]1 - Fgb _ t t - 1 - a ix G ]1 - FgcF d p y G ]1 - Fge F f p z (3)
= + 1 + 2
dt 0

where:
F is the fraction reacted
t is time
0 is initial density
is the current density
p is the pressure
I, b, a, x, G1, c, d, y, G2, e, f and z are constants for the ignition
term and the two growth terms
Three more constants are added to the model: Fmxig, Fmxgr
and Fmngr, which limit the contributions of the three terms
respectively: a maximum reacted fraction Fmxig for the first
term, a maximum fraction Fmxgr for the second term and a
minimum fraction Fmngr for the last term. For detonation, FIG 1 Scheme for acquiring velocity of detonation.
the first reaction term represents the ignition of the explosive (A) Detonation front at T2; (B) detonation front at T1.

586 11TH INTERNATIONAL SYMPOSIUM ON ROCK FRAGMENTATION BY BLASTING / SYDNEY, NSW, 2426 AUGUST 2015
 PARAMETER DETERMINATION OF AN IGNITION AND GROWTH MODEL FROM EMULSION EXPLOSIVE TESTS

The curvature radius of the detonation front is a function of TABLE 1

R (or Z) and is calculated by the following equation: Percentage of contents and density of the emulsion explosive.
Rm AN (%) EM (%) 7514oil (%) Water (%) MB (%) 0 (g/cm3)
k ]Rg = (5)
^1 + Rl2h3 2
77.49 0.748 3.926 14.76 3.075 1.165
As mentioned earlier, the calibration process of reaction AN ammonium nitrate; EM Anfomul S6 emulsifier from Croda Oleochmicals; 7514oil ES 75-14
rate parameters is often time consuming. The values of the mineral oil from Mobil Oil; MB B23/500 mircoballoons from 3M.
reaction rate parameters should be adjusted iteratively.
The calculated VoD and the radius of the detonation
front curvature needs to be iteratively compared with the
experimental results. A set of suitable values for the rate
parameters that fit predicted data well to all experimental
results can be calibrated after these iterative simulations and
comparisons. To improve the efficiency of calibration, a code
written in the Perl programming language was developed
to perform the modelling. The flow chart of the calibration
process is shown in Figure 2.
LS-Prepost in Figure 2 is the pre- and post-processor
for the LS-DYNA code, while Gnuplot is a command-line
program that can generate two- and three-dimensional plots
of functions, data and data fits. The processes in the dashed
frame can be processed by the Perl code. All the results, such
as VoD, the parameters of an ellipse for a fitted detonation
front, the figure (shape) of a fitted detonation front, the figure
(shape) of the detonation front curvature radii and the value of
the radius of curvature at the vertex of the detonation front, can
be output within two minutes after the LS-DYNA calculation.
The code can also vary the values of the parameters and run
LS-DYNA with each value as different cases. FIG 3 LS-DYNA computation set-up.

parameters. After a certain computation time, the first part

(PETN 1 in Figure 3) is deleted so that the computation can be
CALIBRATED RESULTS FOR AN EMULSION
carried out smoothly to the end. Axi-symmetric geometry is
EXPLOSIVE used in the computation.
Several sets of cylinder expansion experiments with different The parameters for the JWL EoS for reacted and unreacted
charge diameters in a polymethyl methacrylate (PMMA) explosives are from Deng, Nie and Chen (1995) and are
tube were carried out by Nyberg, Deng and Chen (1994) to listed in Table 2. The dynamic behaviour of PMMA is
investigate the VoD and detonation front curvature of an modelled with the mat_elastic_plastic_hydro model in
emulsion explosive. Eight of them recorded both VoD and LS-DYNA combined with a Grneisen EoS, and the
detonation fronts successfully. The composition and density
parameters for PMMA are listed in Table 3. The PETN that
of the emulsion explosive are listed in Table 1.
was treated as an ideal explosive was modelled with the mat_
The LS-DYNA set-up for the computations is shown in high_explosive_burn model in LS-DYNA combined with
Figure 3. In order to avoid numerical instability due to large a JWL EoS. The parameters for PETN are listed in Table 4.
deformations, the pentaerythritol tetranitrate (PETN) booster After iterative simulations, a set of suitable parameters for
consists of two parts, with both parts having exactly the same the reaction rate function in the I&G reactive model were
calibrated for the emulsion explosive (Table 2) (the columns
under Reaction rate parameters for three terms model). The

TABLE 2
Jones-Wilkins-Lee (JWL) parameters and the calibrated reaction rate parameters.

JWL parameters JWL parameters Reaction rate parameters

of the unreacted of the reaction for three terms model
explosive products
A (GPa) 4 269.3 A(GPa) 444.23 I 40.0 G2 200
B (GPa) -0.0601 B(GPa) 15.47 b 0.667 e 0.667
R1 22.628 R1 6.1045 a 0.09 f 1.00
R2 1.7907 R2 2.2733 x 7.00 z 2.0
3.0908 0.9353 G1 41.7 Fmixg 0.01
Cv(GPa/K) 1.0885e-3 Cv(GPa/K) 1.5557e-3 c 0.667 Fmxgr 0.80
T(K) 298.0 E0 (GPa) 3.348 d 0.333 Fmngr 0.80

FIG 2 The flow chart for the calibration process. y 1.12

11TH INTERNATIONAL SYMPOSIUM ON ROCK FRAGMENTATION BY BLASTING / SYDNEY, NSW, 2426 AUGUST 2015 587
C YI, U NYBERG AND D JOHANSSON

TABLE 3
Parameters for polymethyl methacrylate (Deng, Nie and Chen, 1995).

Density Shear modulus Yield stress C S1 S2 S3 Gama0 a

1185 kg/m3 1.2 GPa 76 MPa 2 572 m/s 1.536 0 0 0.97 0

TABLE 4
Parameters for pentaerythritol tetranitrate (Deng, Nie and Chen, 1995).

Density VoD CJ pressure A B R1 R2 E0

1500 kg/m 3
7 450 m/s 22 GPa 625.3 GPa 23.29 GPa 5.25 1.6 0.28 7.78 GPa

results of the simulations, together with the experimental VoD and the predicted radii of detonation front curvature at
results, are shown in Figure 4. It can be seen that both VoD and different charge diameters using the calibrated parameters are
the radii of detonation front curvature are in good agreement in good agreement with the experimental data. The calibrated
with the experimental data, which indicates that the calibrated reaction rate parameters can be used for numerical modelling
reaction rate parameters are suitable to simulate the non- for the emulsion explosive.
idealities of the emulsion explosive at density 1.165g/cm3.
ACKNOWLEDGEMENTS
CONCLUSIONS The authors would like to thank Magnus Bergh of the Swedish
The reaction rate parameters in the I&G reactive model were Defence Research Agency for his useful discussion and help
calibrated for an emulsion explosive from the experimental with the Perl programming.
detonation velocities and the detonation front curvature
radii at different charge diameters. A code written in the REFERENCES
Perl programming language was developed to make the Aldis, D F, Quirk, W and Breithaupt, R D, 1991. The effect of
parameter calibration process efficient. Both the predicted detonation curvature on cylindrical wall motion, LLNL report.
Bdzil, J B, 1981. Steady-state two-dimensional detonation, J Fluid
6.0
Mech, 108:195226.
Bdzil, J B, Aslam, T D and Short, M, 2002. DSD Front Models: nonideal
5.5 Experiment explosive detonation in ANFO, presented to 12th International
Fitting line Symposium on Detonation, San Diego.
LS-DYNA Bdzil, J B and Stewart, D S, 2007. The dynamics of detonation in
VOD (Km/s)

5.0 explosive systems, Annu Rev Fluid Mech, 39:263292.

Braithwaite, M, Cunningham, C V B and Sharpe, G J, 2006. Modeling
and numerical simulation of steady detonations in highly non-
4.5
ideal explosives, presented to 13th International Detonation
Symposium, Norfolk.
4.0 Braithwaite, M, Farran, T, Gladwell, I, Lynch, P M, Minchinton, A,
Parker, I B and Thomas, R M, 1990. A detonation problem posed
as a differential/algebraic boundary value problem, Math Engng
3.5 Ind, 3:4557.
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Christiansen, T, Wall, L and Orwant, J, 2012. Programming Perl:
Inverse charge diameter 1/ (1/mm) Unmatched Power for Text Processing and Scripting, 4th edition,
200 1184p (OReilly Media).

180 Deng, J H, Nie, S L and Chen, L, 1995. Determination of burning rate

parameters for an emulsion explosive, Swedish Rock Engineering
Radius at charge axis R (mm)

160 experiment Research (SveBeFo) report 17.

140 fitting curve Deng, J H, Nie, S L, Nyberg, U and Ouchterlony, F, 1999.
0

LS-DYNA A burning model for five emulsion explosives and some

120 applications, Swedish Rock Engineering Research (SveBeFo)
100 report 43.

80
Esen, S, 2004. A statistical approach to predict the effect of confinement
on the detonation velocity of commercial explosives, Rock Mech
60 Rock Engng, 37(4):317330.
40 Esen, S, 2008. A non-ideal detonation model for evaluating the
performance of explosives in rock blasting, Rock Mech Rock Engng,
20 41(3):467497.
0 Forbes, J W and Lemar, E R, 1998. Detonation wave velocity
10 20 30 40 50 60 70 80 and curvature of a plastic-bonded, non-ideal explosive PBXN-
Charge diameter (mm) 111 as a function of diameter and connement, J Appl Phys,
84(12):66006605.
FIG 4 Comparisons between experimental results and Forbes, J W, Lemar, E R, Sutherland, G T and Baker, R N, 1992.
numerical results. (A) Velocity of detonation versus inverse charge Detonation wave curvature, corner turning and unreacted
diameter; (B) radii of curvature versus charge diameter. Hugoniot of PBXN-111, Naval Surface Warfare Center, report.

588 11TH INTERNATIONAL SYMPOSIUM ON ROCK FRAGMENTATION BY BLASTING / SYDNEY, NSW, 2426 AUGUST 2015
 PARAMETER DETERMINATION OF AN IGNITION AND GROWTH MODEL FROM EMULSION EXPLOSIVE TESTS

Hamashima, H, Kato, Y and Itoh, S, 2004. Determination of JWL Persson, P, Holmberg, R and Lee, J, 1993. Rock Blasting and Explosives
parameters for non-ideal explosive, in Proceedings Conference of Engineering, 560 p (CRC Press: Florida).
the American Physical Society Topical Group on Shock Compression of
Price, M A and Ghee, A H, 2009. Modeling for detonation and energy
Condensed Matter (eds: M D Furnish, Y M Gupta and J W Forbes),
release from peroxides and non-ideal improvised explosives,
vol 706, pp 331334 (American Institute of Physics: Portland).
Central European Journal of Energetic Materials, 6(34):239254.
Kennedy, D, 1995. The challenge of non-ideal detonation, Journal de
Schoch, S, Nikiforakis, N, Lee, B and Saurel, R, 2013. Multi-phase
Physique IV, 05(C4):C4191-C4207.
simulation of ammonium nitrate emulsion detonations,
Kennedy, D L, 1998. Multi-valued normal shock velocity versus Combustion and Flame, 160:18831899.
curvature relationships for highly non-ideal explosives, in
Souers, P C, Garza, R and Vitello, P, 2002. Ignition and growth and
Proceedings 11th International Detonation Symposium, pp 181192
JWL++ detonation models in coarse zones, Propellants, Explosives,
(Office of Naval Research: Colorado).
Pyrotechnics, 27:6271.
Kirby, I J and Leiper, G A, 1985. A small divergent detonation theory
Takahashi, K, Murata, K, Kato, Y, Fujita, M and Itoh, S, 1999.
for intermolecular explosives, in Proceedings Eighth International
Non-ideal detonation of emulsion explosives, Journal of Materials
Symposium on Detonation, Albuquerque, pp 176186 (Office of
Processing Technology, 85:5255.
Naval Research: New Mexico).
Tarver, C M, Hallquist, J O and Erickson, L M, 1985. Modelling short
Lee, E L and Tarver, C M, 1980. Phenomenological model of shock
pulse duration shock initiation of solid explosives, in Proceedings
initiation in heterogeneous explosives, Phys Fluids, 23:23622372.
Eighth Symposium (International) on Detonation, Albuquerque,
Livermore Software Technology Corporation, 2014. LS-DYNA pp951960 (Office of Naval Research: New Mexico).
Keywords Users Manual, volume I, version 971.
Tarver, C M, Kury, J W and Breithaupt, R D, 1997. Detonation waves
Lu, J P and Kennedy, D L, 2003. Modelling of PBXW-115 using kinetic in triaminotrinitrobenzene, J Appl Phys, 82:37713782.
CHEETAH and the DYNA codes [online], technical report,
Tarver, C M and McGuire, E M, 2002. Reactive flow modeling of the
Defence Science and Technology Organisation (Australia), DSTO-
interaction of TATB detonation waves with inert materials, in
TR-1496, 34 p. Available from: <http://dspace.dsto.defence.gov.
Proceedings 12th International Detonation Symposium, San Diego,
au/dspace/handle/1947/4021> [Accessed: 29 June 2015].
pp 641649 (Office of Naval Research: San Diego).
Nyberg, U and Deng, J H, 1994. Determination of Hugoniot
Watt, S D, Sharpe, G J, Falle, S and Braithwaite, M, 2012. A streamline
parameters for an emulsion explosive matrix, Swedish Rock
approach to two-dimensional steady non-ideal detonation:
Engineering Research (SveBeFo) report 11 (in Swedish).
the straight streamline approximation, Journal of Engineering
Nyberg, U, Deng, J H and Chen, L, 1994. Measurement of detonation Mathematics, 75(1):114.
velocity and curvature of the front of the burning model
development emulsion explosive K1, Swedish Rock Engineering
Research (SveBeFo) report 6 (in Swedish).

11TH INTERNATIONAL SYMPOSIUM ON ROCK FRAGMENTATION BY BLASTING / SYDNEY, NSW, 2426 AUGUST 2015 589