Development of a 3D numerical model for simulating a blast wave propagation system considering the position of the blasting hole and in-situ discontinuities

Blasting operations are one of the most important parts of geotechnical and mining projects. Most rocks naturally have a series of discontinuities that significantly affect their responses to blast waves. In this paper, the propagation of a blast wave in one intact rock and four rocks with different joint conditions are simulated by a 3-dimensional distinct element code. The results showed that the joint in the model acted as a wave barrier and passed part of the waves, absorbed a portion, and reflected the remaining part into the model. In other words, a discontinuity reduces the energy of the wave and causes more wave attenuation. In addition, a shorter distance between the joint and the hole causes slower wave propagation and greater damping. Moreover, the results showed that the smaller the angle between the discontinuity and axis of the blast holes, the more stress occurs in the rock bench.


Introduction
Blasting is one of the most important parts of the geotechnical and mining projects (Hudaverdi, 2012; Kecojevic and Radomsky, 2005; Monjezi and Rezaei, 2011). Blasting in holes drilled in rocks are common in underground and open pit mining, construction of roads, tunnels, wells, and other underground spaces. The quality of a rock mass has a strong influence on the level of dynamic forces to which the rock can be subjected and sustain damage. Most rock masses naturally have a series of discontinuities, fissures, joints, cracks, bedding planes, layering, foliation and faults that significantly affect their physical and mechanical properties, and consequently their behaviour under blast waves (Bakar et al., 2013;Jaeger et al., 2009). It has been reported that blasting in a homogenous isotropic medium naturally results in a different fragmentation pattern in comparison to when the medium is permeated with discontinuities. In most rock materials, fissures occur, thus reducing the explosive induced stresses due to shock wave reflections (Persson et al., 1994). Discontinuities may be closed, open, or filled, so they can pass deferent fractions of explosion energy. The faces of the discontinuities are the surfaces onto which the strain wave is reflected, damped, or propagated (Priest, 2012).
The blasting process in rocks is so sophisticated that it is not yet well known due to certain limitations. Many researchers have studied the theories governing wave propagation. One of the first attempts for presenting an analytical explanation of blasting was made by studying the explosive loading in a hole near the free surface based on reflection theory. This simple and easy theory is based on the fact that the tensile strength of rock is lower than its compressive strength (Yilmaz and Unlu, 2013). De Lemos employed a 2D distinct element model to study the attenuation of a transient transverse wave passing across a single joint with the Coulomb slidingfriction behaviour. The effects of the fracture on wave attenuation are determined as magnitudes of the reflection and transmission coefficients (de Lemos, 1987). Brady et al. showed that the results of numerical modelling based on the distinct element method on the slip of a fracture due to an explosive line source are in agreement with analytical solutions (Brady et al., 1990). Pyrak et al. studied wave propagation in a single joint and a joint set. They found that rock joints act as high frequency filters that only the low frequency waves can pass through (Pyrak-Nolte et al., 1990). Kulatilake et al. used the distinct element method to perform a 3D numerical study on the effect of intermittent joints on the deformation of a jointed rock mass (Kulatilake et al., 1993). Chen verified the capability of the distinct element code to model the behaviour of jointed rock masses in an explosion operation (S. G. Chen, 1999 & Katsabanis developed a dynamic blasthole expansion model which considered the effects of blast geometry (blasthole shape, angle, and location), the physical properties of the intact rock and existing discontinuities, the distribution and orientation of pre-existing discontinuities, and the blasthole pressure on the processes of burden breakage, fragment throw and muck pile formation (Mortazavi and Katsabanis, 2000). Singh reports that the underground controlled blasting techniques, which work well in the massive deposits, bring poor results while adopted in jointed rock formations, especially in controlling over break (Singh, 2005). Zhou  As reviewed above, numerical studies on the dynamic responses of jointed rock masses are mostly performed using 2D models by the distinct element method. In this paper, the propagation of a blast wave in a jointed rock mass is simulated by a 3-dimensional distinct element code, and the position of the joints relative to the blast hole on the blast wave propagation is studied. For this aim, five different models with different conditions consist of: one intact rock without a joint, two models with a single joint, and two models with two joints are studied.

Methodology
The interaction of a blast wave and a discontinuity is generally studied in three ways: experimentally, analytically, and with numerical methods. The numerical methods in comparison with laboratory and experimental studies provide easier and more economical conditions for studying the wave propagations in a jointed rock mass, especially for complex cases where experimental solutions seem impossible. The numerical methods used in rock medias are classified into three categories: continuous, discontinuous, and hybrid methods. The hybrid methods are a combination of the continuous and discontinuous methods. The hybrid methods have been de-veloped in recent years due to the limitations of both the continuous and discontinuous methods. The hybrid methods include: the Numerical Manifold Method (NMM), which is a combination of Discontinuous Deformation Analysis (DDA), and the Finite Element Method (FEM); the Finite Element/ Distinct Element Method; and the Particle Manifold Method (PMM) that was developed by combining the Distinct Lattice Spring Model (DLSM) and NMM methods. The most important distinction and key difficulty of different numerical modelling of wave propagation in jointed rocks is how the joints are presented. In the Distinct Element Method, the rock and joints are introduced as separate blocks and interfaces between blocks, respectively (Zhou and Zhao, 2011). In this paper, the DEM is implemented for analysis of the propagation of a blast wave in a jointed rock mass. Distinct Element Methods (DEMs) use an explicit scheme to evolve the equations of motion of discrete bodies directly. The bodies may be rigid or deformable (by subdivision into elements). When using the DEM, contacts are always deformable and can employ

Three-dimensional Numerical analysis
The main steps required to perform a numerical analysis are: build model geometry, mesh generation, assign constitutive models and then material properties, apply initial and boundary conditions, estimate the initial equilibrium, apply changes to the model, solve the model, evaluate the results and, if necessary, adjust the model. In addition, to solve the dynamic problem, some other steps must be considered, such as evaluation of the transmission wave quality, assign mechanical damping, define dynamic boundary conditions, and apply dynamic loading (Raffaldi and Loken, 2016).
In this research, five different conditions have been studied: one model with intact rock containing two blast holes, and four models with rock containing one joint and two blast holes. The geometry of all four jointed rock models is the same and they differ only in the way the joint is located relative to the axis of the holes. The geometry of the model is contained within a 30×26×20 m 3 rock block. The model has a bench with a 5m burden and a slope of 85 degrees in x direction. Two blast holes with a 6m distance from each other are located in the bench with a 165mm diameter and 15m height. The stemming of blast holes is 4m, and so the length of the explosive charge is 11m. The main explosive used in several mines is ANFO with dynamite as primer. So in this research, the characteristics of the dynamic wave of blasting with ANFO have been used in the numerical analyses. The parameters of the blasting including: hole length, stemming, burden, and hole spacing, are determined based on the theory of blasting in rocks presented in (Lopez Jimeno et al., 1995). The geometry of the numerical model and the details of the blast holes are shown in Figure 1. The maximum mesh size is considered 50cm, which is less than 0.1 of the wavelength emitted due to blasting using ANFO as suggested by Kuhlemeyer and Lysmer (Kuhlemeyer and Lysmer, 1973).
The rock in the first model is assumed intact without any discontinuity. The next four models each have a joint, but differ in the location and direction of the joint relative to the hole. In the second model, one vertical joint is located between two holes in equal distance, and parallel to the blast holes (see Figure 2a). In third model, the joint is parallel to the bench face and is located in the middle of the burden. In other words, the distance of the joint from both the face and the two blast holes is 2.5 m (see Figure 2b). The joint in model #4 is perpendicular to the bench face and is located beside the left blast hole (see Figure 2c). In model #5, the distance of the joint from the blast holes is similar to model #4, but the joint and the bench face make a 45-degree angle (see Figure 2d). The dip of joint in all model is 90 degrees. The three-dimensional view of the situation of the joint in models 2 to 5 is shown in the Figure 2.
The mechanical behaviour of the rock and joint is adopted to be elastic-plastic according to the Mohr-Coulomb criterion. According to this criterion, when a sample of rock is subjected to stresses, if the applied shear stress in the given plane exceeds the shear strength of the rock, failure will occur. The properties of the rock and joint are presented in Table 1 and Table 2, respectively. The rock used in this paper is limestone and the properties of the joint are considered clay.
The model is confined by the appropriate boundaries along the x, y and z in static or dynamic conditions. In dynamic analysis, the load of the blasting is applied as a load history. This means that the applied load is specified as a time function F(t). The blasting load is generally defined as simulated time histories or mathematical functions by an ascending graph that rapidly reaches a maximum pressure at peak time and then slowly decreases to zero. In this research, the blasting load is applied on the wall of the blast holes as a compression wave that is evaluated based on the explosion function proposed by Starfield  (1) (2) Where: t -time; t r -the rise time (time of peak pressure); P P -the peak pressure at the blasthole wall. Figure 3 shows the applied explosion pulse shape which resulted from the explosion Starfield and Pugliese function. It also shows that in this study, the two blast waves propagate without a time difference and at the same time.

Results and Discussion
After numerical analysis on the aforementioned five models, the propagation of the blast wave in the rock and also the effect of the situation of the joint in relation to the blast holes are investigated.

Comparison blast wave propagation in the intact and jointed rock
The propagation of the blast wave in intact rock (model #1) from 18 to 800 microseconds in the longitudinal and transverse sections are shown in Figure 4 and Figure 5, respectively. As seen from 18 microseconds onwards, the velocity contours start to spread in the model in the form of concentric circles in the horizontal section (see Figure 4), and in the vertical section in the form of the cylinders (see Figure 5). The time of the wave propagation is exactly parallel to the timing of the dynamic loading as a compression pulse with a maximum pressure of 1.8GPa (see Figure 3). The wave front becomes wider and larger over time until it reaches the model boundaries. In the viscos boundaries of the model, a portion of the wave is absorbed to prevent the waves from reflecting into the model (see Figure 5g, 5h, and 5i). Furthermore, when the wave reaches the free surface, part of it is transmitted and another part of the wave is reflected as a tensile wave into the model, which is consistent with the theory of blast wave reflection (Kleine et al., 1903).
The propagation of the blast wave in jointed rock of model #3 from 450 to 750 microseconds is shown in Figure 6. It is clear in Figure 6 that the pressure wave propagated in model #3, in the right, left and top boundaries of the model is absorbed due to the presence of absorbent boundaries. In other words, boundaries assigned to the model prevented the emission of the reflected waves into the model. However, the lower boundary (free surface) reflects the wave into the rock mass and the compressive stresses are converted to tensile stresses. Moreover, according to Figure 6, the joint in the model acts as a wave barrier and passes part of the waves, absorbs a portion, and reflects the remaining part into the model. In other words, each discontinuity reduces the energy of the wave and causes more wave attenuation. If the energy of the reflected wave is strong enough (more than the dynamic resistance of the rock), it may cause fractures or even scaling of the surface of the rock fragments.

Effect of joint situation on the wave propagation in jointed rock
The blast wave propagation in models #2, #4, #5 are presented in Figure 7.
Considering the fact that the diagrams of Figure 7 show the same time of wave propagation, it can be concluded that wave propagation in the fourth model is slower than in model #2 and model #5. The discontinuity in the fourth model is next to one of the holes, and so, it prevents the proper symmetrical propagation of the wave front. The results presented in Figure 7 prove that a shorter distance between the joint and the hole causes slower wave propagation and greater damping. A diagram of the peak particle velocity of models #2, #4, #5 is shown in Figure 8.
According to Figure 8, the further the discontinuity from the axis of the hole, the slower the wave damping, and the discontinuities that are close to the hole, such as the second and third models, accelerate the damping of the wave front. Figure 8 consists of three parts: one part, the negative distance to the centre of the hole (shows the rock behind the blast hole), one part in the middle of the graph, and the wave is propagated in the rocky environment (0<distance<5, in front of the hole), and in the final part, the wave has spread in the free surface. For both behind and in front of the hole, the PPV of rock media decreases with an increase in the distance from the axis of the hole.
It is understood from Figure 7 that in model #4, for which the distance between the joint and the stemming is very short, when the explosion occurred, a high stress spread in the model, while in model #3, where the joint is relatively far from the top and middle of the model, less stress is applied to the rock and the blast wave damped later. Figure 9 shows the maximum horizontal stress (S XX ) in a point between the blast's holes. The results showed that the smaller the angle between the discontinuity and the axis of the blast holes, the more stress occurs in the rock bench.

Conclusion
One intact rock and four rocks with different joint conditions are analysed by a 3D DEM model. The propagation of the blast wave in the rock and also the effect of the situation of the joint in relation to the blast holes on wave propagation are investigated. The results showed that the joint in the model acts as a wave barrier and passes part of the waves, absorbs a portion, and reflects the remaining part into the model. In other words, each discontinuity reduces the energy of the wave and causes more wave attenuation. In addition, a shorter distance between the joint and the hole causes slower wave propagation and greater damping. The outcomes also showed that for both behind and in front of the blast hole, the peak particle velocity (PPV) of rock media decreases with an increase in the distance from the axis of the blast hole. Moreover, the results showed that the smaller the angle between the discontinuity and the axis of the blast holes, the more stress occurs in the rock bench. It should be noted that due to the lack of access to information of a real project, in this research, only numerical modelling has been done and field experiments are needed. It is better to combine these results with field validation to achieve a more comprehensive method.