Modeling non-ideal velocity of detonation in rock blasting

The performance of commercial explosives is an important subject in rock blasting modeling and simulation. As a result of its non-ideal behavior, these explosives usu-ally react below their ideal detonation velocity. In these cases, the multi-dimensional effects, heterogeneities and confinement conditions become important for properly quantifying the detonation state. In this sense, an engineering approach to model two-dimensional steady non-ideal detonations for cylindrical stick explosives is used to quantify the expected detonation velocity for given reaction rate parameters and confinement conditions. Founded on an ellipsoidal shock shape approach (ESSA), the proposed model combines the quasi-one-dimensional theory for the axial solution with the unconfined sonic post-flow conditions at the edge of the explosive. A mechanistic confinement approach is coupled with the ESSA model to estimate the effect of the inert confiner on the detonation flow. Finally, the proposed model is used to estimate the expected detonation velocity of two typical commercial explosives in a number of different confinement conditions.


Introduction
Since the ability of modeling is inherent to any optimization strategy, a reliable explosive energy release quantification is required for a more realistic rock blasting simulation. In this regard, the interest of studying non-ideal detonations is considered to be a fundamental step for further downstream mining modeling. The explosive performance must be carried out using realistic approaches, such as those based on the Euler reactive flow analysis, and including Direct Numerical Simulations (DNS), slightly divergent flow theory, quasi-unidimensional analysis or streamline approximations, and others, rather than those based on ideal thermodynamic codes, like CHEETAH, W-DETCOM, ATLAS-Det, and others.
The detonation process in highly non-ideal explosives and its interaction with the inert confiner material is still a matter of research and discussion (Sharpe & Braithwaite, 2005;Sharpe & Bdzil, 2006;Sellers, 2007;Esen, 2008;Sharpe et al., 2009;Sellers et al., 2012;Braithwaite & Sharpe, 2013). Their interaction is complex, especially in cases where the acoustic velocity of the confiner is higher than the velocity of detonation. However, when combined with a mechanistic model for the confiner material, these non-ideal detonation models can predict, within the experimental scatter levels, important properties of the detonation, such as detonation velocity, pressure, specific volume, and others.
Mining explosives are strongly dependent of blasthole diameter, densities, reaction rates, confinement and others. These characteristics require the use of some classes of non-ideal detonation models to properly quantify the explosive's performance. These methods must be able to describe the reactive flow solution of the problem, including pressure profiles, densities and others. In this sense, one practical indicator of the explosive characteristics is the detonation velocity. This election is a consequence of the relatively simple method of recording the confined in-hole detonation velocity in real scale shots. On the other hand, it is a good indicator of the explosive's non-ideality (Bilgin & Esen, 1999) and can be normally associated with a large set of factors, such as explosive type, density, temperatures, reaction rates, blasthole diameter, confinement, primer size and many others (Sanchidrián & Muñiz, 2000).
Thus, the prediction of the in-hole confined detonation velocity becomes particularly important to quantify the performance of a given non-ideal explosive and confinement conditions. Herein, a non-ideal detonation model based on the Ellipsoidal Shock Shape Approach (ESSA) is used to properly describe the in-hole confined detonation velocities of two different typical mining explosives in a large set of experimental and simulated data.

Axial flow solution
Modeling non-ideal detonations are often based on the reactive Euler equations for the conservation of mass, momentum and energy The ESSA model requires a complete description of the axial flow solution to construct the ellipsoidal extension of the detonation shock front to the twodimensions. Although any axial solution could be used, such as those based on the slightly divergent flow theory (Kirby & Leiper, 1985), the Q1D model (Sharpe & Braithwaite, 2005) is adopted in this study due to its solid improvements over other axial models. One of the most important improvements is the substitution of the divergent term -which is unknown -by the axial curvature. Since its development, Q1D has been demonstrating excellent results regarding the axial flow solutions, even in highly non-ideal detonations (Sharpe et al., 2009).
The Q1D theory (Sharpe & Braithwaite, 2005) assumes that the axial solution depends only parametrically on the detonation velocity or the shock front curvature, suggesting that it is governed by a simple D n -κ law. When combined with the quadratic pseudopolytropic equation of state (Sharpe & Braithwaite, 2005) where ρ is the density; Q is the heat of reaction and the quadratic gamma γ* is given by: where ρ o is the initial density; γ o , γ 1 and γ 2 are coefficients obtained from an ideal thermo-dynamic detonation code, the following set of ordinary differential equations: where u is the velocity; ρ is the density; P is the pressure; E is the internal energy; λ is the reaction progress (λ = 0, for an unreacted product and λ = 1 for a complete reaction process); and W the is reaction rate, and the operator D / Dt = ∂ / ∂t + u•∇.
The governing equations are closed by defining the equation of state E(P, ρ, λ) and the reaction rate W(P, ρ, λ). In the Detonation Shock Dynamics (DSD) theory, a change of variables is necessary to describe the behavior of the detonation process along the normal direction n, at any point of the shock front. Thus, the Euler reactive flow equations must be transformed from the cartesian coordinate system (x,y,t) into a shock-attached coordinate system (n, ξ, t). A robust mathematical foundation of this transformation is discussed elsewhere (Stewart, 1993;Yao & Stewart, 1996;Sharpe & Braithwaite, 2005).A typical unconfined non-ideal detonation structure is presented in Figure 1. When particularized to the axis of the charge, these equations can describe the complete axial flow solution through a unique relationship between the normal detonation velocity and axial shock curvature, constituting the basis of the quasi-one-dimensional Q1D theory (Sharpe & Braithwaite, 2005). (4) can be obtained. Where u n is the normal particle velocity; ρ is the density; D n is the normal velocity of detonation; c is the sound speed; Q is the heat of explosion; γ is the adiabatic gamma; κ*=κ(1+nκ) -1 ; and W is the reaction rate, given by a pressure dependent equation where n, m and τ are fitting parameters; P ref is a reference pressure. This type of pressure-based reaction rate law presents a maximum at the shock (Cartwright, 2016).
Therefore, once defined the equation of state, the sound speed c can be obtained as where ρ=1/v. Depending of the "form" of the reaction rate equation, the number of differential equations needed to solve the problem can be reduced. In this article, since W is dependent on the pres-sure, and consequently on the specific density or volume, the pressure can be expressed as by assuming the strong shock approximation.
Thus, the set of ordinary differential Equations (4), (5) and (6) form a boundary-value problem, in other words, an eigenvalue problem in κ (D n ) for a given D n (κ). The most common resolution method is the shooting method, subjected to the jump shock conditions at the shock front and generalized CJ conditions at the sonic locus (Stewart & Yao, 1998;Sharpe, 2000;Sharpe & Braithwaite, 2005).
In the case of a highly non-ideal detonation, experimental shock measurements indicate that the detonation shock front can be well represented by an elliptic arc (Kennedy, 1998;Sharpe & Braithwaite, 2005). This observation is also supported by Direct Numerical Simulations (DNS). Additionally, Watt et al., (2009), using the Maximum Effective Entropy Theory (Byers Brown, 2002) in the construction of the detonation driving zone (DDZ), obtained better results when an ellipsoidal shock shape function was used to represent the detonation shock front. These evidences suggest that an ellipse of the form There is a unique relationship between the axial flow solution, shock shape parameters and radial dimension of the charge. Since the axisymmetric curvature κ is the solution of the set of differential Equations (4), (5) and (6) for a given velocity of detonation D n and reaction rate parameters n, m and τ, a further relationship between κ, α and β can be addressed in order to ensure a proper dependence between the shock shape parameters and axial solution. This is a necessary condition for the problem.
The shock front shape function z f can be differentiated twice in order to express the curvature at any point along of the shock. Thus, at r =0, z f '' can be related to the axisymmetric cylinder shock front curvature κ axis by the following expression which relates the axial solution with the shock shape parameters.
Once the interdependence between the shock shape parameters and axial solution is established, the first boundary condition at the edge of the unconfined charge can be defined in terms of the shock slope. This condi-tion establishes that the post-shock flow must be exactly sonic at the charge edge, r = R (Cartwright, 2016). That is, given by the shock polar analysis as: can be assumed to well represent the shock shape. Here, z and r are the axial and radial directions, respectively; and α and β are the semi-minor and semi-major axis of the ellipse. In this study, the semi-axis of the ellipse is called shock shape parameters. The origin of the adopted coordinate system is located at the central axis of the shock front, where r =0 and z =0. An assumed shock shape function allows one to calculate the slope angle and curvature κ at any point along of the shock z f .

Two-dimension model expansion
(10) where R is the charge radius, z' f =dz f /dr and γ is the polytropic gamma. Thus, the first derivate of Equation (10) must be equal to (12) in order to comply with the first boundary condition.
However, under an ellipsoidal shock shape hypothesis, many possible sonic solutions could be found at different charge radii, if Equation (12) is the only boundary condition. Consequently, in order to calculate the charge radius, an additional second and complementary boundary condition must be formulated by considering the interdependence of the shock shape parameters and axial solution.
Exploring the ellipsoidal structure of the problem, for a given set of shock shape parameters, it is expected to find a charge radius R smaller than the semi-major axis of the ellipse z f , so that R max =β. Hence, a relationship between the unconfined charge radius and the semi-major axis of the ellipse can be established as f n =R/R max , where f n is an expression which relates to the degree of non-ideality of the explosive. Comparable to other relationships such as D o / D CJ or λ, the R/R max ratio should present a similar behavior when the detonation approaches its ideal performance. Thus, f n is defined as a dimensionless expression of the following form where λ is the reaction progress at the axis; m is a reaction rate parameter; D o is the velocity of detonation; and D CJ is the thermodynamic ideal velocity of detonation; f a is a function dependent on the pressure and adiabatic coefficient.
Thus, both the charge radius R and velocity of detonation D o can be related to the shock shape parameters by combining the Equations (10), (11), (12) and (13). From these equations, one can observe how R increases when D o and λ increases. In the limit case, when the detonation approaches to its ideal behavior λ→1 and D o →D CJ , we have that R→R max . This behavior is in line with Sharpe & Braithwaite (2005) findings about how the shock locus becomes flatter at the axis when the detonation speed increases.
The strategy adopted in the present article is based on the ideas initially proposed by Eyring et al., (1949) and continued by Souers et al., (2004). A simple inspection of the diameter-effect curves resulting from unconfined and confined detonations reveals that for a given diameter, the detonation velocity increases as the confinement increases. Similarly, for a given detonation velocity, it shows how the charge radius decreases when the confinement increases. Thus, once the unconfined detonation state is known, the equivalent confined state is found to be dependent on the confinement, which leads to In order to explore some features of the ESSA model and its potential application in rock blasting simulations, several detonation cases were modeled to observe the influence of different explosives, rock confinements and blasthole diameters upon the degree of the detonation velocity. The evaluation is centered on the performance of two different explosives, an ANFO (Kirby et al., 2014) and blended Emulsion 70/30 (Sujansky & Noy, 2000), covering blasthole diameters from 165mm to 311mm. The properties of the different rock masses, where several in-hole detonation velocities were measured, are presented in Table 1. These data were presented by Esen (2008) as part of the DeNE code validation study.
where f c is the confinement factor; R u and R c are the unconfined and confined radius, respectively.
Sellers (2007)  f z is the artificial pre-compression factor due to the subsonic coupling (f z =1 for supersonic; and f z >1 for subsonic).
In the case of rock blasting dimen-sions, the confining thickness is far bigger than the critical thickness, corresponding to the case of infinite thickness confinement.  (Esen, 2008).

Explosive characterization
Before proceeding with the confined detonation cases, it is necessary to characterize the explosives involved in the study. In the ESSA model, and in most of non-ideal detonation approaches, the characterization is normally carried out in two parts.
The first step is the thermody-namic characterization of such explosive, which includes its ideal detonation velocity, heat of explosion and the quadratic coefficients of the isentropic gamma. These parameters are calculated with ideal thermodynamic codes, such as IDEX and ATLAS-Det, and are presented in The second step is quantification of their non-ideality behavior. Since these ammonium nitrate-based materials are non-ideal explosives, multi-dimensional effects and mixing heterogeneities become extremely important to properly quantify their detonation performance, which requires some knowledge about their reaction kinetics. The chemical reaction rate is not known analytically and requires the calibration of the reaction rate parameters against a set of unconfined detonation velocities, distributed in a set of different diameters.
The ESSA model presents an attractive fitting capability for unconfined detonation data. This is possible, since it is coupled with the Q1D model (Sharpe & Braithwaite, 2005) for the axial flow solution. The fitting process can be performed when experimental data -diameters and unconfined detonation velocities -are available. In this study, the ANFO data set is obtained from Kirby et al. (2014) while the EM 70/30 is coming from Sujansky & Noy (2000). The fitting strategy is straightforward. It consists in minimizing the sum of residues formed by the square of the difference between the experimental and calculated inverseradius of the charges by varying the reaction rate parameters n, m and τ. Thus, once both sets of information are available, the ESSA model can reliably reproduce the structure of any unconfined detonation within the set of experimental diameters, including the full mapping of the diametereffect curve of the explosive.
The values of n, m and τ resulting from the minimization process are presented in Table 3. In addition, it is noted that a weaker state dependency was required to experimentally adjust both explosives, since n=1.78 and n=5.77 were obtained. Cowperthwaite (1994) was the first to observe the existence of an inflection point in diameter-effect curves when n≥1.5 (Watt, et al., 2012;Cartwright, 2016). Figure 2 shows the modeled diameter-effect curves for the ANFO and EM 70/30 together with the experimental data (Kirby et al., 2014;Sujansky & Noy, 2000). The presence of an inflection point is evident in both curves. In the ANFO case, it indicates a critical diameter of 78.8 mm, which is very close to the critical experimental diameter of 79mm. On the other hand, the EM 70/30 experimental data suggest a critical diameter less than 40mm. In turn, the ESSA model identifies a critical diameter of 36mm, in perfect coherence with the experimental data. It is further noted that as the diameters increase, the non-ideal detonation velocities approximate to the ideal detonation velocity, as expected.

ESSA model application in rock blasting
Since the reproduction of real-scale detonation experiments in the laboratory are restrictive, the validation of non-ideal detonation models is normally carried out against experimental in-hole detonation velocity taken in production blasts. In real applications, the detonation velocity can be a good indicator of the degree of nonideality of a given explosive. In this study, a set of published in-hole detonation velocities (Esen, 2008), together with results obtained with DeNE code, were used to compare the prediction capability of the ESSA model. The DeNE code, developed by Esen (2008), is based on the slightly divergent flow theory, combining the pseudo-polytropic equation of state, pressure-based reaction rate law and statistical expressions to model the effect of confinement.
Thus, the ESSA model was applied in similar circumstances of application, such as blasthole diameters, rock types and explosive properties as defined by Esen (2008). The corresponding in-hole confined velocities of detonation and their comparison with the experimental and DeNE values were evaluated. The results are presented in Table 4. According to Esen (2008), the mean error in the experimental measurements was 3.5%, presenting a variation between 2.2% and 7.2%. As can be seen in Table  4, the results achieved with the ESSA model were very good, presenting errors consistent with the experimental tests. The total mean error of the ESSA model was 2.4% while the DeNE was 3.8%, as presented in Table 5. Analyzing by groups of explosives, the simulations with DeNE for the ANFO showed an average error of 4.7 %, while ESSA was 2.9%. For the blended emulsion, the errors were 3.3% and 2.0%, respectively.
The ESSA model shows a good response to the sensitivities imposed by different confinement conditions. This sensitivity characteristic is important to obtain a reliable prediction of the in-hole detonation velocity, and the fully axial reactive flow solution, in a sort of conditions where commercial explosives are normally applied. Thus, since the observed errors -shown in Table 5 -are within the range of experimental errors presented by Esen (2008), it is believed that the ESSA model is adequate for its practical application in rock blasting simulation and modeling.

References
A two-dimensional steady nonideal detonation model based on the Ellipsoidal Shock Shape Approach (ESSA) was successfully used to predict the in-hole detonation velocity under several confinement conditions. Once the reaction rate parameters were fitted, the ESSA model could provide an excellent prediction of the unconfined diametereffect curves for both ANFO and EM 70/30 explosives, including their failure diameters. The experimental in-hole detonation velocity data was presented alongside with the predictions of the DeNE (Esen, 2008) and compared with the ESSA model outcomes. The overall results indicate a better response of the proposed model to the different confinement conditions, explosives types and blasthole diameters, within the set of analyzed experimental data, endowing the ESSA model as a reliable information source for a more realistic rock blasting analysis and simulation.