Barbara ?timac,Vinko ?krlec,Mario Dobrilovi?,Muhamed Su?eska

Faculty of Mining,Geology and Petroleum Engineering,University of Zagreb,Pierottijeva 6,10000,Zagreb,Croatia

Keywords: ANFO Non-ideal detonation Reaction rate model Numerical modelling EXPLO5

ABSTRACT Ammonium nitrate and fuel oil(ANFO)based explosive is a classic example of non-ideal high explosives.Its detonation is characterized by a strong dependence of detonation parameters on explosive charge diameter,presence and characteristics of confinement,as well as incomplete consumption of explosive at the sonic point.In this work we propose a detonation model based on the Wood-Kirkwood(WK)theory coupled with the thermochemical code EXPLO5 and supplemented with reaction rate models.Our objective is to analyze the validity of the model for highly non-ideal ANFO explosives,with emphasis on effect of reaction rate models.It was found that both single-step and two-step pressure-based models can be calibrated to reproduce experimental detonation velocity-charge radius data of ANFO at radii significantly above the failure radius (i.e.for D/Did >～0.6).Single-step pressure-based model,with the pressure exponent equal to 1.4,proved to be the most accurate,even in the vicinity of the failure radius.The impact of the rate models is most evident on temporal (and spatial) distribution of flow parameters in detonation driving zone,especially when it comes to the conversion and width of detonation driving zone.

1.Introduction

Detonation is reactive wave phenomena which produces selfsustaining detonation wave (shock wave followed by a narrow chemical reaction zone) traveling at supersonic velocity(approaching 10 km/s),and reaching extremely high pressure and temperature(up to 40 GPa and 6000 K)in nanoseconds.There are two generally accepted detonation theories,based on conservation laws and hydrodynamic reactive flow models,that describe detonation process:a) Chapman-Jouguet (C-J) theory which assumes the reactions take place instantaneously (which implies that reaction zone does not exist) and b) Zeldovich-von Neumann-Doering(ZND) hydrodynamic theory that considers existence of a finite width reaction zone and a finite reaction time [1].Advances in experimental techniques and in understanding energy transfer and chemical kinetics in detonation reaction zone,have led to improved model of detonation called the non-equilibrium Zeldovich-von Neumann-Doring (NEZND) theory [2].Unlike the ZND model that assumes reactions are initiated immediately behind the shock front,the NRZND considers existence of non-equilibrium processes that precede and follow exothermic chemical reactions,i.e.considers existence of an induction period.

The C-J theory is successfully implemented in several thermochemical computer codes (e.g.Cheetah,EXPLO5,TDS,etc.) [3-5]which predict performance of explosives on the basis of formula,heat of formation and density of explosive.Explosives that behave according to the C-J theory are called“ideal explosives”.However,commercial explosives are often poorly modelled by the CJ theory[3].Detonation velocity (VOD) and pressure of these explosives,calculated theoretically applying the CJ theory,are significantly higher than experimentally measured.In addition,detonation velocity is strong function of explosive charge diameter and existence of confinement [6-10].Explosives that exhibit such behaviour are called“non-ideal”.

ANFO,widely used commercial explosive,is a typical representative of non-ideal explosives.In itself it is a mixture of 94%ammonium nitrate and 6% fuel oil.Non-ideality of ANFO is characterized by a wide detonation reaction zone width(several tens of mm),curved shock front and partial reaction at the end of detonation driving zone.In addition,ANFO also shows strong dependence of detonation velocity on charge radius and existence and characteristics of confinement.Numerous studies are devoted to the study of the impact of ANFO explosives in rock blasting application [6,11-13],as well as in metal pushing applications [14-16].

The main reason for non-ideal behaviour is a relatively long reaction time in chemical reaction zone (microseconds timescale)comparing to“ideal”explosives where reaction time is several tenths of nanoseconds.More specifically,non-ideal detonations have,due to slower reaction:a) a curved shock front;b) a nonlinear VOD vs initial density dependence;c) are charge diameter and confiner dependent;d) undergo partial reaction;e) have a limiting critical detonation failure diameter.To describe the nonideal detonation in a satisfactory way,kinetics of reaction in chemical reaction zone and radial expansion of the products(due to a long duration of reactions) must be considered.The C-J theory does not take into account existence of reaction zone;thus,it is not applicable to non-ideal detonation.The ZND theory,on the other hand,considers reaction zone and kinetic of reaction,but it does not consider radial expansion since it was developed for onedimensional flow.Thus,it is also inapplicable to non-ideal detonation.

Non-ideal detonation has been studied extensively both experimentally and theoretically.Researchers from Los Alamos National Laboratory [17-20] published a series of papers on experimental measurements of ANFO charge diameter effects on the detonation velocity of unconfined and confined charges.They have developed so-called Detonation Shock Dynamics (DSD) non-ideal detonation theory based on the motion of curved detonation and explosive geometry obtained experimentally [21].

Researchers from Lawrence Livermore National Laboratory[3,22]dealt mostly with non-ideal detonation of military explosives and they incorporated the Wood-Kirkwood(WK)slightly divergent detonation theory [23] in their thermochemical code Cheetah in order to model both ideal and non-ideal detonation[3,24].Sharpe and Braithwaite[10],Schoch and Nikiforakis[25]from University of Cambridge dealt with numerical modelling of non-ideal detonation and developed hydro-code COBRA based on quasi-polytropic equations of state and pressure-based reaction rate model.

Despite such numerous studies,there is still no generally accepted non-ideal detonation theory and model that can reliably predict the behaviour of highly non-ideal explosive,such is ANFO.In practice,detonation velocity of ANFO can be 50% of the ideal values predicted by thermochemical codes.Among the theories proposed so far,the most reliable and generally accepted is the Wood-Kirkwood (WK) slightly divergent axial flow detonation theory [23],Kirby and Leiper’s [26] slightly divergent flow theory which is an extended version of the WK theory,and the aforementioned detonation shock dynamics theory of Bdzil and Stewart[21].The WK theory,and its variants,are incorporated into several non-ideal detonation codes,e.g.CPeX,Vixen-n,DeNE,and COBRA[8].The common feature amongst these codes is the use of polytropic or quasi-polytropic equations of state and pressure-based reaction rate models.In addition,they require experimental data such as the unconfined velocity of detonation versus charge diameter data as input to calibrate the kinetic parameters and for validating the model.

One main issue in numerical modelling of non-ideal detonation of commercial explosives is the significant reliance on empirical information,use of incomplete equations of state of detonation products,overly simplistic or complex reaction rate models with a large number of adjustable parameters,and no sound theory-based models to account for the effect of confinement.

Despite its shortcomings(e.g.the lack of reliability in describing detonation of close-to-failure diameters of highly non-ideal explosives and empirical models for divergence term)the WK theory can successfully explain qualitative features of non-ideal detonation.This work aims to present a model of non-ideal detonation model based on the WK theory coupled with the EXPLO5 thermochemical code and to analyze the choice of reaction rate model on the calculation results.

1.1.Description of non-ideal detonation model

Our non-ideal detonation model is based on the Wood-Kirkwood (WK) slightly divergent detonation theory [23] coupled with EXPLO5 thermochemical code [4].Unlike the Zeldovich-von Neumann-Doering (ZND) detonation theory which considers onedimensional flow and always predicts the same value of detonation velocity as the Chapman-Jouguet (CJ) theory,the WK theory considers radial expansion of products and predicts detonation velocity as a function of relative rates of reaction and radial expansion.When the radial expansion tends to zero,the WK theory approaches the ZND theory.

The WK theory [11] starts with the Euler hydrodynamic flow equations coupled to chemical kinetics and solves the flow equations along the central streamline of the cylindrical explosive charge of infinite length.The radial expansion of products is treated as a first order perturbation to a perfect one-dimensional flow along the streamline.The theory results in a set of ordinary differential equations that describe hydrodynamic variables and chemical concentrations of reactant and products along the centre of the cylinder [3].

The Euler flow equations describe the flow in a chemically reacting medium and they express the conservation of mass,momentum,and energy.Taking the coordinate system attached to the shock front and moving with the shock velocity,the flow parameters become independent of time.Under such condition the Euler’s flow equations can be transformed into Lagrangian coordinates.The distance behind the shock (x) is related to the time after a Lagrangian particle reaches the shock (t) by equation:dx=u·dt,whereuis axial particle velocity in the moving frame.According to the WK theory,the Lagrangian equations for the steady state flow,with radial expansion treated as a first order perturbation,can be expressed in the form given by Eqs.(1)-(4)[1,3,24]:

where ρ is density,uis axial particle velocity in the coordinate system moving with shock velocity(u=D-up),pis pressure,λ is chemical concentration vector,Ris chemical reaction rate vector and ωris the rate of radial velocity(divergence term).Parameter η(so-called sonic parameter) and ψ (so-called pressure production parameter) are defined by Eqs.(5) and (6) [1]:

where σ is the termicity given by Eq.(7) [12]:

Sound velocity (C) is calculated using the energy derivatives as shown in Eq.(8) [1]:

The WK equations (Eqs.(1)-(4)) are integrated in a separate subroutine which is coupled with thermochemical code EXPLO5[4].EXPLO5 calculates the concentration and thermodynamic parameters of reaction products,the energy,the energy and pressure derivatives,the sonic and the pressure production parameters for a given conversion,pressure and density along the Rayleigh line(Fig.1).

By the integration of the WK equations,the flow properties behind the shock front are obtained (flow parameters,thermodynamic parameters and the state variables,reacted fraction of explosive and composition of products,etc.as a function of time or distance).The initial conditions for the WK equations are the state variables (pressure,density,energy,and reacted fraction of explosive)at the shock front wheret=0.The state variables at the spike and behind the shock front down to the sonic point,and the selfsustaining steady-state detonation velocity,are calculated as follows:

? The detonation velocity is treated as a known (specified)parameter.Intersection of the Rayleigh line for specified detonation velocity and the pressure on shock Hugoniot of unreacted explosive gives the von Neumann spike and corresponding pressure and volume (pVNS,vVNS) -initial conditions for specified detonation velocity.

? Once initial conditions are specified,EXPLO5 calculates temperature and energy of unreacted compressed explosive (atpVNS,vVNS,and λ=0).Then,for eachp,v,and λ state along the Rayleigh line,EXPLO5 calculates the energy and the concentration of products:

where subscripts u and p mean unreacted explosive and reaction products.

The WK subroutine integrates the WK equations and reaction rate equation to give flow parameters and reacted fraction of explosive (conversion),while EXPLO5 performs thermochemical equilibrium calculation to determine the concentration of individual products,the energy,the energy derivatives,the sonic and the pressure production parameters,and the sound velocity.

? Consumption of ANFO is kinetically controlled,determined by specified reaction rate equation for each reactant,while concentration of reaction products generated by reacted fraction of explosive is thermodynamically controlled and determined by the state of instantaneous chemical equilibrium between the products at a given,p,v,Tstate.EXPLO5 treats a mixture of unreacted explosive and reaction products assuming pressure and thermal equilibrium between all phases of the mixture.Consumption of ANFO is given by Eq.(10):

wheren1-n8are mol numbers of individual products,which are calculated by EXPLO5 for eachp,v,Tstate.

Fig.1.Schematic representation of the model.

? From the definition of the sonic parameter (Eq.(5)) and the pressure production parameter (Eq.(6)) it follows that the WK equations are finite when η=0 exactly when ψ=0.This solution corresponds to self-propagating flow.The self-sustaining steady-state detonation velocity is obtained by varying detonation velocity until two conditions are satisfied simultaneously:a) the flow is sonic (i.e.η=0) and b) the release of energy by chemical reactions is balanced by energy loss by the radial flow (ψ=0).

? We applied the minimization method described by Refs.[3,8,24],to find the sonic point.The method involves varying detonation velocity between some minimum and maximum values,specified by user,until the sonic condition is not satisfied:

whereY(D) is the merit function.

In addition to the initial condition,the non-ideal detonation model is solved simultaneously with:

-equations of state of unreacted explosive and detonation products

-thermodynamic functions of unreacted explosive and products as a function of temperature

-reaction rate model

-rate of radial expansion

1.2.Equation of state of unreacted ANFO and its detonation products

The gas phase products equation of state (EOS) is described by Becker-Kistiakowsky-Wilson (BKW) EOS incorporated in EXPLO5[4]:

where:pis the pressure,Vis the volume occupied by gaseous products (molar volume of gases)

ki-covolume ofi-th detonation product (values taken from Ref.[4])

xi=ni/nT(mol fraction ofi-th detonation product),

α,β,κ and θ adjustable constants (values taken from Ref.[4]).

The state of the unreacted ANFO and condensed detonation products is described by Murnaghan EOS given by Eq.(13)[11,27]:

where:

-V0is the molar volume of a unreacted ANFO whenp=0:V0=Mw/ρ0

-κBis the inverse of the bulk modules:

-nis the pressure derivative of the bulk modulus:n=4s-1

-C0andsare parameters in shock Hugoniot equation:Us=C0+s up

It should be emphasized that there is no published experimental shock Hugoniot data for ANFO mixtures was found [28,29],however experimental shock Hugoniot data for porous ammonium nitrate [30] and for paraffin [31] and fuel oil [32] are available.The shock Hugoniot for ANFO is usually constructed from a mass weighting of the constituents of ANFO and their shock Hugoniots[28].Dremin et al.[30] reported the shock Hugoniot for AN at 0.86 g/cm-3asUs=2.2+1.96up(in mm/μs),and shock Hugoniot of fuel oil is reported in Ref.[32]:Us=1.775+1.725up(in mm/μs)for a density of 0.837 g/cm.Reported data refers to a fixed densities and extrapolation of these data to other densities is not trivial task and requires use of phenomenological porosity model [29].

We used values of parametersnand κBreported by Esen et al.[11]:n=4.6 and κB=1.46·10-9Pa-1,valid for ρ0=0.8 g/cm.The authors derived these values from ANFO’s shock Hugoniot given by equation:Us=0.92 +1.4up.The same parameters ofnand κBare used for ANFOs having densities in the range 0.8-0.91 g/cm,while molar volume(V0)is calculated from density and molecular weight of ANFO (V0=Mw/ρ0）.

1.3.Thermodynamic functions of unreacted ANFO

Thermodynamic functions (cv,E,S,G) of unreacted ANFO are derived from the enthalpy,where the enthalpy dependence on temperature is described by a four-degree polynomial [33].The polynomial coefficients (c1-c5) are derived based on the heat capacity data [34,35] for pure ammonium nitrate atT<415 K,and taking constant heat capacity above 1000 K:c1=-13674,c2=369.28,c3=0.066585,c4=-1.24E-5,c5=8.279E-10.

1.4.Model of radial expansion

One of the key difficulties when using WK theory is the unknown axial flow divergence term.The exact nature of the radial flow can be obtained only by solving the fully coupled 2D hydrodynamics problem,which is complicated for many reasons.However,its magnitude can be estimated by semi-empirical relationships between the detonation front curvature and the explosive charge diameter [36].Such treatment is probably only adequate in the cases where the radial expansion produces a small perturbation to planar wave propagation,i.e.not at charge diameters close to the failure radius [3].

Eq.(14)proposed by Wood and Kirkwood[23]was incorporated in our non-ideal detonation model to estimate the rate of radial expansion (ωr）along the centre streamline:

which relates radial expansion rate with the shock front curvature radius (Rc),the particle velocity in the shock frame (u)and the detonation velocity (D).

The radius of curvature for a given charge diameter can be determined experimentally [37,38] or estimated by empirical equations(e.g.Refs.[3,39]which take into account the relationship between the radius of curvature,the charge radius (R0),and the failure radii (Rf) of explosive.Based on experimental data for different explosives,Souers proposed the following relationship applicable for various composite explosives[3]:

wherea=0.0316 andb=0.178.

1.5.Reaction rate models

Rate of chemical reaction in solid state is a function of concentration,temperature,and pressure (dλ/dt=f (λ,T,p).Because chemical reactions in detonation process take place under unusually high pressures and temperatures,experimental determination of reaction rate is a hard task.Thus,one of the key difficulties when implementing WK theory is that most of the relevant detailed reaction rate laws are unknown [3].

Several dozens of reaction rate models have been incorporated in hydro-codes in order to model initiation and propagation of detonation [40].Some of them are simple and insufficient,others are too complicated with numerous adjustable constants.In addition,some models are purely empirical,while some are theoretically based.For heterogeneous explosives,such as ANFO,in which initiation is known to be governed by the hot-spot mechanism,various pressure dependent models (single-step and multi-step models) have been mooted.

We have analysed the capability of single-step pressure-based(PB) reaction rate model to describe experimentally observed behaviour of ANFO explosive.This reaction rate model is most frequently used in modelling the non-ideal detonation of ANFO[10,18,27,41].The model assumes reaction rate is a function of conversion and pressure:

where λ is the mass fraction of reacted explosive(conversion),kis the rate constant,pis the pressure andb,c,anddare the rate constants.For ANFO most authors assumeb=0 andc=1,whiledandkvary.For example,for ANFO having density 0.8 g/cm3Sharpe and Braithwaite[10]usek=0.036 μs-1GPa-1.5andd=1.5 while Short et al.[18]usek=0.02 μs-1GPa-2andd=2.For density 0.88 g/cm3Wang et al.[41]usek=0.176 μs-1GPa-1.3andd=1.3.

To describe the initiation and detonation of ANFO via the hot spot mechanism,Kirby and Chan [42] proposed a two-step pressure-based (K&C) reaction rate model in which hot-spot reactions are assumed to start above some critical pressure:

where the first term in the brackets represents hot spot behaviour(when λ=0 andah=1) and the second pressure dependent bulk burning(when λ→1 andah→0).phrepresents the critical pressure below which reaction cannot proceed,constants τhand τsrepresent the burn time parameters of hot spots and bulk,respectively,andNpthe pressure power dependence constant.The switching parameter between the two terms is given by:

where ωhrepresents the fraction of explosive reacted in the first term andNadesribes how fast the transition between the two processes takes place.Schoch and Nikiforakis [25] adjusted the reaction rate parameters for ANFO having density 0.8 g/cm3to reproduce the experimentally obtained unconfined detonation velocity as a function of charge diameter to obtain the following:

? τh=31 μs GPa,τs=60.0 μs GPa,ph=0.4 GPa,

? Np=1.38,Wh=0.7,Na=9.0,C=0.667

Based on the original work of Lee and Tarver [43],Souers et al.[44]and Kittel et al.[29],proposed a two-term ignition and growth(I&G) pressure-based rate model to describe the shock initiation and detonation of ANFO bland and generic ANFO formulations.The model assumes that ignition starts at the hot-spots and propagates beyond the sites,and is given by Eq.(19):

0<λ<λigmax|0<λ<λGmaxwhere ρ is the current density,ρ0is the initial density of explosive,Iis the coefficient of ignition,G1is the coefficient of growth,ais a compression,b,c,d,x,andyare the rate constants.The first term in Eq.(19) describes hot-spot reactions and second term growth of reactions from hot-spots.According to Kittel et al.[29],ANFO reactions take place through spherical hole burning with late time grain burning behaviour.Based on small scale experiments (charge pressed into either 6.52 or 11.25 mm inner diameter and 5.72 mm long stainless steel tube to density 0.826 g/cm3and shock wave velocity is measured by microwave interferometer),and numerical modelling the authors proposed the following set of constants for ANFO:

?I=10 μs-1,a=0.2,G=6·10-11μs-1Pa-y,

? b=2/9,x=4,y=0.9,c=2/9,d=2/3,

? λigmax=0.3,λGmax=1

Various approaches are applied to calibrate reaction rate constants,but all of them rely on some experimental data.To calibrate his ignition and growth rate model,Tarver [45] used embedded pressure-gauge experiment to record pressure-time profile,and time of arrival of shock wave to a gage.He varied the rate constant until hydrodynamic calculation described well both pressure-time profile and arrival time of shock wave.Kirby and Leiper [26] used experimental electromagnetic particle velocity gauge measurements to obtain pressure-time profile in the reaction zone and burn rate-pressure measurements to calibrate their rate model.Yi et al.[46] used experimental detonation velocity and shock front curvature and hydrodynamic simulation to calibrate rate model by varying rate constants until satisfactory agreement is obtained between calculated detonation velocity and curvature radius vs.charge radius,experimentally obtained.Similar approach is used by Park et al.[47].Different approach is used by Souers [22] who develop the method of evaluation of the rate constant in simple pressure-based model from the shock front curvature radius.

2.Results and discussion

A key influence on the accuracy of non-ideal detonation models have the rection rate model,the radial expansion model,and the equations of state of unreacted ANFO and detonation products.To study effect of reaction rate models on calculation results,we calculated radial expansion rate using the Wood and Kirkwood equation (Eq.(14)) and an empirical equation based on literature experimental data to estimate radius of shock curvature (Eq.(22))and Murnaghan equation of state (Eq.(13)) for unreacted ANFO,while reaction rate models are varied.Reaction rate models are calibrated using literature experimental detonation velocity -charge radius data.

To validate the non-ideal detonation model described above,we used various experimental results on detonation velocities of unconfined ANFO charges,as well as the results of numerical modelling reported in literature.It should be emphasized that our model predicts self-sustaining steady-state detonation velocity as a function of unconfined charge diameter,and it does not model initiation phenomena nor interaction of detonation products with surroundings.

2.1.Estimation of shock front curvature radius

In order to describe the experimental detonation velocityvs.inverse unconfined charge radius profile,our non-ideal detonation model requires,as an input,an empirical relationship for the calculation of shock front curvature radius as a function of unconfined charge radius.To derive such relationship,we used Jackson and Short’s [20] relationship between detonation velocity and shock front shape obtained from rate-stick experiments.The relationship relates normal component of detonation velocity(Dn)and the shock curvature (κ) at the centre axis:

where firing parameters (B,C1-C5) are taken from Short [48]:B=37.52 mm,C2=0.05715 mm,C3=2417 mm2,C4=74 mm2andC5=2.74 mm2.The parameters are derived for ANFO having density 0.88 g/cm3andDid=5.2 mm/μs.

From the relationship between shock front curvature radius(Rc)and shock front curvature (κ) given by equation [36,49]:

where α=2 for the axisymmetric cylinder geometry,we establishedD-Rcrelationship.Then,using experimental detonation velocityvs.charge radius data from Refs.[17,20,49]we established the relationship between curvature radius and explosive charge radius(R0) (Fig.2b).

Our initial idea was to use Souer’s equation(Eq.(15))to estimate shock front curvature radius.However,we realized that Souer’s original coefficients (a=0.0316,b=0.178) predict considerable higher values of the shock curvature radii for the same charge radius comparing with experimentally observed by Jackson and Short(Fig.2b).We fittedRC-R0data derived from Jackson and Short to Souers’ relationship (R0/2Rc-Rf/R0,whereRf=38.5 mm [20])and obtaineda=0.0847 andb=0.158.These constants describe much better theR0/2Rc-Rf/R0relationship.However,at larger radii predicted shock curvature radii are still slightly higher than experimental.From bothR0/2Rc-Rf/R0plot (Fig.2a) andRc-R0plot (Fig.2b),we found that the power law describes the best the entire range ofRc-R0data,while Souer’s model overestimatesRcat larger charge radii.

For calibration of reaction rate models,we used a power law dependence of shock front curvature radius on charge radius (Eq.(22)) given by:

According to literature sources,the failure radius of generic ANFO ranges from 31.5 mm [50] to 38.5 mm [20].In our calculations we usedRf=38.5 mm.

2.2.Calibration of reaction rate constants

The effect of reaction rate model is strongly coupled with the radial expansion model and equation of state of unreacted explosive.Given that some researchers used different equations of state of solid ANFO and different radial expansion models,we recalibrated the reaction rate constants reported by other authors for pressure-based models [10,18,41] and I&G model [30].The calibration is done by varying the values of constants in selected reaction rate model until the calculatedD-R0curve satisfactorily reproduces the experimental findings.

For calibration purpose,we considered experimental detonation velocity-charge radius data of ANFO explosives reported by Catanach and Hill [19],Jackson et al.[20] and Souers and Vitello [50].Catanach and Hill’s data refer to ANFO charges which densities varying between 0.88 and 0.93 g/cm3(average 0.905 g/cm3) and charge diameters between 77 and 205 mm.Jackson et al.’s results refer to density 0.88 g/cm3and charge diameters between 77 and 205 mm,while Souers and Vitello’s data are obtained for ANFO having density 0.8 g/cm3and charge diameters between 63 and 241 mm.These experimental data used are summarized in Fig.3.

A certain degree of agreement of the results from different sources is noticeable from Fig.3.Given that different densities of explosives are used in the above-mentioned experiments,the difference in detonation velocity is expected.However,we noticed some discrepancies that cannot be attributed to density of explosive only.For example,Catanach and Hill’s detonation velocities forR0=45-100 mm are the lowest,although density of explosive they used is the highest.The discrepancies may also be related to differences in characteristics of AN prills used,to density uncertainty,as well as to experimental error.From the experimental results reported in Refs.[19,20]we noticed that for the same charge radius error in detonation velocity can reach 9% (particularly at charge radii close to the failure radius).At the same time difference between maximum and minimum detonation velocities,for densities ranging between 0.8 and 0.93 gcm-3,at radii ranging between 45 and 100 mm equals about 8% -which is within the limits of measurement uncertainty.

Fig.2. R0/2Rc vs.Rf/R0 relationship (a) and radius of shock front curvature vs.explosive charge radius (b).

Fig.3.Experimental detonation velocity from literature vs.the inverse of unconfined ANFO charge radius.

The relationship between detonation velocity and inverse charge radius may be described satisfactorily by Eyring’s equation[23]:

whereDidis the ideal detonation velocity at infinite radius,Dis the detonation velocity at charge radiusR0,AandBare fitting constants.

Since ideal detonation velocities of ANFO are achieved only at infinite diameters,they cannot be measured.Typically,they are derived by extrapolating experimentalD-1/R0data or by thermochemical calculation.This is the reason why sometimes data reported in literature are significantly different.We obtained ideal detonation velocities for studied densities by thermochemical calculations using EXPLO5:Did=4.78 km/s at ρ0=0.80 g/cm3,Did=5.08 km/s at ρ0=0.88 g/cm3,andDid=5.26 km/s at ρ0=0.905 g/cm3.That means that ideal detonation velocity increases for 380 m/s when density increases for 0.1 g/cm3.The fitting constantsAandBin Eyring’s equation are:A=15.11 andB=13.35 for ρ0=0.80 g/cm3;A=14.52 andB=17.19 for at ρ0=0.88 g/cm3,andA=23.85 andB=21.51 for at ρ0=0.905 g/cm3.

Since quite reliable equation of state of unreacted ANFO[11]and pressure-based reaction rate model [10] are available for ρ0=0.80 g/cm3,we decided to use Souers and Vitello [50] experimental detonation velocities-charge radius data for calibration of reaction rate constants.The calibration is done by varying the values of constants in a reaction rate model until the calculatedD-1/R0curve satisfactorily reproduces Eyring’s curve(Did=4.78 km/s,A=15.11 andB=13.35) derived from Souers and Vitello’s experimental data (Fig.3).An additional criterion used for calibration was that the calculated failure radius corresponds to the experimental(Rf=38.5 mm).The values of the constants obtained by the calibration are given in Table 1.

Table 1 Reaction rate models and values of constants used in calculations.

It should be mentioned that the constants in I&G_K model are taken from Ref.[29]and should be treated with caution as they are calibrated to describe shock initiation and detonation of ANFO.We have adjusted the growth term and found that experimental data can be reproduced better by usingG=2.5·10-10μs-1Pa-0.9(I&G_M).

The calculation predicts almost linear decrease of detonation velocity and pressure with charge radius above approximately 50 mm for all reaction rate models except I&G_K model with rate constants proposed by Ref.[29] and PB_2.0 model (with pressure exponent equal 2).The recalibrated I&G_M model reproduces the experimentalD-1/R0data much better (difference between detonation velocities predicted and calculated from Eyring equation is less than 5% aboveR0=45 mm).A larger discrepancy between predicted and experimental detonation velocities exists in the vicinity of the failure radius (forR0<45 mm difference goes above 10%),where detonation velocity drops below 2 km/s (Fig.4).The best agreement between predicted detonation velocity and velocity calculated by Eyring’s equation is obtained for PB_1.4 model-the discrepancy is less than 4% forR0>39 mm.

Two-step hot-spot models (K&C and I&G) display a turningpoint behaviour in the steady detonation solution as the charge radius is decreased.This means that two steady detonation velocities are possible for the same charge radius.The upper branch of the curve corresponds to a stable detonation velocity,the turning point corresponds to the failure radius,while the lower branch of the curve represents the threshold initiation conditions in which a steady detonation can be maintained.However,any fluctuation in velocity can lead either to failure of detonation or increase of detonation velocity[26].The turning point also exists for PB models if the pressure exponents is larger than 1.3.

The failure radius predicted by our calculations strongly depends on reaction rate.The PB_2.0 and I&G_K rate models overpredict the failure radii(60 and 100 mm,respectively),while PB_1.0 underpredict the radius(22 mm).These three models are excluded in further analysis as they are deemed least accurate.The most accurate failure radius is predicted by PB_1.4 and I&G_M models(predictedRf=39 mm,which is only 0.5 mm larger than experimental).

The larger deviation of predicted detonation velocities from experimental near the failure radius deserves special mention in light of the WK theory.First,it is related to the fact that the WK theory is based on slightly divergent flow and as such its use is applicable for larger charge diameters for whichD/Did→1,and where the divergence term is not large.Near the failure radius(forD/Did≈0.42),the divergence term is large and the validity of the WK theory becomes questionable.However,the fact that some reaction rate models can quite accurately describeD-1/R0near the failure radius,suggest that,with an adequate reaction rate model,the WK theory can describe behaviour of ANFO near the failure radius reasonably well (e.g.PB_1.4 model predicts detonation velocity aboveR0>40 mm with error less than 4% comparing to values calculated from Eyring’s equation).In addition,the equation of state of solid ANFO and radial expansion model both play significant roles in the accuracy of our model near the vicinity of failure radius.

Fig.5(a) and (b)illustrates the reaction rate and conversion(fraction reacted) histories obtained from different reaction rate models.The calculation is done using large charge radius (Ro=1000 mm,close to ideal behaviour).Note,the constants in models are calibrated to reproduce the detonation velocity-inverse charge radius experimental data.As expected,the hot spots based models(I&G and K&C) are characterized by very high reaction rate in the beginning of reaction(t→0,λ→0),followed by slow decomposition as the growth term is activated.For illustration,I&G predicts reaction rate of 32 μs-1,K&C model 8.8 μs-1while single-step pressure-based models predicts～1 μs-1.As a consequence of different reactions happening at different rates in the early stage,almost 30%of explosive reacts in the first 10 ns in the case of I&G model,K&Cmodel predicts 80% conversion in first 200 ns,both followed by slower reactions and ultimately wider Detonation Driving Zone(DDZ) i.e.the region between the shock front and the sonic point,and longer duration of reactions.On the other hand,single-step pressure-based model,which does not consider the hot-spot mechanism,predicts much slower reactions at the beginning.All models predict a drop in reaction rate as reactions goes to completion.As will be discussed later,differences in reaction rates and temporal rate profiles are responsible for different distribution of flow parameters within DDZ,as well as width of DDZ.

Fig.4.Detonation velocity vs. inverse charge radius for single-step (a) and two-step (b) pressure-based reaction models.

Fig.5.Effect of reaction rate model on reaction rate vs.time (a) and fraction reacted vs. time relationship (b) (unconfined charge, R0=1000 mm).

2.3.Size effect on detonation parameters and reaction zone structure

While the effect of reaction rate model onD-1/R0andpsonic-1/R0curves is not too pronounced,the effects on λsonic-1/R0andtsonic-1/R0are much more obvious(Figs.2 and 3).λsonic-1/R0andtsonic-1/R0show a turning-point behaviour as the charge radius decreases and the lower,unstable,branch is characterized by a smaller reacted fraction.AtR0→∞,λsonictends to 1 (ideal behaviour) and gradually decreases with decrease of charge radius for all rate models except I&G_M model.As shown in Fig.6,the fraction of ANFO reacted predicted by different models near the failure radius are very different;I&G_M model predicts λsonic=0.999,K&C predicts λsonic=0.72,and PB_1.4 predicts λsonic=0.52.It is obvious that the choice of reaction rate model has a significant effect on the calculated reacted fraction of ANFO.

Fig.6.Charge radius and reaction rate model effect on fraction reacted at sonic point.

The agreement exists to some extent between various rate models in terms of calculated sonic times forR0→∞;sonic times lie between 5 and 14 μs (this corresponds to 33-55 mm of DDZ width) (see Fig.7).With a decrease in charge radius,different rate models show different findings.The K&C model predicts a decrease in sonic time with decrease of charge radius up to the turning point,followed by an increase.The PB_1.4 and I&G_M models predict an increase of the sonic time up to the sonic point,with I&G_M model predicting decrease after the turning point.Hence,the use of different reaction models result in very different predicted sonic times near the failure radius;times may range from 8.5 μs(xsonic=15 mm)by K&C model to 65 μs(xsonic=68 mm)by PB_1.4 model.

The results show thatD-1/R0experimental behaviour can be reproduced reasonably well by different reaction rate models if the rate constants are properly adjusted.However,reacted fraction of ANFO at the sonic point,the sonic time,and the width of detonation driving zone depends strongly on the reaction rate model used(see Figs.8 and 9).This points to the conclusion that to develop an accurate reaction rate model with reasonably calibrated rate constants,experimental data on the structure of DDZ are needed.Unfortunately,such data are either unavailable or not reported in common literature.Instead,the usual practice is to calibrate the rate constants which describe specific phenomenon,such as for shock initiation,for steady state detonation,etc.(e.g.using embedded pressure gauges technique [45]),The rate models and rate constants obtained in this way are typically inadequate for addressing the entire problem related to non-ideal detonation.

Fig.7.Charge radius and reaction rate model effect on sonic time.

Fig.8.Comparison of spatial distribution of conversion for three different reaction rate models (R0=1000 mm,PB_1.4 rate model).

Fig.9.Comparison of spatial distribution of pressure for three different reaction rate models (R0=1000 mm,PB_1.4 rate model).

2.4.Comparison with experimental results and hydro-code calculation

Validation of our non-ideal detonation model is done comparing calculation results with literature reported experimental results and with the results from hydro-code calculations for DDZ structure (particularly duration time and width,pressure-time profile,etc.).For this purpose,we used published experimental results for confined cylindrical ANFO charge having 146 mm radius and for unconfined ANFO stack having around 2.5 m radius and containing 109 tons [16].The calculation for confined charge is done using non-ideal detonation model and two different values of shock front curvature radius,while calculation for 109 tons of ANFO is done using both ideal detonation and non-ideal detonation model.The calculation is done using the PB_1.4 reaction rate model,which has shown to be most accurate in reproducing the experimentalD-1/R0relationship (Fig.4).The results are summarized and compared in Table 2

Detonation parameters of confined ANFO charge are calculated in two ways;usingRc=549 mm (obtained by Eq.(22) for unconfined charge 146 mm in diameter) andRc=1200 mm(adjusted so to reproduce experimental detonation velocity of 4.55 km/s).Since confined curvature radius is always larger than unconfined[13],the first way may serve only as an approximation for larger charge radii,where detonation velocity approaches ideal detonation velocity.In this case detonation velocity is 94%of ideal(Did=4.85 km/s at ρ0=0.82 g/cm) and calculated detonation velocity is close to experimental(difference 4.8%).The second way is appropriate,but it requires an accurate model of prediction of radius of shock curvature for confined charges.Calculated detonation pressure is much lower than experimental (16.7%) and consequently particle velocity and sound velocities derived using Rankin-Hugoniot relationships differ about 15-20%comparing to experimental.

The comparison of calculation and experimental results for detonation of 109 tons of ANFO,2.5 m in radius,shows that almost ideal detonation is achieved.Calculation using ideal detonation model predicts detonation velocity of 4.93 km/s,which is 4%higher than experimental.However,calculation done using nonideal detonation model predicts detonation velocity of 4.83,which is 1.9% lower than experimental.Calculated detonation pressure is again about 13% lower than calculated.The fact that detonation velocities are reproduces fairly well (error below 5%)while detonation pressures differ much more(13-17%)should be viewed in the light of fact that detonation velocity can be measured with an error of less than 1%,while measured detonation pressure bay various methods span a range of 10-20% [51].

Validation is also done comparing calculation results with results of Sharpe and Braithwaite’s [10] hydro-code modelling.The authors reported calculations for the planar (ZND) detonation using the hierarchically adaptive grid hydro-code COBRA.The code uses pseudo-polytropic EOS and pressure-based reaction rate model.Our calculation is done using ideal detonation model(C-J model) and using non-ideal detonation model with large charge radius and shock curvature radius (R0=5 m mandRc=81.2 m).Using largeRc,radial expansion term in flow equations (Eq.(2)) approaches zero and non-ideal detonation model approaches the ZND ideal detonation model.The results of comparison are summarized in Table 3.

The calculations are in good agreement with Sharpe and Braithwaite’s hydro-code COBRA calculations.For illustration,calculated detonation velocity differs 0.4%when ideal detonation model is used,and 2.2% when non-ideal detonation model is applied.Detonation pressure differs 3.5% for ideal detonation model and 1.2% for non-ideal detonation model.

Good agreement is also seen for reacted fraction-distance and pressure-distance profiles obtained by hydro-code COBRA and our calculations(Figs.8 and 9).Different rate models predict different dependence of reaction rates on reacted fraction and pressure.As a result,the calculated DDZ width varies from 14.7 to 42.2 mm.It is obvious that for accurate calibration of reaction rate models,experimental pressure-time profile is necessary,along withD-1/R0.

In order to gain insight on the predictability of effect of charge radius on width of DDZ,we have compared calculated detonation velocities and parameters in DDZ(pVNS,psonic,xsonicand λsonic)for unconfined charges having different radii,with those predicted by hydro-code COBRA.The results are summarized in Table 4.

Table shows fairly good agreement for the detonation velocity and pressure forR0=75 and 100 mm (bellow 2%),while the results forR0=50 mm differ significantly (Dfor 7% andpsonicfor about 15% for PB_1.4 rate model).However,rate model PB_1.4 gives the best agreement with experimental detonation velocitywith a discrepancy less than 4.5% for all three charge radii.Some difference is observed in calculatedxsonicand λsonic.The biggest difference between our calculation and hydro-code for λsonicis about 22%(atR0=50 mm,PB_1.4),whilexsonicdiffer 29%for PB_1.4 rate model and by 75% for K&C rate model atR0=100 mm.Generally,hydro-code calculations result in slightly lower conversion and narrower DDZ at the sonic point.Since COBRA uses pressure-based reaction rate model,better agreement between Sharp’s calculation and our calculation with PB_1.4 model is expected.

Table 3 Comparison of calculated detonation parameters of ANFO with hydro-code calculations.

Table 4 Comparison of hydro-code calculations and our calculations for unconfined ANFO charges(ρ0=0.8 g/cm3).

Such results do reflect the impact of reaction rate model on DDZ,particularly onxsonicand λsonic.Besides,the differences between Sharpe’s hydro-code and our calculations could be partly attributed to the differences in equations of state used,which in turn affects reaction rate.

Fig.10.Calculated time profiles of parameters within DDZ of ANFO(R0=1000 mm,PB_1.4 rate model);particle velocity,sound velocity and temperature(a)and pressure,density,and reacted fraction (b).

2.5.Calculated reaction zone profile

Calculated temporal distribution of detonation and flow parameters within ANFO’s DDZ for nearly ideal detonation(R0=1000 mm,D/Did=0.97) is illustrated in Fig.10.

Certain aspects of Fig.10 are worth highlighting.Firstly,sonic point is attained after approximately 11 μs (32 mm) and fraction reacted at the sonic point equals nearly one,which means detonation is nearly ideal.Secondly,a rapid drop in pressure,density,particle velocity,and fast increase in conversion and temperature,is observed within the first 5 μs,followed by slower reactions towards the sonic point.It should also be noted that temperature of compressed ANFO at the von Neumann spike is relatively low(TVNS=875 K) and increases to 2860 K at the end of reaction.Relatively lowTVNSmay explain why decomposition of ANFO in the reaction zone is much slower comparing to some high ideal explosives where temperature in von Neumann spike may be two times higher.Thirdly,Fig.10a clearly illustrates the sonic condition:near the shock frontC0>(D-up),i.e.flow is subsonic,and over the time sound velocity decreases and flow velocity increases,reaching the point at which (D-up)=C0,i.e.sonic point.Beyond the sonic point flow becomes supersonic and communication of products with shock front becomes impossible,i.e.products do not contribute to the steady state shock front.

3.Conclusions

In this paper we describe a non-ideal detonation model based on the Wood-Kirkwood detonation theory coupled with thermochemical code EXPLO5 and supplemented with radial expansion and reaction rate models.The radial expansion is derived using the Wood and Kirkwood model (Eq.(14)),with the radius of shock curvature estimated by an empirical relationship based on experimental data in the literature (Eq.(22)).Three different types of reaction rate models are incorporated in the non-ideal detonation model(single-step and two-step pressure-dependent models).The reaction rate constants are calibrated using experimental data on detonation velocityvs.unconfined charge radius reported in literature.

The validation of described non-ideal detonation model is done by comparing calculated detonation parameters with experimental data,as well as with the results of reported hydro-code calculations.

It was found that the calculation results are in good agreement with both experimental data and hydro-code calculations onD-R0relationship,provided that the reaction rate models are properly calibrated.All three types of rate models used can reproduce experimentalD-R0data fairly well for unconfined cylindrical charges with radius above 50 mm(which is significantly larger than the failure radius of 38.5 mm).For radii close to the failure radius,the deviation is larger.This finding is usually attributed to the validity of the slightly divergent flow approach for highly non-ideal explosives and for radii close to the failure radius (where flow is“highly divergent”).However,the results show that the choice of the reaction rate model also plays a very important role in the accuracy of calculation in the region near the failure radius.

We found that changing pressure exponent in single-step pressure-based model our model can reproduce quite wellD-R0data in the vicinity of the failure radius and predict accurately the failure radius.For example,taking pressure exponent to be 1.4(PB_1.4 model)predicted failure radius to be 39 mm,which closely matches the experimental determined radius.Recalibrated twostep ignition and growth model (I&G_M) also reproduces very well the experimentalD-R0data and the failure radius(39 mm),but it overpredicts the conversions for larger radii compared to the other models.

While all three reaction rate models tested produce similar results forD-R0andpsonic-R0for radii significantly above the failure radius,there is a big difference between the models in terms of λ-R0andtsonic/xsonic-R0dependence.The Kirby and Chan (K&C) model predicts the smallest conversion at the same radius(eg.λ=0.836 atR0=60 mm),while I&G_M model predicts the highest conversion(λ=0.999 atR0=60 mm).This indicates a weak effect of reaction rate model on calculated detonation velocity and pressure forD/Did→1,but a very strong effect on temporal(and spatial)distribution of flow parameters within the DDZ.This is clearly visible from Fig.4 which shows a large difference in reaction rate-time curves for different rate models.The models that mimic hot-spot behaviour(K&C and I&G) predict very fast consumption of explosive at the beginning of reaction,while single-step pressure-based model predicts much slower consumption.The obtained results suggest that the calibration of reaction rate models based solely onD-R0experimental data may lead to erroneous results regarding the structure of the DDZ.Therefore,it would be necessary to calibrate reaction rate models considering also the reaction zone,unfortunately such information is difficult to obtain experimentally.

The results also demonstrate that our non-ideal detonation model can predict a broad range of characteristics of ANFO and other non-ideal explosives,including the key flow properties(D,p,V,λ,tsonic,xsonic) within the DDZ and at the sonic locus on the charge axis.In addition,the model can calculate the expansion isentrope of partially reacted explosive and the energy of detonation products available to perform mechanical work (detonation energy).It should also be noted that when large charge diameters are used,the radial expansion term in flow equations becomes negligible and the model tends towards the ideal ZND detonation theory.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

This research was supported by the Croatian Science Foundation(HRZZ),Croatia,under the projects IP-2019-04-1618 and I-2243-2017.