Simulation of Dynamic Fracturing of Continuum Rock in Open Pit Mining ()
1. Introduction
The term “rock (or blast) fragmentation” is an index that is used to estimate the effect of bench blasting in the mining industry. Knowledge of the fragmentation mechanisms in explosively loaded rock is critical for developing successful methods for excavating rock rapidly for a variety of purposes, and has advanced considerably in the last twenty years. In rock blasting, it is generally understood that both the stress wave and the gas pressurization make significant contributions to rock fragmentation. The importance of shock and gas in fragmentation has been debated for the last 50 years. Recent studies tend to support the view that stress waves generated by the detonation of an explosive charge are responsible for the development of a damage zone in the rock mass, and for the subsequent fragment size distribution, while the explosion gases are important in separating the crack pattern that is formed after the passage of the stress wave, and in throwing the fragments.
Drilling and blasting is the main and most economical procedure to extract valuable mineral resources from the earth. Researchers still make different hypotheses in explaining the fundamental operative mechanisms responsible for rock fracturing by explosive energy in spite of its prevalent use at a large scale. Enormous experimental research efforts have been made over the last five decades to understand the rock fracturing mechanisms.
The process of rock fragmentation by blasting is a complicated phenomenon which is controlled by many variables and parameters. Considering all this parameters in a single analysis is not possible at the present time, especially when some of them are not clearly understood yet and the effect of others is difficult to quantify. In most blasting practices, empirical or semi-empirical techniques are used for blast design and fragmentation analysis. These techniques are based on information obtained for certain range of rock types and blasting conditions and cannot be generalized for all blasting conditions. With regard to the limitations of empirical methods, numerical methods are viable tools to further understand and illustrate the fragmentation process. Application of numerical methods in blasting allows for consideration of complex boundary conditions, material non-linearity, dynamic material behavior, geometric non-linearity and complexity associated with blasting operations. The energy release characteristics in the former are highly variable, depending on the prevailing field parameters such as borehole diameter, density gradient, and sympathetic pressures among explosive columns from detonation of a neighboring hole. Similarly, the response of the target rock to high dynamic loadings, which may last only for a few milliseconds, remains largely unknown. Under this scenario, the best approach to the study is to first generate an extensive experimental database on these properties. At the same time, it is essential to investigate the fracture and fragmentation mechanism through numerical models so as to obtain a better understanding of the dominant parameters that control blast results.
In rock blasting, it is generally agreed that two types of loading operate on the surrounding rock: stress wave (or shock wave) loading and explosion gas pressure loading. The stress wave loading arises out of detonation of the explosive column in the borehole. For typical commercial explosives, the detonation pressure exerted on the borehole wall at the moment of initiation could easily exceed 10 GPa. This high pressure on the borehole wall sets off a shock wave in the adjacent rock mass, but it soon decays to a high amplitude stress wave propagating at the velocity of longitudinal wave in the rock mass. It is immediately followed, albeit at a much reduced velocity, by the longer duration gas pressure loading. This loading due to gas expansion continues long after the stress wave has dissipated, as its expansion rate is considerably lower than that of the propagating stress waves [1-3]. The former initiates cracks around the borehole, and the latter penetrates into these cracks and causes their further extension and propagation. Therefore, both loadings play a very important role in efficient use of the explosive energy. Theoretical studies on rock blasting [3-5] generally agrees with the crack extension data; however, the main focus has been on the propagation of the pre-existing cracks under the gas pressure loading and lesser attention has been paid to the initial cracks resulting from stress wave loading. Numerical methods have been employed by several researchers [6-11], using various numerical codes to simulate fragmentation process in blasting. Although many significant results have been published, it is far from complete for the numerical study of rock fragmentation.
The objective of first part of this numerical study is to analyze rock fracture mechanisms and explain how the crushed zone, severely fractured zone, and incipiently cracked zone are created under stress wave loading. The second part of the study discuses about fault effect on wave propagation. The nature and degree of heterogenity of the rock mass is very important in blast design. That is, discontinuities such as joints, bedding planes, faults, and soft seams can allow the explosive’s energy to be wastefully dissipated rather than perform the work intended. In some cases, the discontinuities can dominate the fracture pattern produced by the explosive, and the influence of the structural geology often overshadows that of the rock’s mechanical and physical properties. Best fragmentation is usually obtained where the face is parallel to the major discontinuity set [12]. The last few decades have seen a variety of numerical studies on the blast-induced waves and their propagation in rock masses with much efforts being placed on the study of dynamic responses of continuous rock masses under blast loading [13,14]. However, rock masses encountered in reality generally contain geological discontinuities (e.g. joints, faults and bedding planes). The properties of rock masses are determined by both the properties of the intact rock and the discontinuities [15,16].
There are some numerical tools available for discontinuous rock masses at present, the popular ones being the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM), etc. For example, joints or faults are modeled as a kind of special joint elements in FEM [17], whereas in FDM joints or faults can be simulated using slide-lines [18]. However, these treatments are inefficient when numerous discontinuities exist. Although the equivalent continuum hypothesis maybe employed at certain circumstances, it is incapable of handling the responses under dynamic loadings [19,22].
The DEM, which was developed to model discontinuous problems, has been recognized to be a more superior alternative for studying the mechanical behaviors of fractured rock masses [21]. The two-dimensional numerical code can simulate the responses of discontinuous media subjected to either static or dynamic loadings. The discontinuum is represented as an assemblage of discrete blocks. The internal discontinuities are treated as boundary conditions between blocks, and large displacements along discontinuities and rotations of blocks are allowed. Individual blocks behave as either rigid or deformable material. Deformable blocks are subdivided into constant strain triangle-shaped finite difference zones, and each zone responds according to a prescribed linear or nonlinear stress-strain law [21].
It is intention of this paper to take advantage of the powerful numerical tool and simulate the significance of the presence of a faulted rock mass subjected to blast loading. The rock mass failure pattern and stress distribution through the rock mass in the presence of the fault is examined.
Numerical codes are useful tools to build models of complex problems, which have complex geometry, loading conditions and boundary conditions. The rock-explosive interaction in faulted rock masses is a perfect example of such problems. Laboratory and field blasting experimentation is very difficult, expensive and not easily doable in the actual field scale. On the other hand, sophisticated codes enable handling of dynamic behavior and can include complex geometries and non-elastic material behaviors. Numerical methods, once calibrated with practical experiments and observations, can be used for parametric studies aimed at analyzing the effect of critical parameters on the structure response.
2. Modelling Strategy and Input Data
As pointed out in the previous sections the goal of this work was to use 2D distinct element code to look into stresses role on fracture initiation and propagation and the effect of stresses on blast-induced wave propagation in surface mining. Figure 1 illustrates an overall view of the model.
The model consists of a single vertical blasthole of 76 mm in diameter with a burden of 5 m. Since the objective of the study was to look at the problem from a stress/ failure mechanism point of view, the Mohr-Coulomb material model was used to model the rock mass behavior. The problem can be treated as a plane-strain case, in which the x-axis and y-axis origin is at the centre of the model. In order to simulate the wave propagation with accuracy a uniform discretization was used. Moreover, to prevent the unwanted boundary reflections, absorbing boundaries were employed in all sides with the exception of upper (northern) boundary which is the free face. The rock type was modeled using limestone material properties. Limestone is the dominant rock type in most highway cuts in northern Iran as well as in the most mines. The materials properties used are given in Table 1.
Natural dynamic systems contain some degree of damping of the vibration energy such that the system would not oscillate indefinitely when subjected to driving forces. The distinct element code uses a dynamic algorithm for the solution of mechanical problems. The damping in the numerical simulation should attempt to reproduce the energy loss in the system when subjected to dynamic loading. Rayleigh damping is commonly used which is approximately frequency-independent over a certain range of frequencies [21]. Damping parameters are very important for DEM dynamic analyses. An important characteristic of DEM dynamic analysis is the reproduction of frequency-independent damping of materials at the correct level. For geo media, natural damp-
Figure 1. Overall plan view of the UDEC single blasthole model.
Table 1. Rock mass properties used as input.
ing is mainly hysteretic and difficult to decide, but it commonly falls in the range 2.0% 5.0% of critical damping [16]. Herein, the fraction of critical damping of 0.02 is specified in the present study for better results.
Blast pressure is exerted on the blasthole’s wall such that a force term, i.e. pressure multiplied by area, is eventually assigned to the rock elements. In response to the force, the rock mass elements expand outward in a radial direction, colliding with the surrounding elements and eventually produce cracks. In order to estimate the generated pressure from the charge, experimental methods can be helpful. The magnitude of shock wave pressure is a function of velocity of detonation, density and charge’s ingredients [12]. The present study adopts the empirical formula proposed by Lopez [22]. Although this relation is very complicated, but the following equation can estimate blast load:
(1)
where, 
is blast pressure (MPa), 
explosive density (gr/cm3) and VD velocity of detonation (m/s). Putting the dynamite properties in the above equation gives:
Gas pressure usually is considered half of the blast pressure, e.g.:
(2)
If the diameter of the explosive is equal to blasthole’s diameter, i.e. coupled detonation, then there is no gap between blasthole and explosive and the related pressure can be calculated as follow:
(3)
where 
is hole radius (mm), b explosive radius (mm), 
specific heat coefficient, and q shape factor of explosive 2 for cylindrical charges and 3 for spherical charges) [23]. Therefore:
On the other hand, applied dynamic pressure on blasthole’s wall is a function of time due to interaction between rock and generated shock wave. Many experimental equations have been presented to calculate this parameter, but presented equations are widely used. According to Starfield’s equation [24], generated dynamic pressure on the wall (P(t)) is a function of rock density
, explosive density
, P-wave velocity
, velocity of detonation (VD) and PW. The following equation gives P(t):
(4)
The explosive density 
is 1.45 (gr/cm3). Therefore for limestone:
Figure 2 shows the graph based on above equation.
As can be seen from Figure 2, the peak pressure is about 820 MPa. It was assumed that the explosive is of a shocky type and delivers most of its energy in the form of stress wave.
Figure 2. Dynamic pressure applying on the blastholes’s wall.
3. Blasting Induced Fracture in Continuum Media
To investigate wave and fracture initiation and propagation in rock masses, a 20 × 10 m block was constructed and a 38 mm radius borehole was located at the center of the model (Figure 1). Figure 3 shows wave propagation in different time steps. It should be noted that in Figure 3 all boundaries are viscous boundary.
Due to the high radial compressive stresses near the borehole, an intense shear stress field develops near the borehole. In the shear failure zone, the rock is severely crushed. The crushed zone is followed by a severely fractured zone. In this zone the crack density is very high and the damage is severe. Beyond the severely fractured zone is the incipiently cracked zone. In this zone, the radial crack density is lower than that in the severely fractured zone. If the upper boundary of model is free (free surface), then the rock mass failure mechanism differs in the free boundary. Figure 4 illustrates the material status at a free face boundary as well as stress wave reflection.
The action of the reflected stress waves from the free surface causes circumferential cracks, which are more commonly known as spalling cracks. It should be noted that the extent of spalling would depend on the dynamic tensile strength of the target rock and the distance traveled by the stress wave. When stress wave collides the
Figure 3. Stress wave front and rock mass status as a function of time.
free face in compression, reflects back into the rock mass in the form of tensile wave and leads to tensile stress concentration at tips of initiated cracks from the borehole and results in further propagation of fractures. Therefore, linking fractures to the free face causes rock fragmentation in area between borehole and free face and formation of overbreak. To further evaluate the issue, some history points (recording points) were located at different points of the model, for example, in the crushed zone, in the severely fractured zone, in the radial fractures zone and in the free face. Recorded parameters are Sxx, Syy, and Sxy. Figure 5 shows the variations of these stresses at a point located in the crushed zone.
As shown in Figure 5, the main stress caused crushed zone around blasthole is Sxx. The maximum value of xx-stress at crushed zone is about 326 MPa in compression. This value of Sxx leads to crushing the rock mass in a radius about 40 cm around the blasthole. On the other hand, Syy and Sxy also play significant roles in crushing the rock mass at crushed zone. The maximum value of yy-stress is about 102 MPa in compression and for xy-stress is about 60 MPa. At the time of 65 µs after detonation, Sxx hits rock elements at the immediate vicinity of the blasthole and begins to compressive failure of rock material. After 200 µs, Syy activates and reaches to its maximum valu. Although at this time the Sxx reaches to its residual value about 30 MPa, but Syy, on the other hand, starts to continue crushing the rock to further radius.
Figure 4. Reflected stress wave front from the free face and generation tensile failure due to wave reflection.
Figure 5. The variation of Sxx, Syy and Sxy at a point located in crushed zone.
Some history points were also located in severely fractured zone. Severely fractured zone contains a number of short cracks initiated immediately after crushed zone. Figure 6 illustrates a view of severely fractured zone and its intense fractures. It can be seen that just a few fractures from severely fractured zone can propagate in longer distances and results in generation of incipiently cracked zone.
The variations of three mentioned stresses is plotted and drawn in Figure 7. As can be observed in Figure 7, the main cause of generation severely fractured zone is yy-stress component with the maximum value of 54.5 MPa in compression state. This magnitude of compressive stress leads to further propagation of fractures from crushed zone in short lengths (maximum 50 cm in this model) and excessive breakage of rock mass in a radius of about 1 meter. Moreover, xx-stress component contributes in rock breakage but with secondary role with the maximum compressive stress about 6.5 MPa and after 0.45 ms reaches to 10 MPa in tension. Shear stress almost does not have determinant role on rock failure at severely fractured zone.
The similar history points were located in incipiently cracked zone and in a point out of breakage zone. Figure 8 illustrates three graphs of Sxx, Syy and Sxy at two points: one located in incipiently cracked fracture (Target 1) and another one in non fragmented zone (Target 2) (shown in Figure 8).
As can be seen from Sxx variations in two target points, the maximum value of xx-stress a two points is almost similar about 25 MPa in compression, but the main difference between two trends is in residual stresses. The
Figure 6. Illustration of crushed and severely fractured zones around blasthole (a magnified view).
Figure 7. Variations of Sxx, Syy and Sxy at severely fractured zone.
residual xx-stress at target point 1 is about 11 Mpa, while this parameter in target point 2 is about 1.9 MPa. The difference between residual stresses at two target points is the main reason for more propagation of fractures to reach to target 2, while low residual stress at areas around target 2 is unable to break the rock mass at the points. This trend is also visible in Syy and Sxy stresses. In the case of Syy, the maximum magnitudes of compressive and tensile stresses at target point 1 are about 5.2 and 8.2 MPa, respectively.
While these values for target point 2 are about 3.8 and 3.2 MPa, respectively. On the other side, the residual stresses at target points 1 and 2 are about 3.6 and 0.3 MPa,
Figure 8. Illustration of variations of Sxx, Syy and Sxy at two target points.
respectively. Although these differences at stresses at two target points are not so sensible, but contributes in stress concentration at cracks tip in associate with Sxx to fail the rock mass in far distances. The xy-stress follows the same trend as two previous stresses. The Sxy reaches to its maximum value of about 6.2 Mpa at target point 1, while the maximum value of Sxy for target point 2 is about 0.3 MPa. The Sxy in target point 1 leads to, once again, shear stress concentration in propagated fractures from severely fractured zone and contributes in further propagation of fractures in far distances in associate with Sxx and Sxy.
4. Summary and Conclusions
The 2D dynamic distinct element commercial code was employed to study of the failure mechanism of rock under blast loading and the effect of a fault on blast wave propagation. In the study, a continuum rock mass was considered as the basic model and blast load was numerically modeled into the model. Wave and fracture initiation and propagation in the rock mass were plotted as the outputs. Moreover, the histories of Sxx, Syy and Sxy were drawn in some graphs in different points of model for evaluating stresses role on rock breakage. The conducted numerical study showed that in three different failure zones of around a blasthole (crushed zone, severely fractured zone and incipiently cracked zone) the role of stresses is different and every individual stress has different trend and role while failure the rock mass. The rock mass was modeled as a medium strength limestone typical of host rock in highway cuts in northern Iran. A Mohr-Coulomb material constitutive law was used to model the rock mass deformation and failure. Important stress and rock failure patterns were observed as critical parameters to illustrate the stresses effect on blasting induced wave propagation.
NOTES