Numerical study of the fracturing mechanism around a blasthole and investigating the effect of discontinuities on the fracture pattern

The mechanism of rock fragmentation around blastholes is of prior importance in the evaluation of drilling and blasting performance in open pit and underground mines. This paper aims to numerically investigate the crack propagation mechanism around a single blasthole using the distinct element method (DEM). In this study, two specimens with different borehole diameters were considered and the effects of the stress waves on their cracking mechanism were simulated. To validate the numerical model, the length of the radial cracks around each blasthole was measured and compared against an analytical fracture mechanics model. The fractured zones around the blasthole were also compared against previous experimental tests and good agreement was observed. The effect of a single discontinuity on the crack propagation mechanism was also studied and it was found that the discontinuity normal stiffness plays a significant role in the fractured zones around the blasthole. For low values of normal stiffness, the discontinuity surface acted as a free surface, and the shock wave was significantly reflected, while at high values of normal stiffness, cracks propagate across the discontinuity surface.


Introduction
Blasting is a common practice for rock mass breakage which is commonly used in mining and civil industries such as open-pit mines and underground tunnels. Although blasting is an efficient technique for rock breakage, it has several drawbacks, such as over-breakage, fly rock, and vibration-induced instabilities (Boris et al., 2019; Gharehgheshlagh and Alipour, 2020). The size of blast-induced fragments is an important parameter in the evaluation of the efficiency of rock blasting in openpit mines. To minimize the noted drawbacks and reach an optimal rock breakage, it is essential to improve our understanding regarding the sophisticated and complex physical and mechanical dynamic breakage of rock mass in the blasting process (Bakhshandeh Amnieh and Bahadori, 2017; Ma and An, 2008; Starfield and Pugliese, 1968; Zhu et al., 2007).
In the blasting process, two types of stress waves and explosion gas loading areas are applied to the surrounding rock mass. The stress wave and gas pressurization theories have been employed (either separately or in a combined form) to study the complex rock fragmentation mechanism (Božić, 1998 However, a comprehensive model which can consider all aspects of the rock blasting process is not available at the moment. Fracturing patterns, crack length estimation, and the number of radial cracks around the blastholes are the main concern in the rock fragmentation analysis. Several experimental, analytical, and numerical studies have been devoted to explore the radial crack propagation mechanism around the blastholes drilled in some typical rocks. Some experimental studies have been performed in both laboratory scale The preparation of a sample that is representative of the rock mass is a challenging task in laboratory scale experiments. In field-scale experiments, special difficulties exist in the interpretation of the results due to the influence of several uncontrolled parameters on the results. Some efforts have been undertaken to solve explosion-based problems using various analytical solutions (Cho, 2003;Ouchterlony, 1983 Rocks are frequently accompanied by joints and discontinuities. During the dynamic fracturing of rocks by shock waves, the stress wave travels into the rock mass and encounters joints. Therefore, partial transmission and reflection of the encountering stress wave may occur while the stress wave crosses the joint plane. The joint properties will definitely affect the transmission and reflection characteristics of the stress wave as well as intact rock fragmentation in both sides of the rock defects A review of previous studies clearly shows that it is impossible to fully quantify the simultaneous effects of a shock wave and explosion gas flow in the jointed rock mass based on our current knowledge. As the shock wave development plays an important role in the initial phase of the rock fragmentation, it would be helpful to study the behavior of a rock mass exposed to shock wave stresses. This paper aims to investigate the effect of a blastinduced stress wave on rock fragmentation in intact rock and jointed rock mass. To this end, two different cylindrical samples with a central borehole were numerically simulated using a universal distinct element code (UDEC) and the results were compared against physical experiments undertaken by Dehghan Banadaki (2012). The results of the numerical models were compared against analytical models based on the fracture mechanics theory. Then, using the validated numerical models, the effects of mechanical and geometrical properties of single discontinuity on rock fragmentation and the fracturing pattern were investigated.

Rock fragmentation analysis using the linear elastic fracture mechanics (LEFM)
Understanding the concepts of linear elastic materials and the mechanism of rock fracturing has been a key element in the solution of many engineering problems related to rock structures (Broek, 1989;Whittaker et al., 1992). The mechanism of crack propagation in a rock mass is a complex problem and often needs some modern theoretical and numerical methods for the prediction of fracture patterns. Fracture initiation is the beginning of the crack extension process (Broek, 1989;Hoek and Bieniawski, 1965). In the analytical and theoretical rock fracture mechanics, it is common to consider rock as an ideal material i.e. rock is considered as homogeneous, isotropic, and a linear elastic material. This asumption is based on Irwin's suggestion of the critical strain energy release rate which is found to be useful to study the mechanism of the rock fracturing process (Bažant, 1984;Feng, 2017;Hoek and Bieniawski, 1965;Whittaker et al., 1992). In the linear elastic fracture mechanics (LEFM) analysis of cracks, the stress intensity factor K at the crack tip plays the most important role and it is evaluated by assuming that Hooke's law is valid without limitations to stresses and strains. However, stress intensity factors are defined based on stresses rather than loads in engineering analyses (Broek, 1989;Whittaker et al., 1992).
Irwin originally defined the concept of stress intensity factor K, as a measure of the stress of singularity at the crack tip. He explained that controls local stresses, where a is half of the crack length which is subjected to a stress field of σ (Irwin, 1957).
Considering a rock material of arbitrary shape with a crack of arbitrary size subjected to a specified loading condition as shown in Figure 1  (1) Where: r and θ -polar coordinates used to describe a specified point from the crack tip, K I and K II -stress intensity at the tip of a crack in tensile (mode I) and shear (mode II) modes, respectively. Based on the concept of LEFM, when the stress intensity K reaches a critical value, called fracture toughness K IC , the crack propagates. Conceptually, the rock fracture toughness represents the material resistance to crack propagation under a given stress state (Whittaker et al., 1992).
Rock material can be divided into two or more segments under the applied stresses, which are called rock fragmentations. The strength of a rock structure can be considerably reduced due to the existence of a crack in a rock material. The applied stress at a rock body may in-crease several times at the tip of a crack (discontinuity) due to the stress concentration phenomenon. When this crack tip stress reaches a critical level (i.e. the fracture toughness of the rock), the crack propagates ahead of its tip in a particular direction (Broek, 1989; Hoek and Bieniawski, 1965).
The LEFM theory has been effectively employed in dynamic rock fracture mechanics as well as the crack branching process, which occurs in some dynamic phenomena such as rock blasting and rock fragmentation by disc cutters (Bažant, 1984;Feng, 2017;Irwin, 1957). In this study, the LEFM theory is used to assess the dynamic crack propagation process around a single blasthole.

Analytical Solution for rock blasting
The explosive charge detonation in blastholes generates shock waves in the rock mass. These shock waves may insert high compressive stress far beyond that of the strength of the rock material at the vicinity of the blasthole in a radial direction. The radial expansion of the rock just behind the wavefront relaxes the stresses. Then, the produced intermediate tangential stress component soon tends to that of the tensile which initiates the radial cracks around the blast hole (Bažant, 1984;Marji et al., 2006). The generation of a radial crack system in PMMA (Poly-Methyl MethAcrylate) for a model scale blast is illustrated in Figure 2.
A circular blasthole with radial cracks of alternating lengths is shown in Figure 3, which resembles the blasthole situation. The star crack with arms of alternating lengths resembles the situation at a later stage when the cracks have grown large compared to the hole radius (Rossmanith, 2014). The expansion load configuration was analyzed by Ouchterlony (Ouchterlony, 1983). Due to the symmetry of the model, these cracks will most likely grow in the radial direction. The normalized crack of length m can be represented by Equation 2: (2) Where: a -the length of a radial crack having n arms, R -the radius of blasthole.
Where: σ x , σ y and τ xy -the stress state, u x and u y -displacements in x and y directions, μ -the shear modulus, ν -the Poisson's ratio, k = 3-4ν in the plane strain condition. Based on the Ouchterlony's analytical solution explained above, Figure 4 is adopted for the radial crack propagation around a typical blasthole in a rock mass (Ouchterlony, 1983).
Stresses at the crack tip and the corresponding radial crack growth around a blasthole can be studied by using a normalized crack of length m. When the radial cracks become comparatively long i.e. m>>1, the pressure acting on the blasthole between the two adjacent radial cracks can be expressed concerning the resultant force, F as in Equation 5 (Ouchterlony, 1983): Where: p -pressure around the blasthole. Therefore, it is possible to deduce Equation 6 for the Mode I (tensile) stress intensity factor, K I considering the radial cracks of equal lengths.

(6)
Equation 6 shows that the value of K 2 (n) decreases with an increase in the number of radial cracks n as shown in Figure 5.

Methodology
In this research, the universal distinct element code (UDEC) code was employed to simulate the dynamic crack propagation mechanism around a blast hole. Numerical models with similar geometry of experimental tests undertaken by Dehghan Banadaki (Dehghan-Banadaki, 2010) were generated. In these experiments, copper tubes were installed inside the blastholes to remove the penetration of explosive gases inside the rock  Table 1. Due to the similarity of geometry along the vertical axis of rock samples, a two-dimensional plane strain condition was used in this research. Two samples having borehole diameters of 9.65 and 6.45 mm, copper tube thicknesses of 0.8 and 0.6 mm and coupling air media of 3.98 and 2.38 mm, similar to experimental tests (Dehghan-Banadaki, 2010) were simulated in UDEC. The geometry of these samples is shown in Figure 6.

Where:
Δl -the mesh size (or element size) in the numerical modeling, f -the input frequency, λ -the wavelength, C p -the speed of p-wave propagation in rock media, ρ -the material density, K and G -bulk and shear modulus, respectively. Based on the above equations and preliminary studies, the element size of 2 mm was found appropriate for this study.

Simulation of the Dynamic load
There is no direct way to measure the blast pressure behind the blasthole wall during the explosion process. In this study, the empirical relation suggested by National Highway Institute (1991) was used to simulate the blast pressure. Based on this empirical model, the maximum wall pressure P m (kbar) can be estimated as Equation 9 (Konya and Walter, 1991): Where: ρ e -the explosive density (g/cm 3 ), VOD -the velocity of detonation (ft/sec). This pressure is the pressure of the fully coupled explosive (when the explosive completely fills the blasthole). Since the blasthole may not be fully filled with the explosive, there might be an empty space between the explosive and the blasthole wall which can reduce the blasting pressure. Niv and Olsson proposed Equation 10 for the calculation of the effect of decoupling on the blasthole wall (Nie and Olsson, 2001): Where: P b -the modified pressure behind the blasthole wall, r c -the coupling ratio (the ratio of the explosive diameter to the blasthole diameter), λ -the adiabatic constant of the gas explosion, which can be approximately set to 1.2. Therefore, the maximum modified pressures on the blasthole wall for the first and the second samples can be estimated as 200 and 527MPa, respectively.
The dynamic pressure acting on the blasthole wall varies over time. Therefore, the dynamic pressure-time history is required to simulate the time domain behavior of the rock material. The blasting pressure-time function introduced by Starfield and Pugliese (Starfield and Pugliese, 1968)  The size of the numerical elements is of great importance in the simulation of blast-induced wave propagation in the rock medium. As suggested by Kuhlmeyer and Lysmer,the element size Δl in the DEM modeling should be lower than one-tenth to one-eighth of the wavelength λ. The following relations (Equations 7 and 8) were used to find the appropriate element size Where: t r -the rise time.
Considering all above equations and the proposed approach, time histories of the wall pressure for two samples of this study can be obtained as Equations 13 and 14: For the first sample (13) For the second sample (14)

Numerical simulation of radial cracks propagation around a blasthole and validation study of the numerical model
The fracturing mechanisms for two numerical models were recorded by numerical modeling and are compared against physical experiments in Figure 7 and Figure 8. It can be seen that three fractured zones are generated in the samples. In the vicinity of the blasthole, a crushed zone is generated where the rock material is extremely damaged. With a distance from the crushed zone, a ra- dial crack zone is generated, which is represented by long cracks. The third zone is the spalling zone where the stress wave strikes the free outer boundary surface and is reflected as tensile stress into the medium. A comparison of numerical and experimental fractured zones in Figure 7 and Figure 8 shows that the numerical model can reproduce failure zones around the blast hole very well. Some asymmetric radial crack propagation can be observed in experimental fractured zones of the first sample which may refer to the heterogeneity of the first sample. However, this asymmetrical crack propagation cannot be observed in the numerical model due to the homogenous nature of the numerical medium. To further investigate the validity of the numerical model, the lengths of radial cracks in numerical models were compared against the analytical model explained in Section 2. According to Eq. 6, the crack length can be obtained based on the fracture mechanics theory as follows: The analytical crack length for different samples was calculated and compared against the average numerical crack length. The results are presented in Table 2. These results clearly show that there is a good agreement between the results of numerical models against experimental and analytical results. Therefore, the numerical model provides a valid approach for studying the frac-ture mechanism of rock materials. In the following section, the effect of rock discontinuities on the fracturing mechanism will be investigated.

Effect of rock discontinuities
on the rock fragmentation

Effect of mechanical properties of rock discontinuities
In rock engineering, the interaction between blast wave and rock discontinuities is crucial in the rock fragmentation process. When the stress wave encounters a rock joint, some part of the stress wave is transmitted across the joint plane and, some part of the stress wave is reflected in the medium. Therefore, the rock joints can significantly alter the fracture pattern around a blasthole Two cases were considered in this study to investigate the effect of a joint's normal stiffness on the crack pattern around a blasthole for both blast tests noted above. For case A, the joint's normal stiffness was set at 50 MPa/mm while for case B the joint's normal stiffness was changed to 5000 MPa/mm. The ratio of joint's normal stiffness to shear stiffness was set 0.6 and the cohesion and friction angle of both joint cases were set to 3 MPa and 35 degrees, respectively. The fracture patterns of both samples for these cases are illustrated in Figure  9 and Figure 10.
At the beginning, the blasting induced radial cracks which are initiated from the blasthole and gradually propagate towards the sample edge. The induced stress waves travel into the rock and strike the joint plane. For low normal stiffness joints, the fracture acts as a free face so that the tensile waves will be reflected from the joint plane (as shown in Figure 9(a) and Figure 10(a) for Joint-1). This reflected tensile stress is usually higher than the tensile strength of the rock; therefore, rock spalling occurs gradually parallel to the joint surface towards the blasthole. On the other hand, for high normal stiffness joints (as shown in Figure 9(b) and Figure  10(b) for Joint-2), the induced stress waves cross the joint plane (i.e. there is no reflection and no rock spalling phenomenon in this case).
The joint's normal stiffness can be representative of a joint's normal deformability or joint aperture. For case A, the value of a joint's normal stiffness is low which means that greater normal deformation is required to transmit the applied stress wave. So, it reacts like an open joint and the stress wave is mainly reflected like encountering a free surface. However, for case B the stress wave is transmitted with lower normal deformation and the joint behaves like a closed and stiff fracture. Therefore, the radial crack propagates across the fracture.

Effect of Joint orientation on the fragmentation pattern
The effect of joint orientation was studied for two samples by changing the dip angle from 15 degrees to 45 degrees as shown in Figure 11. The mechanical properties of rock joints were chosen similar to the previous section having joints with a normal stiffness of 5000 MPa/mm. The fracture patterns of both samples are depicted in Figure 12 and Figure 13.
These results clearly show that as the blast wave propagates radially from the blast hole, the vertical distance from the blasthole is a controlling parameter. For the case of 15 degrees, the joint is at the closest distance from the blasthole and significant stress wave reflection and stress transmission can be observed. This results in the generation of spalling cracks parallel to the joint surface as well as propagation of radial cracks. The increase of the joint dip angle results in an increase of distance from the blasthole. Due to the decay of the stress wave with propagation in rock material as well as a high joint›s normal stiffness, the amount of reflected stress wave reduces, less tensile spalling cracks are observed and the fragmentation is mainly caused by radial cracking.

Conclusions
In this paper, the effect of stress wave loading on crack initiation and propagation was numerically studied using distinct element modeling. The nature of all subsequent crack extension, branching, and coalescence would be largely governed by the initial crack patterns generated by the stress wave loading. Two samples with different geometry and blasting conditions were numerically simulated and the corresponding fracturing pattern was compared against physical and analytical models. Three zones of the highly damaged zone, radial cracking and spalling were observed in numerical models and good agreement was found between numerical and physical experiments. The lengths of numerical radial cracks were also found in good agreement with the analytical model based on the fracture mechanics. The effects of joint mechanical and geometrical parameters on fracture patterns were also investigated and it was found that two parameters of a joint's normal stiffness and vertical distance from the blasthole are significant parameters on the fracture pattern. It was found that at low values of a joint's normal stiffness, joints react as free surfaces and significantly reflect the stress waves, while for high values of a joint's normal stiffness, the stress waves mainly transmitted across the joint planes.