Estimation of Damage Induced by Single-Hole Rock Blasting: A Review on Analytical, Numerical, and Experimental Solutions

Abstract

This paper presents a review of the existing models for the estimation of explosion-induced crushed and cracked zones. The control of these zones is of utmost importance in the rock explosion design, since it aims at optimizing the fragmentation and, as a result, minimizing the fine grain production and recovery cycle. Moreover, this optimization can reduce the damage beyond the set border and align the excavation plan with the geometric design. The models are categorized into three groups based on the approach, i.e., analytical, numerical, and experimental approaches, and for each group, the relevant studies are classified and presented in a comprehensive manner. More specifically, in the analytical methods, the assumptions and results are described and discussed in order to provide a useful reference to judge the applicability of each model. Considering the numerical models, all commonly-used algorithms along with the simulation details and the influential parameters are reported and discussed. Finally, considering the experimental models, the emphasis is given here on presenting the most practical and widely employed laboratory models. The empirical equations derived from the models and their applications are examined in detail. In the Discussion section, the most common methods are selected and used to estimate the damage size of 13 case study problems. The results are then utilized to compare the accuracy and applicability of each selected method. Furthermore, the probabilistic analysis of the explosion-induced failure is reviewed using several structural reliability models. The selection, classification, and discussion of the models presented in this paper can be used as a reference in real engineering projects.

1. Introduction

To fracture in-situ rock mass and prepare it for subsequent drilling and transport, explosion is widely used in the mining industry. In such conditions, run of mine fragmentation is assumed to be good when it is fine and loose enough to provide an efficient digging and loading operation [1]. Thus, significant attention has been drawn to estimating explosion-induced damage size in rock mass. The primary objective of research in this area has been to tailor blast fragmentation as well as to optimize mineral extraction and recovery cycle [2].

It should be noted that large amounts of fine materials are also produced by the crushed zone induced around the blast hole [3]. Thus, increasing the amount of fines multiplies handling and processing costs and, in many cases, reduces product value. Additionally, in some cases, such as quarry production, generated fines are recognized as waste. The volume of such wasted fines in Europe alone has been estimated to be about 500 million tons per year [4]. Thus, determining the size of the crushed zone and produced fines appears to be necessary [5]. Moreover, damage size at the perimeter of an excavation is of importance once the so-called drill and blast techniques are used for rock excavation [6]. In this area, minimizing explosion-induced damage is the main objective. This principle also needs to be considered, for example, in walls of drifts and other underground openings as well as slopes of surface mines. The damage penetrated through the walls and slopes is taken into account as unwanted damage or overbreak. This type of damage caused by explosion can thus have a direct impact on the stability and performance of the main structure [7]. Accordingly, diminishing such damages can:

• prevent possible damage to adjacent structures [8],
• enhance the stability of roof and side walls [9],
• improve excavation rate,
• reduce manufacturing expenses, and
• cut down operating costs [10].

In summary, by controlling the size of the crushed zone, one can optimize blasting fragmentation and minimize produced fine materials and recovery cycle. At the same time, optimizing the cracked zone can lead to a reduction in damage beyond the expected excavation boundary, control undesirable damage, and fit the explosion scheme to the geometric design [11]. That is why one main objective in rock blasting operations has been to keep unwanted damage under control [12]. To meet this objective, it is essential to understand and predict destruction caused by explosion [13].

In this paper, various methods associated with measuring and estimating explosion-induced rock damage are carefully studied and classified into different groups. To this end, this review is organized as follows: (1) the whole process of single-hole blasting is described from detonation initiation to wave propagation and rock vibration, and details are separately provided for each step (Section 2); (2) rock explosion damage is classified into different groups and illustrated schematically (Section 3); (3) different models are classified for assessing the size and severity of rock damage and presented in three different groups of analytical, numerical, and laboratory methods (Section 4); (4) the most commonly used methods are adopted in 13 case studies, and their results are compared (Section 5.1); (5) probabilistic methods are examined for calculating failure probability caused by rock explosion, and their potential differences are described compared to deterministic methods (Section 5.2).

2. A Review on Explosion Mechanism

In this section, the single-hole blasting process is examined and then the impact of induced detonation wave is described on surrounding rock medium. For this purpose, consider a single blast hole containing an explosive charge, as shown in Figure 1. Assume that a gauge is placed at point B to record explosion history. The results would be similar to those plotted in Figure 2a,b, wherein pressure–time (p–t) and pressure–distance (p–y) graphs induced by the explosion are schematically presented [4].

Detonation begins from the bottom of the borehole, i.e., point A, which corresponds to point O in the p–t graph. At this moment, the detonation pressure is still zero because a portion of the explosion is yet to be recorded at point B. Then, the detonation wave travels from the bottom of the hole to point B. This part corresponds to $OE$ on the p–t graph, wherein the pressure is still zero. Once the wave front arrives at point B, the p–t graph encounters a sudden peak and the induced pressure reaches its maximum value. This peak point is called the Von-Neumann spike. Next, the detonation wave passes point B toward point C. Consequently, the p–t graph drops sharply from E to F. Once the detonation wave touches point C, some part of the wave goes through stemming, with the remaining part reflecting back into the blast hole. Afterwards, the detonation wave moves through stemming, reaches the collar, and subsequently moves back from point D to point C. In the meantime, the pressure at point B decreases from F to G. After that, the induced wave travels toward point B and consequently the pressure drops from G to H. Following this step, the detonation wave moves toward the bottom of the blast hole and then gradually dissipates through the surrounding cracks and damages and leaks away from the system. This process continues until the borehole pressure reaches the atmospheric pressure. Correspondingly, the p–t graph is slowly reduced from point H to zero.

Pressure changes for different positions along the borehole from point B (wherein the p–t graph is at the maximum point) downward are illustrated in Figure 2b. As demonstrated, the pressure decreases sharply as it goes farther from point B due to borehole expansion, crushing and cracking of surrounding rock medium, explosion-induced gas leakage, etc.

It should be also noted that the process described above is more accurate for cases, in which the blast hole is long enough and the reflected wave in the borehole is neglected. In practice, multiple reflections of the wave from the bottom and the collar as well as the interaction between the wave and the lateral boundaries of the borehole can produce some fluctuations in pressure history. Thus, the actual p–t and p–y graphs are not perfectly smooth and exhibit some fluctuations. The next point to highlight is the role of stemming in the extension of the p–t graph. In fact, stemming causes the detonation wave to trap into the blast hole, making the detonation energy focused on fragmentation and breakage; this issue extends pressure history and consequently enhances explosion efficiency. More precisely, denoting stemming wave velocity and stemming length by $C_s$ and $L_{cd}$, respectively, the time duration of pressure history increases as $t=2L_{cd}/C_s$, provided that stemming is correctly placed. Without stemming, the explosion-induced gas tends to escape from the collar and the described t cannot be saved anymore. It will then waste energy and decrease explosion efficiency [4,14,15]. However, it should be noted that calculating the exact optimal stemming length is very difficult and challenging as the flow of energy during the explosion cannot be modeled accurately and usually accompanies with a margin of uncertainty. What can be mentioned with certainty is that the more the so-called energy leakage gaps are reduced, the better the energy released by the explosion can be used for fragmentation purposes [16].

3. Damage Pattern

Following the explosion, the pressure waves are rapidly released and strongly vibrate the rock environment [11,17]. The resulting vibration, occurring in a fraction of a second, stimulates mechanical and dynamic characteristics of rock mass [18]. This stimulation correspondingly generates a series of tensions and stresses on existing discontinuities, and also contributes to the opening and expansion of joints, depending on rock toughness [19,20].

First, the blast hole is relatively expanded [21,22]. Then, discontinuities increase and lead to formation of an unstable crushed zone due to the growth of fine cracks [23]. On the other hand, cracks affected more by the blast shock go beyond the crushed zone and penetrate radially into the surrounding environment [24,25]. Beyond the crushed zone and radial cracks, the effects of the explosion are observed in the form of ground vibration [26,27]. These three sections are shown in Figure 3.

Thus, the effects of single-hole explosion can be summarized in four steps:

• The blast hole is expanded.
• A crushed zone is formed surrounding the blast hole.
• Radial cracks penetrate through the rock, causing a cracked zone.
• Explosion-induced waves affect the surrounding environment, producing some ground vibrations.

The pattern of damage generated around the explosion point is initially found in practical projects but later proved in experimental models. For instance, Olsson and Ouchterlony [28] showed the pattern generated in experimental models. However, it should be noted that, in practice, the zones mentioned above are interconnected without any sharp boundaries that help in distinguishing them from one another. However, the definition provided for damage pattern can greatly help in establishing models and calculating results.

4. Estimation of Induced Damage

A shock wave is initially generated once an explosion charge is fired. Subsequently, a stress wave affects the surrounding environment, creating two damage zones near the blast hole, namely crushed and cracked zones. As discussed, sizes of these two damage zones are of importance for an optimal blast design. As far as the project client can estimate the size of damage zones (i.e., crushed and cracked zones) as a function of input parameters such as rock properties and explosive characteristics, the optimal values of input parameters can be obtained using a blast design optimization. This optimization can be done through a try-and-error process to obtain the optimal values of target parameters or can be mathematically implemented in an optimization algorithm [29,30,31].

Different researchers have proposed various methods to approximate the induced damage. In a general case, the size of a damage zone can be assumed as a function of input parameters such as:

$$r=f({\theta }_1,{\theta }_2,\dots ,{\theta }_n),$$

(1)

where r is the damage zone radius and ${\theta }_1$ to ${\theta }_n$ represent the input parameters (either rock or explosion characteristics). The most important input parameters are outlined in

Table 1. Main parameters involved in explosion-induced rock damages.

No.	Sources	Parameters	Description
1	Rock characteristics	$E_d$	Young’s modulus of rock
2		${\nu }_d$	Poisson’s ratio of rock
3		${\sigma }_c$	Uniaxial compressive strength of rock
4		$F_c$	Confined compressive strength of rock
5		T	Tensile strength of rock
6	Explosive characteristics	${\rho }_0$	Unexploded explosive density
7		$D_{CJ}$	Ideal detonation velocity
8		$r_0$	Blast hole radius
9		$Q_{ef}$	Effective energy of explosive

.

This form of damage estimation can be simply applied to a particular case study with virtually no complexities. However, it is not always easy to provide a close-form function such as $f(•)$ since several sources of uncertainties are included in explosion-associated problems [32,33,34]. Numerous research studies are also available in the related literature to estimate damage in rock and soil media [35,36,37]. In an overall view, these methods can be categorized into three main groups: analytical, numerical [38], and experimental [39].

Approaches toward problems in each of these three methods are not the same. In analytical techniques, a parameter such as peak particle velocity (PPV) [40,41,42], borehole pressure [43,44], or explosion pressure [45] is generally presented as a critical factor to estimate the size of a damage zone. Next, two different approaches are used to provide a solution; either an analytical calculation is employed to directly determine both the critical parameter and the damage zone, or the relationship between the parameter and the damage zone is estimated and the rest of the problem remains unsolved for the reader. In numerical methods, however, an algorithm such as finite element method (FEM), finite difference method (FDM), discrete element method (DEM), etc. is used to evaluate changes in the stress field surrounding an explosion point and examine consequent issues including induced cracks and damages [46,47,48]. In experimental approaches, some laboratory or in-situ tests are utilized to develop an empirical relationship to estimate the size and dimension of damages [49].

In the following, each of these three methods is separately addressed, and then related previously developed research works are listed and explained in more detail.

4.1. Analytical Approach

As previously described, in analytical approaches, a feature of a model is selected as a main parameter, and it is determined how this parameter is distributed around the blast hole. The relationship between the parameter and rock damage is then examined so that damage size can be measured for each parameter value. Peak ground velocity and borehole pressure are assumed as two parameters widely used for this purpose [50]. Analytical approaches based on PPV and borehole pressure are discussed in Section 4.1.1 and Section 4.1.2, respectively.

4.1.1. Damage Prediction Using PPV

PPV is known as one of the critical parameters, used by several researchers, to estimate damage zones [51]. Accordingly, damage rate can be predicted if PPV is estimated in different areas in rock environments. Some of the PPV-based criteria for blast-induced damage in rocks are presented in

Table 2. PPV-based criteria for blast-induced damage in rock (adopted from Bauer and Calder [52]).

PPV (mm/s)	Effects of Damage
<250	No fracture of intact rock
250–635	Occurrence of minor tensile slabbing
635–2540	Strong tensile and some radial cracking
>2540	Complete break-up of rock mass

and

Table 3. PPV-based criteria for blast-induced damage in rock (adopted from Mojtabai and Beattie [53]).

Rock Type	Uniaxial Strength (MPa)	RQD (%)	PPV (mm/s)
			Minor Damage	Medium Damage	Heavy Damage
Soft schist	14–30	20	130–155	155–355	>355
Hard schist	49	50	230–350	305–600	>600
Shultze granite	30–55	40	310–470	470–1700	>1700
Granite porphyry	30–80	40	440–775	775–1240	>1240

.

Two requirements need to be met when using PPV in damage estimation and structural control: (1) the PPV at the desired location should be predicted; (2) the relationship between the PPV and the damage state (i.e., fragility curve) should be provided. In practice, PPV vectors surrounding a blast hole are difficult to be developed, and there are not many sources available in this area. To further explain this issue, the proposed model by Holmberg and Persson can be noted [40,54], since they offered the following equation to estimate PPV:

$$V=K\frac{W^{\alpha }}{R^{\beta }},$$

(2)

where V shows PPV, K, $\alpha$, and $\beta$ are the empirical constants, W is the charge weight unit, and R denotes the distance unit from the charge. However, this equation is developed for areas far from the explosion point. In fact, the given equation is applicable for areas, in which R is so large that it makes blast hole dimension negligible, and cannot be used for areas close to the explosion point, where a blast hole dimension should be considered.

To solve this problem, Holmberg and Persson [54] assumed that the explosive charge has a cylindrical shape. Accordingly, they divided the explosive charge into small pieces so that each of the pieces has a length of $dx$ and unit charge concentration of q (kg/m). Then, they stated that PPV in an arbitrary point such as ($r_0$, $x_0$) can be calculated as follows:

$$V=K{\left(q{\int }_T^{H+J}\frac{dx}{{(r_0^2+{(x-x_0)}^2)}^{\frac{\beta }{2\alpha }}}\right)}^{\alpha },$$

(3)

where T is the stemming depth (m) [55], H is the charge length (m), and J is the subdrill (m) [56], as demonstrated in Figure 4.

The parameter $\alpha$ was assumed as follows:

$$\alpha =\frac{\beta }{2}.$$

(4)

After integrating from Equation (3), PPV was explicitly presented as follows:

$$V=K{\left(\frac{q}{r}({\tan }^{-1}\left(\frac{H+J-x_0}{r_0}\right)-{\tan }^{-1}\left(\frac{T-x_0}{r_0}\right))\right)}^{\alpha }.$$

(5)

The values of K, $\alpha$, and $\beta$ for hard rock mass are 700, 0.7, and 1.4, respectively. Thus, by having q, the PPV amount can be calculated for any desired distance (r). In this respect, Changshou Sun [6] presented

Table 4. Damage type based on induced PPV in a Scandinavian bedrock.

PPV (m/s)	Tensile Stress (MPa)	Strain Energy (J/kg)	Typical Effect in Hard Scandinavian Bedrock
0.7	8.7	0.25	Incipient swelling
0.1	12.5	0.5	Incipient damage
2.5	31.2	3.1	Fragmentation
5	62.4	12.5	Good fragmentation
15	187	112.5	Crushing

for Scandinavian bedrock to approximate rock damage based on the induced PPV.

Although this appears to be a very simple and practical method, later, Hustrulid and Lu [57] reiterated that the main integral suffered from a basic mistake and the accurate form of the integral was:

$$V=K{q}^{\alpha }{\int }_T^{H+J}{\left(\frac{dx}{{(r_0^2+{(x-x_0)}^2)}^{\frac{\beta }{2\alpha }}}\right)}^{\alpha },$$

(6)

which could not be solved obviously in an analytic manner. Moreover, the proposed method has other problems; for instance:

• Only the magnitude of the PPV is considered and the direction of the PPV is neglected.
• Only the explosive weight is taken into account, and other characteristics are ignored.
• To determine the parameters K, $\alpha$, and $\beta$, further laboratory or in-situ tests are required, which are difficult to conduct.

For such reasons, these types of approaches have not been highly welcomed by research communities.

4.1.2. Damage Prediction Using Borehole Pressure

Borehole pressure is known as one of the common parameters in the estimation of rock damage. Accordingly, it is believed that the initiation and propagation of cracks in rocks are due to severe stress caused by explosion waves. Thus, borehole pressure is being used by many researchers to directly estimate the crushed zone size in rock environments. The next sections present a summary of such studies.

Mosinets’ Model

Mosinets et al. [58] expressed the radius of damage zones (either cracked or crushed ones) surrounding the blast hole by the following equation:

$$r=k\sqrt[3]{q},$$

(7)

where r is the damage zone radius, k represents the proportionality coefficient, and q refers to the charge weight in the TNT equivalent. Each damage zone also has its own coefficient k. For the crushed zone, the coefficient k is as follows:

$$k=\sqrt{\frac{V_s}{V_p}},$$

(8)

where $V_p$ shows the longitudinal wave velocity and $V_s$ denotes the transverse wave velocity.

Drukovanyi’ Model

Assuming an isotropic and incompressible granular medium with cohesion (derived from Il’yushin’s model [59]), Drukovanyi et al. [60] examined behaviors of rocks in the zone of fine crushing. Considering a plane strain for detonation of a column of explosive material in rock mass, they theoretically developed the following relation to determine the crushed zone radius ($r_c$) close to the blast hole:

$$r_c=r_0{\left(\frac{P_b}{-\frac{C}{f}+({\sigma }_c+\frac{C}{f})L^{\frac{f}{1+f}}}\right)}^{\frac{1}{2\gamma }}\sqrt{L},$$

(9)

where $r_0$ stands for the borehole radius (mm), $P_b$ is the borehole pressure (Pa), C shows cohesion (Pa), f refers to the internal friction coefficient ($f=\tan \left(\phi \right)$), ${\sigma }_c$ is the uniaxial (unconfined) compressive strength (Pa), $\gamma$ denotes the adiabatic expansion constant of explosive, and L represents a constant defined by:

$$E=\frac{\frac{E}{1+\nu }}{{\sigma }_c\left(1+\ln \left(\frac{{\sigma }_c}{T}\right)\right)},$$

(10)

where E is the Young’s modulus (Pa), $\nu$ refers to the Poisson’s ratio, and T shows the tensile strength (Pa) of the rock. Drukovanyi’s model assumes that the borehole pressure is calculated as follows:

$$P_b=\frac{1}{8}\rho D^2,$$

(11)

where $\rho$ shows the explosive density and D is the detonation velocity. They also referred to Il’yushin’s model and considered $\gamma =3$, $\rho =0.9$ $\mathrm{gr}/{\mathrm{cm}}^2$, and D = 4000–6000 m/s.

It was noted by Drukovanyi et al. [60] that damage predicted using this model can give values higher than reality to rocks by the compressive strength less than 100 MPa. Moreover, it has been reported in the related literature that this approach is limited to cases, where the main mode of failure is compression [6].

Senuk’s Model

Senuk [61] developed the following relation to estimate the cracked zone radius in the vicinity of a cylindrical explosive charge:

$$r_c=k{r}_0\sqrt{\frac{P_b}{T}},$$

(12)

where k is the stress concentration factor in sharp cracks and joints, which is assumed to be approximately equal to 1.12, $P_b$ shows the blast hole pressure approximated by $P_b={\rho }_0{D}_{CJ}^2/8$, and T represents the rock tensile strength. Thus, Equation (12) can be rewritten as follows:

$$r_c=1.12r_0\sqrt{\frac{{\rho }_0{D}_{CJ}^2}{8T}}.$$

(13)

Szuladzinski’s Model

According to the hydrodynamic theory for rock explosion, Szuladzinski [62] proposed a model to predict the crushed zone radius around a blast hole. In this model, the rock environment in the vicinity of the blast hole was assumed as an elastic medium with cracking and crushing capabilities. The effective explosion energy in the model was also roughly assumed as two-thirds of the total explosive energy, and decoupling effects were ignored. The relation presented in this study for the crushed zone radius is as follows:

$$r_c=\sqrt{\frac{2r_0^2{\rho }_0{Q}_{ef}}{F_c^{\prime }}},$$

(14)

where $r_0$ (mm) is the borehole radius, ${\rho }_0$ ($\mathrm{g}/{\mathrm{mm}}^3$) shows the explosive density, $Q_{ef}$ (Nmm/g) represents the effective explosive energy (assumed to be two-thirds of the complete reaction heat), and $F_c^{\prime }$ (MPa) refers to the confined dynamic compressive strength of rock mass (assumed to be approximately eight times of the unconfined static compressive strength, ${\sigma }_c$ [63,64,65,66]).

SveBeFo Model

The Swedish Engineering Research Organization (SveBeFo) conducted several studies on the initiation and propagation of cracks in rock environments under explosion load. Based on the results of these studies, Ouchterlony [67] presented the following relation to calculate the length of induced radial cracks:

$$2\frac{r_{co}}{d_h}={\left(\frac{P_b}{P_{b,cr}}\right)}^{2/(3{\left(\frac{D}{c}\right)}^{0.25}-1)},$$

(15)

where $r_{co}$ denotes the unanchored radius of the cracked zone, $d_h$ shows the borehole diameter, $P_b$ is the borehole pressure, D represents the velocity of detonation (VOD), c is the sound speed in a rock environment, and $P_{b,cr}$ denotes an experimental parameter showing critical borehole pressure which can be estimated by the following equation:

$$P_{b,cr}=3.3\frac{K_{IC}}{\sqrt{d_h}},$$

(16)

where $K_{IC}$ is the fracture toughness of rock mass and $d_h$ shows the borehole diameter. Accordingly, Ouchterlony [67] introduced the following equation to estimate the borehole pressure:

$$P_b=\frac{{\gamma }^{\gamma }{\rho }_0}{{(\gamma +1)}^{\gamma +1}}{D}^2{\left(\frac{d_e}{d_h}\right)}^{2.2},$$

(17)

where $\gamma$ is an isotropic exponent for a specific explosive (1.254–2.154), ${\rho }_0$ denotes the explosive density, D stands for VOD, and $d_e/d_h$ represents the ratio of explosive diameter to borehole diameter called the decoupling ratio. Later, Ouchterlony et al. [68] provided the following correction coefficients to improve the given method:

$$r_c=R_{co}{F}_h{F}_t{F}_r{F}_b,$$

(18)

where $R_{co}$ stands for the corrected damage zone radius, $F_h$ is correction for hole spacing, $F_t$ shows correction for time spread in initiation, $F_r$ stands for correction for wet holes, and $F_b$ is correction for fracturing. These coefficients could be determined based on the concept of fracture mechanics, which is complicated for conventional engineering design. Moreover, it is not easy to determine rock fracture toughness, i.e., $K_{IC}$, especially for weak rocks. These issues led to limited applications of the given method in research works and consequently no extensive use of the method in practical designs.

Quasi-Static Model

Assuming a balance between borehole pressure and stress distributed in the surrounding rock medium, Sher and Aleksandrova [69,70] provided a model to predict crack size around a cylindrical borehole. In this model, a dynamic process was approximated by a quasi-static method, and the following equations were then proposed to estimate the radius of the cracked zone:

$$(\frac{Y}{E\alpha }-\frac{q}{E}){\left(\frac{r_d}{r}\right)}^{\alpha /(1+\alpha )}-\frac{Y}{E\alpha }-\frac{P_h}{E}=0,$$

(19)

$$\frac{q}{E}=-\frac{\frac{{\sigma }_c}{E}+2(\frac{\sigma }{E}\times \frac{{\sigma }_c}{T})}{1+\frac{{\sigma }_c}{T}},$$

(20)

$$\frac{u_b}{r_h}=-(1+\nu )\frac{r_d}{r_h}\frac{q+2\sigma (1-\nu )}{E},$$

(21)

$${\left(\frac{r}{r_h}\right)}^2-1={\left(\frac{r_d}{r_h}\right)}^2-{(\frac{r_d}{r_h}-\frac{u_b}{r_h})}^2,$$

(22)

$$\frac{P_h}{E}=\frac{P_{CJ}}{E}{\left(\frac{r}{r_h}\right)}^{-2{\gamma }_1},r\le r^{*},$$

(23)

$$\frac{P_h}{E}=\frac{P_{CJ}}{E}{\left(\frac{r}{r_h}\right)}^{-2{\gamma }_1}{\left(\frac{r^{*}}{r_h}\right)}^{-2{\gamma }_2},r>r^{*},$$

(24)

where $r_h$ denotes the initial hole radius, r shows the final hole radius, $u_b$ is the elastic deformation of rock, $r/r_h$ represents the ratio of final radius to initial radius of hole, $r^{*}$ refers to the radius at which the adiabatic constant changes, ${\gamma }_1$ is the initial adiabatic expansion constant (${\gamma }_1=3$), and ${\gamma }_2$ signifies the final adiabatic expansion constant (${\gamma }_2=1.27$). Both $\alpha$ and $\beta$ parameters can be calculated as follows:

$$\alpha =\frac{2\sin \varphi }{1-\sin \varphi },$$

(25)

$$\beta =\frac{2c\times \cos \varphi }{1-\sin \varphi },$$

(26)

where c is the cohesion and $\varphi$ refers to the internal friction angle. Based on the equations presented above, Hustrulid [71] provided the following process to determine the cracked zone radius:

Step 1 

Calculate $q/E$ from Equation (20)

Step 2 

Approximate a value for $r_d/r_h$ (this value is approximated in this step and later modified in a cyclic process)

Step 3 

Substitute $q/E$ and $r_d/r_h$ in Equation (21) and calculate $u_b/r_h$

Step 4 

Substitute $u_b/r_h$ in Equation (22) and determine $r/r_h$

Step 5 

Select one of Equations (23) or (24) based on a comparison between $r/r_h$ and $r^{*}/r_h=1.89$ and then calculate $P_h/E$

Step 6 

Substitute $P_h/E$ in Equation (19) to assess if equality is achieved (if so, $r_d/r_h$ is the final answer. Otherwise, the steps 2–6 should be repeated until the final answer is reached).

Djordjevic’s Model

Based on the Griffith’s failure criterion, Djordjevic [72] developed a model mostly applicable for brittle rocks [2]. The crushed zone radius proposed in this study is:

$$r_c=\frac{r_0}{\sqrt{24T/P_b}},$$

(27)

where $r_0$ (mm) is the blast hole radius, T (Pa) shows the tensile strength of rock material, and $P_b$ (Pa) represents the borehole pressure.

Kanchibotla Model

Kanchibotla et al. [73] considered damage around a blast hole as a function of blast hole radius, explosion pressure, and uniaxial compressive strength of rock mass. Using these parameters, they proposed a relation to estimate the crushed zone radius as follows:

$$r_c=r_0\sqrt{\frac{P_d}{{\sigma }_c}},$$

(28)

where $r_0$ (mm) is the borehole radius, $P_d$ (Pa) refers to the detonation pressure, and ${\sigma }_c$ (Pa) denotes the unconfined compressive strength of rock.

Johnson’s Model

Johnson [74] considered four different zones near a blast point including borehole, crushed zone, cracked zone (also named as transition zone), and no-damage zone (also labeled as seismic zone). To calculate the crushed zone radius in this study, the following formula was proposed:

$$Q_{cd}=P_b\sqrt{\frac{r_0}{r_c}}{e}^{-(r_c-r_0)\lambda },$$

(29)

where ${\sigma }_{cd}$ is the dynamic compressive strength of rock, $P_b$ represents the borehole pressure, $r_c$ shows the crushed zone radius, $r_0$ stands for the borehole radius, and $\lambda$ is the crush damage decay constant determined by laboratory experiments.

Modified Ash’s Model

Hustrulid [71] improved the Ash’s model [63,64,65,66] and provided the following equation to approximate the size of the cracked zone around a blast hole:

$$r_c=25r_0\left(\frac{d_e}{d_h}\right)\sqrt{RBS}\sqrt{\frac{2.65}{{\rho }_r}},$$

(30)

where $r_0$ is the borehole radius, $d_e/d_h$ shows the decoupling ratio, $RBS$ stands for the relative bulk strength (compared with ammonium nitrate/fuel oil (ANFO)), and ${\rho }_r$ refers to the rock density. $RBS$ can be calculated as follows:

$$RBS=\frac{{\rho }_0{S}_{ANFO}}{{\rho }_{ANFO}},$$

(31)

where $S_{ANFO}$ is the weight strength of explosive relative to ANFO, ${\rho }_{ANFO}$ stands for the density of ANFO, and ${\rho }_0$ shows the explosive density.

4.2. Numerical Approach

When a rock is idealized as a continuous medium, continuum-based approaches such as FEM [75,76,77,78,79,80,81] or FDM [82,83,84,85,86,87] are usually employed to simulate the model. However, the rock is modeled as a set of structural units (such as springs, beams, etc.) or as separate particles bonded at contact point if modeling of discontinuities is required. In such cases, DEM [88,89,90,91,92,93], the bonded particle model [94] or hybrid methods [95,96,97,98,99] are used for numerical analysis. Typically, due to the high volume of calculations required in such methods, calculations are performed using software packages to simulate the model. Various packages have thus been developed for geotechnical modeling. In this regard, two famous DEM codes, UDEC and 3DEC, are widely utilized for modeling two- and three-dimensional cases of jointed rock, respectively.

Wang and Konietzky [46] used UDEC to study the initiation and propagation of damage around the blast hole in rock. They first used general-purpose multiphysics simulation software package (LS-DYANA) to model an explosion in intact rock and calculate the explosion load imposed on the borehole wall. Then, they returned to UDEC and modeled a jointed rock environment. They assigned the obtained load as a radial velocity history to the borehole wall and studied the damage evolution. Their coupled model is schematically shown in Figure 5.

To model the jointed rock environment in UDEC, Wang and Konietzky [46] used two sets of joints: a randomly-generated series of polygonal joints (e.g., Voronoi joints) and two series of orthogonally aligned joints dipping at $-{15}^{\circ }$ and ${75}^{\circ }$. Having the load exerted on the blast hole, they began to study the growth of damage in rock. It was observed that some small damage formed around the blast hole and then penetrated through outer layers. By passing the time, the damage size continuously increased until $t=5$ ms that it stopped growing further and reached at its maximum size. Later, it was noted that, at the time $t=5$ ms, the maximum length of cracks reached to 2.5 m while the crushed zone radius varied between 7 to 9 times of the blast hole radius. This estimation is moderately higher than that in similar studies approximating the ratio below 5 between the crushed zone radius and the borehole radius ($r_c/r_o$).

Using the analytical code of the Realistic Failure Process Analysis 2D (RFPA2D), Liu et al. [100] analyzed rock mass behavior under explosion assuming different geometrical properties for joints including average distance from blast site to joints, length of joints, number of joints, and relative angles of joints. These properties are schematically shown in Figure 6.

The results of this study showed that the effect of joints on blast-induced damage is slowly reduced as the average distance of joints from the blast hole increases. It was also observed that rock masses with long joints could experience more damage than those having short joints. Their results also revealed that the damage caused by the blast was significantly greater in rock masses with joints than in cases without joints. However, little change might be observed in damage intensity caused by an explosion when the number of joints exceeded 3. It was also concluded that joints could highly facilitate the propagation of cracks and play a significant role in determining their shapes, dimensions, and directions of propagation.

Moreover, Lak et al. [101] used a hybrid finite difference-boundary element method to investigate the propagation of blast-induced cracks around a wellbore. To this end, they planned two separate steps. First, they investigated the formation of radial cracks due to the propagation of shock waves caused by explosion using the dynamic finite difference method. Importing outputs into the second step, they then modeled the propagation of cracks caused by gas expansion using the quasi-static boundary element method.

Two types of cracks were considered in this study: type I and type II. The cracks along the direction of maximum horizontal stress (${\sigma }_{Hmax}$) were considered as type I, and those along the direction of minimum horizontal stress (${\sigma }_{Hmin}$) are marked as type II. The results of this study showed that the ratio of ${\sigma }_{Hmin}/{\sigma }_{Hmax}$ has an obvious effect on the radial cracks propagation. When ${\sigma }_{Hmin}/{\sigma }_{Hmax}$ has a small value near 0, the length of type I cracks is larger than that of type II cracks (roughly about twice, $a_2/a_1=0.5$). With the increase of the ratio of ${\sigma }_{Hmin}/{\sigma }_{Hmax}$, the length of type I cracks decreases and the length of type II cracks increases. Finally, in hydrostatic stress conditions, i.e., cases where ${\sigma }_{Hmin}/{\sigma }_{Hmax}=1$, the length of type I and type II cracks turn out to be equal ($a_2/a_1=1$).

4.3. Experimental Approach

Explosion-induced damage in rocks occurs through a complex process. Therefore, laboratory studies have been exploited as one of the main methods to investigate this problem in the past few decades given the complexity of research in this area [102]. It should be noted that experimental studies of explosion phenomena as well as subsequent failures mainly deal with two main aspects [103]:

• Primary cracks due to the high amplitude of stress waves
• Further development of cracks due to gas penetration

In a study on stress-wave, Lownds [104] conducted a set of reduced-scale experiments in granite. A sensor hole was also placed at a distance of 150 mm away from a borehole with a diameter of 32 mm. A pressure transducer was subsequently placed in a water-coupled sensor hole with the same diameter as the borehole to measure stress waves. They recorded the effect of different coupling media on stress waves using the installed pressure transducer. In a similar study, Talhi et al. [105] measured the relative magnitude of stress pulses by placing a pressure gauge in a water-coupled sensor hole with a diameter of 8 mm to check the variation of stress waves caused by changes in decoupling ratio and borehole length. Through these methods, peak pressure at a relatively close distance to the borehole was measured providing a database of relative measurement for different borehole conditions. In addition, Teowee and Papilon [106,107] utilized piezoelectric sensors to measure the sympathetic pressure of adjacent deck charges and boreholes. For this purpose, they compared peak pressure between these piezoelectric gauges and concluded that the given sensors could measure peak pressure up to 138 MPa. Due to the relative size of the manometer and the applicability of measuring high pressure, piezoresistive gauges, such as carbon composite resistors, have been employed by different researchers to investigate pressure pulse caused by stress wave propagation [44,108].

Measuring damage size and pattern using induced gas pressure and penetration has also been investigated by different researchers. For example, the direct measurement of crack length was performed by Olsson and Bergqvist [109,110], Deghan Banadaki [111], and Nariseti [112] to reflect on blasting-induced rock fractures. Paventi and Mohanty [113] also shed light on the effect of coupling medium on crack pattern through laboratory-scale tests of samples with 10-cm diameter. Boreholes with diameters of 6 mm and 10 mm were thus equipped with pentaerythritol tetranitrate explosive in central position while containing different coupling media such as air, water, and clay. The results showed that fracture intensity was smallest in air-coupled boreholes followed by clay-coupled boreholes. Additionally, they reported that water-coupled boreholes led to intense fracturing near boreholes. Furthermore, Yamin [114] conducted reduced-scale experiments to examine damage extent around boreholes. To this end, pressure sensors were placed in sensor holes on concentric circles at different distances surrounding a borehole to measure gas pressure. It was consequently reported that gas penetration was observed from a 75-mm blast hole charged with a 40-mm emulsion cartridge up to a distance of about 15 borehole diameters. Additionally, Yamin [114], McHugh [115], and Brinkmann [116] investigated rock damage following gas penetration.

It should be noted that laboratory and experimental models to directly determine damage zones in rock under explosion load are difficult and costly [117]. The test site should also be safe and equipped with required facilities to measure rock behavior. Providing such conditions, however, is not simple. Therefore, there are few studies available in this area.

One of the existing approaches in this field is Esen’s model [2]. To measure the crushed zone size around a blast hole, Esen et al. conducted a detailed laboratory test on 92 samples mostly made up of concrete with dimensions of 1.5 m × 1 m × 1.1 m. They consequently defined a crushing zone index ($CZI$) based on the crushing process around the explosion point as follows:

$$CZI=\frac{P_b^3}{K\times {\sigma }_c^2},$$

(32)

where $P_b$ is the blast hole pressure (Pa), K stands for the stiffness of rock mass (Pa), and ${\sigma }_c$ shows the uniaxial compressive strength of rock (Pa). The values of $P_b$ and K are calculated as follows:

$$P_{CJ}=\frac{{\rho }_0\times D_{CJ}^2}{4},$$

(33)

$$P_b=\frac{P_{CJ}}{2},$$

(34)

$$K=\frac{E_d}{1+{\nu }_d},$$

(35)

where $P_{CJ}$ is the ideal blast pressure (Pa), ${\rho }_0$ shows the unexploded explosive density (kg/m³), $D_{CJ}$ represents the ideal detonation velocity (m/s), $E_d$ denotes the dynamic Young’s modulus of rock (Pa), and ${\nu }_d$ refers to the dynamic Poisson’s ratio of rock. These relationships are also used to approximate blast hole pressure and rock stiffness. Where more accurate values of these parameters are available through direct measurement or numerical modeling, they can be used instead of the presented equations [118]. After calculating $CZI$, Essen et al. [2] found a power relationship between this factor and the crushed zone radius as follows:

$$\frac{r_0}{r_c}=1.23\times CZ{I}^{-0.219},$$

(36)

where $CZI$ is the crushing zone index, $r_c$ represents the crushed zone radius, and $r_0$ shows the blast hole radius. The crushed zone radius is then calculated as follows:

$$r_c=0.812\times r_0\times CZ{I}^{0.219},$$

(37)

5. Discussion

5.1. Comparison of Different Models

Various methods have been described to estimate sizes of damage zones around blast holes. In this section, the most commonly used methods are listed and their results are compared through several case studies. For this purpose, 13 rock explosion samples were selected from various studies, in which explosions were carried out in two different types of rocks, including clayey-limestone and basalt, with two types of explosives including ANFO and water-resistant ANFO. The details of these samples including the characteristics of the rocks, explosives, and blast holes are outlined in

Table 5. The characteristics of the rocks, explosives, and blast holes in the 13 studied samples.

Case No.	Rock	Explosive	P (g/cm3)	q (MJ/kg)	${\mathit{D}}_{\mathit{CJ}}$ (km/s)	${\mathit{r}}_0$ (mm)	${\mathit{r}}_{\mathit{c}}$ (mm)	${\mathit{P}}_{\mathit{b}}$ (GPa)
1	CL	ANFO	0.803	3.812	5.016	165	82.5	3.045
2	CL	ANFO	0.803	3.812	5.016	229	114.5	3.477
3	B	ANFO	0.803	3.812	5.016	102	51	2.061
4	B	ANFO	0.803	3.812	5.016	165	82.5	3.148
5	B	ANFO	0.803	3.812	5.016	229	114.5	3.595
6	CL	WR ANFO	0.994	3.918	5.829	51	25.5	2.016
7	CL	WR ANFO	0.994	3.918	5.829	102	51	4.033
8	CL	WR ANFO	0.994	3.918	5.829	165	82.5	4.974
9	CL	WR ANFO	0.994	3.918	5.829	229	114.5	5.44
10	B	WR ANFO	0.994	3.918	5.829	51	25.5	2.085
11	B	WR ANFO	0.994	3.918	5.829	102	51	4.169
12	B	WR ANFO	0.994	3.918	5.829	165	82.5	5.141
13	B	WR ANFO	0.994	3.918	5.829	229	114.5	5.623

.

In the following, five different methods were selected, including Esen’s [2], Il’yushin’s [59], Szuladzinski’s [62], Djordjevic’s [72], and Kanchibotla’s [73] models, to calculate damage size corresponding to each of the 13 samples presented. The results of these models for the introduced case studies are displayed in the columns 2–6 of

Table 6. Results of the selected methods in estimating damage radius for the 13 studied samples.

Case No.	Esen et al. [2]	Il’yushin [59]	Szuladzinski [62]	Djordjevic [72]	Kanchibotla [73]
1	372	1269	379	466	1192
2	564	1761	526	647	1654
3	67	402	108	139	339
4	143	651	175	225	549
5	217	903	242	312	762
6	88	441	132	186	476
7	277	881	264	372	953
8	513	1426	427	602	1541
9	756	1979	593	836	2139
10	34	239	61	90	219
11	107	478	122	179	439
12	198	774	197	290	710
13	291	1074	273	403	985

. The results are also plotted in Figure 7.

As can be observed from Figure 7, in almost all the 13 cases, Il’yushin’s and Kanchibotla’s models yielded relatively larger results than the other models. One possible reason for this issue is that both models used an ideal explosion assumption in their relationships, while the real cases are closer to a non-ideal one. This issue has been pointed out by other researchers as well. It has been argued in the related literature that the purpose of Kanchibotla’s model is to study mine fragmentation and optimize mine excavation [2]. Therefore, it may not be able to provide a precise crack propagation estimation and damage zone determination. The other three models, i.e., Esen’s, Szuladzinski’s, and Djordjevic’s models, correspondingly provide relatively close estimations. However, Esen’s model is preferred to the other models in the literature. The advantage of using Esen’s model compared with other ones is its development based on experimental and real-scale samples, whose results have also been validated by rock blasting real projects. However, the other three models have been developed theoretically through several simplifications and simple assumptions. Moreover, there are several comparisons between these models available in the literature, and almost all of them have addressed Esen’s model as the most accurate approach to predict the crushed zone size. For instance, Amnieh and Bahadori [119,120] conducted several single-hole explosion tests at the Gotvand Olya Dam and compared the results with those of Ash’s, Djordjevic’s, Szuladzinski’s, Kanchibotla’s, and Esen’s models. Eventually, they concluded that the results obtained from Esen’s model were closer to real values compared with the other models in all the cases observed. Additionally, studying the crushed zone radius around a blast hole using laboratory tests, Changshou Sun [6] introduced Esen’s model as the most complete set of data for the crushed zone extent and then applied it to validate their laboratory test results.

5.2. Probabilistic Approaches

The models investigated in this review up to now have been addressing single-hole rock blasting in terms of a deterministic problem. Although the deterministic estimation of damage has been incorporated in a simple algorithm without any certain complexity, some uncertainties have remained unnoticed in the problem. Neither rock medium nor blast load characteristic values are deterministic and consequently cannot be determined with certainty. The deterministic estimation of damage zone is similar to assuming a 100% probability for a specific failure size, which does not appear to be a rational workaround. A better solution in this regard is to match each damage zone radius with a failure probability. In other words, instead of the deterministic estimation of damage zone, an exceedance probability is assumed as the goal of analysis. In fact, determinacy should be left aside and uncertainty needs to be modeled using random variables. As a result, one could estimate the probability required for a crack to exceed a certain length value. This issue has been formulated as a reliability problem in the related literature and investigated in some research works.

Defining the involved parameters as random variables with certain mean and standard deviation and establishing a limit state function for crush and crack zone radii, Shadabfar et al. [121,122,123] incorporated the Monte Carlo method to calculate failure probability. Based on their results, they concluded that exceedance probability was severely reduced following an increase in the crushed zone radius. Accordingly, the probability of the crushed zone radius longer than 0.5 m was reported to be less than 1%. Moreover, the comparison of different probabilistic models developed based on Esen’s, Szuladzinski’s, Djordjevic’s, and Kanchinotla’s models revealed that the results of Esen’s model could exhibit a lower failure probability than those of the other models. In other words, Esen’s model has a more optimistic estimation of the bearing capacity of rock mass compared with the other models.

Additionally, it was observed that Kanchibotla’s model was much more different than the other models, yielding a much higher failure probability. This issue has also been noticed by other researchers [124]. It is reasoned that the objective of Kanchibotla’s model is to study mine fragmentation and optimize mine drilling. Hence, it may not yield a precise crack propagation estimation and damage zone determination.

Furthermore, defining both the crushed and cracked zones in terms of a reliability problem in another study, Shadabfar et al. [125] took advantage of the first-order reliability method and calculated explosion-induced failure probability. The obtained results of this study were then represented in terms of exceedance probability versus damage zone radius (Figure 8).

The diagram was then drawn for all the points in the vicinity of the blast location and presented as exceedance probability contours (Figure 9).

As shown in Figure 9, failure probability is reduced via distancing from the blast location. Thus, failure probability for the crushed zone approaches 0 by exceeding 0.5 m.

In addition, using a parametric study, ShadabFar et al. [125] investigated the effects of decoupling on resulting failures. For this purpose, a number of reliability analyses were conducted and changes in failure probability were recorded accurately through incorporating different ratios of explosive charge to blast hole radii in the governing equation. According to their results, failure probability declined and the exceedance probability diagram showed a lower level of probability as the decoupling ratio was reduced. The reason is that the distance between explosive charge and blast hole wall increased through a reduction in the decoupling ratio, leading to the sharp damping of blast wave before propagating in the rock medium. Additionally, the results revealed that the highest impact of the decoupling ratio occurred in the range of small damage zone radii. More precisely, the decoupling ratio mostly affected failure radii in the range of 300–2350 mm, while its impact was severely reduced for larger failure radii.

Performing a set of reliability sensitivity analyses, ShadabFar et al. [125] also compared the influence of parameters involved in different models. The results of their study showed that the main parameters affecting the crushed zone radius in Esen’s model were the uniaxial compressive strength of rock and the blast hole radius. Furthermore, according to Senuk’s model, two main parameters influencing the cracked zone radius were the blast hole radius and the tensile stress of rock mass. The results of this analysis were then presented in a diagram in terms of sensitivity vector for each parameter involved in different models (Figure 10). This diagram presents the relative importance of each parameter in comparison with other parameters. Load variables (i.e., components whose growth would increase failure probability) are marked with a dark color, whereas resistance variables (i.e., components whose rising trend would decrease failure probability) are marked with a bright color.

6. Conclusions

In this review, the most important existing models for the estimation of explosion-induced crushed and cracked zones were investigated. The models were categorized into three groups, i.e., analytical, numerical, and experimental approaches. First, the rock explosion mechanism was described from the detonation initiation and stress wave propagation to the rock failure. Subsequently, the induced damage was grouped into two forms of crushed and cracked zones and the most important parameters affecting these zones were reported. Then, the most important methods for estimating the dimensions of each damage zone were examined.

More specifically, the analytical models were presented based on the two main parameters of PPV and blast hole pressure. The numerical methods were addressed via the commonly used numerical codes, including the FEM, DEM, and FDM methods. The experimental models were divided into two general categories of primary cracks due to high-stress waves and deep cracks caused by gas penetration and discussed in detail. Finally, a number of empirical models drawn from laboratory results were used in a step-by-step approach.

The most commonly utilized models described in this review were selected to separately calculate the damage dimensions in 13 case studies, which were collected from the related literature and presented in an integrated form. All the results were compared, and their differences or similarities were discussed.

Thereafter, the probabilistic models available for analyzing the failure probability induced by the rock explosion were deliberated, and their advantages over the deterministic models were described. The comparisons were made between different models, and the relative importance of the involved parameters was investigated via the reliability sensitivity analysis.

This review categorized and reported the most important assumptions and key points of the related literature along with their prominent results to allow for more coherent and accurate use of their content. Hence, the results of the present study can be used as a comprehensive and categorized source for estimating the blast-induced damage in rocks. However, for the practical use of the methods presented here, their main sources should be utilized for more details. Finally, this paper only covered the single-hole blasting. In cases with multiple blasting, due to the interaction between the explosion waves released from each blast hole, a system of damage zones is generated, which can potentially overlap and cause more complex failure in the environment. This topic remains as the authors’ concern for future research works.