- Email: [email protected]
Numerical simulation of blasting-induced rock fractures
ARTICLE IN PRESS
International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975 www.elsevier.com/locate/ijrmms
Numerical simulation of blasting-induced rock fractures G.W. Ma, X.M. An School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore Received 3 January 2007; received in revised form 3 December 2007; accepted 4 December 2007 Available online 20 February 2008
Abstract In the present study, the Johnson–Holmquist (J–H) material model is implemented into the commercial software LS-DYNA through user-subroutines to simulate the blasting-induced rock fractures. The J–H model consists of strength models for both intact and fully fractured materials, a polynomial equation of state, and a damage model that represents the material from an intact state to a fully fractured state. Influences of the key parameters in smooth blasting, viz., loading rate, distance from a free face, earth stress, and preexisting joint planes, etc., on fracture patterns are explored. According to the simulation results, the rock fracture pattern is significantly influenced by the loading rate. Fracture control techniques, namely, notched borehole and charge holder with slits are also simulated. Effectiveness of the fracture control techniques is demonstrated. The numerical simulation in the present study reproduces some of the well-known phenomena observed by other researchers. It has the potential to be applied in practical blast control and gas and hydraulic fracturing engineering. r 2007 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; LS-DYNA; Johnson–Holmquist material model; Rock fracture; Blast control
1. Introduction Drill and blast method has been widely used in underground excavation and construction for many years and it is still a popular method of rock breakage. However, a clear description of the fracture process investigated in the rock is still lacking. When an explosive is detonated, a chemical reaction occurs very rapidly and a relatively small quantity of explosive is converted into gas of very high temperature. This reaction results in two types of loadings applied on the borehole wall, namely a stress wave and a gas pressure with longer duration [1,2]. Rock fractures will be subsequently initiated and propagated in the surrounding rock mass. In a coupled borehole blasting, a crushing zone may be generated at a specific depth from the borehole surface. Recent studies have revealed that the stress wave is responsible for initiation of the crushing zone and the surrounding radial fractures, while the gas pressure further extends the fractures [2].
Corresponding author. Tel.: +65 6790 4984; fax: +65 6790 6841.
E-mail address: [email protected] (G.W. Ma). 1365-1609/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2007.12.002
Investigation of rock blasting and its effects by scaled or full size experiment is very expensive and time-consuming. On the other hand, numerical method, derived from sound mechanical principles and validated against experimental data, indicates a promising approach to reveal the blastinginduced rock fracture process. Considerable efforts have been directed towards developing discontinuum [2,3–7], continuum [8–18], and continuum–discontinuum [19] constitutive models of the rock materials. For a discontinuum model, the rock mass is represented by an assembly of jointed elements (e.g. particles and springs) and a fracture is initiated when the stress in the joint between two elements exceeds a critical value. However, it requires large amount of data about the spatial distribution, geometrical and mechanical properties of the rock joints, which are normally not readily available. Continuum damage mechanics is also considered appropriate to describe the blast damage process from a stress wave. The continuum damage mechanics is based on the fact that the brittle rock has pre-existing defects. The material failure is a process of growth and nucleation of these defects. Material strength depends on the strength and distribution of the critical defects. A single variable, or
ARTICLE IN PRESS 967
Intact (D=0)
Damaged (0
HEL
Fully Fractured (D=1.0)
T*
Normalized pressure, P*
D=1
Pressure, P
a set of variables, can be used to describe the changing state of a material. For a continuum damage model, it is critical to define an appropriate damage evolutionary law. Johnson and Holmquist [15,16,20] suggested a constitutive mode for brittle materials (J–H model) that has been applied to high-purity ceramics. The J–H model consists of strength models for both intact and fully fractured materials, a polynomial equation of state, and a damage model that represents the material from an intact state to a fully fractured state. It is pressure- and strain-ratedependent and suitable for large strain, high strain rate and high pressure problems. In the present study, the J–H model is implemented into a commercial software LS-DYNA [21,22], which is an explicit, three-dimensional, dynamic, nonlinear finite element program, through its user-subroutines. Blastinginduced rock fractures are simulated for different loading and boundary conditions of a borehole. It is found that the numerical simulation in the present study reproduces some of the well-known phenomena observed by other researchers. Fracture control techniques, namely, notched borehole and charge holder with slits are also simulated. Effectiveness of the fracture control techniques is clearly demonstrated.
Normalized Equivalent Stress, σ∗
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
0
The J–H material model is schematically illustrated in Fig. 1. The model gives the intact and the totally fractured strengths (Fig. 1a), a polynomial equation of state (Fig. 1b), and a damage model that represents the material from an intact state to a fractured state (Fig. 1c). The intact strength (upper solid line in Fig. 1a) is given as
Relative Volume, μ
si ¼ AðP þ T ÞN ð1 þ C ln _ Þ.
Equivalent Plastic Fracture Strain, εpf
2. J–H material model
(1)
The fully fractured strength (lower solid line in Fig. 1a) is given as sf ¼ BðP ÞM ð1 þ C ln _ Þ.
(2)
The strength of damaged material (dashed line in Fig. 1a) is expressed as s ¼ si Dðsi sf Þ,
(3)
where A, N, B, M, and C are material constants; si , sf , s, P, and T are effective stress of the intact material, effective stress of the fully fractured material, current effective stress, pressure, and tensile strength, respectively, which are normalized by the following expressions: si ¼ si =sHEL , sf ¼ sf =sHEL , s ¼ s=sHEL , P ¼ P=PHEL , and T ¼ T=PHEL , in which P is the current hydrostatic pressure, T is the maximum tensile hydrostatic pressure the material can withstand, sHEL is the effective stress at the Hugoniot elastic limit (HEL) which is defined in the original J–H model to characterize material yielding in a shock wave condition, PHEL is the associated hydrostatic pressure at the HEL, D is the damage variable, s is the
T∗
Normalized pressure P∗
Fig. 1. Schematic illustration of J–H material model: (a) strength model, (b) pressure–volume relationship, and (c) damage model.
effective stress which has the general form of s¼ (
ðsx sy Þ2 þ ðsx sz Þ2 þ ðsy sz Þ2 þ 6ðt2xy þ t2xz þ t2yz Þ 2
)1=2 ,
ð4Þ where sx, sy, and sz are the three normal stresses, and txy, tyz, and tzx are the three shear stresses. It should be mentioned that rate dependence of the J–H model is represented by a dynamic strength increase factor
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
968
(DIF) as follows,
DIF ¼ 1 þ C ln _ ,
(5)
where C is a constant, _ ¼ _=_0 , where _ is the actual equivalent strain rate and _ 0 ¼ 1:0 s1 is a reference strain rate used to normalize the strain rate. The equivalent strain rate is expressed as _ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x _y Þ2 þ ð_x _z Þ2 þ ð_y _z Þ2 þ 32ð_g2xy þ g_ 2xz þ g_ 2yz Þ, 9½ð_ ð6Þ where _x , _y , and _z are the three normal strain rates, and g_ xy , g_ xz , and g_ yz are the three shear strain rates. A polynomial equation of state is employed in the J–H model, the pressure variable is expressed as P ¼ k1 m þ k2 m2 þ k3 m3 þ Dp,
(7)
where k1, k2, k3 are constants (k1 is the bulk modulus), m is compression variable, m ¼ r/r01, r is the current density, r0 is the initial density, and Dp is an additional pressure increment determined from energy considerations to include the dilatation effect after compression failure. For tension with mo0, Eq. (7) is rewritten as P ¼ k1 m.
Table 1 Properties of granite rock mass with good quality Properties
Value
Density, r (kg/m3) Unconfined compressive strength of rock, sci (MPa) Critical tensile strength of intact rock, scri (MPa) Intact rock Young’s modulus, Ei (GPa) Intact rock shear modulus, Gi (GPa) Rock mass deformation modulus, Em (GPa) Poisson’s ratio, v Hoek–Brown constant for intact rock, mi s a m
2650 157 14 70 29 47 0.2 33 0.062 0.5 13.5
criterion [23,24] are listed in Table 1 [23]. The material constants A and N in the J–H model for the intact rock mass are derived from the result of the Hoek and Brown [24]. The Hoek–Brown criterion is an empirical criterion originally developed for applications in underground excavation design, which is expressed as s 1 ¼ s 3 þ sc
(8)
The damage variable is expressed as X Dp D¼ , fp
(9)
where Dep is the effective plastic strain during a cycle of integration and fp ¼ f ðPÞ is the plastic strain to fracture under a constant pressure P and expressed as D2
fp ¼ D1 ðP þ T Þ ,
(10)
where D1 and D2 are two damage constants. The material cannot undertake any plastic strain at P ¼ T, but fp increases with the increase of P. Details of the J–H material model can be found in [15,16,20]. The above J–H model was originally proposed for highpurity ceramics in a shock wave condition. In the present study, it is extended to simulate the fracture process of rock mass subjected to blasting-induced stress wave. Instead of the Hugoniot elastic limit for a shock wave condition, the uniaxial compressive strength of a rock material is used to normalize the stress and pressure variables in the strength models. The uniaxial compressive strength can be easily obtained and in the present loading condition, the loading intensity to the rock mass is not as high as in the shock wave condition produced by a hypervelocity impact. 3. Determination of J–H material constants For application of the J–H material model, a group of material constants should be given. In the present study, rock mass of granite is considered. Typical material constants for granite rock mass based on the Hoek–Brown
a ms3 þs , sc
(11)
where m, s, and a are material parameters, sc is the uniaxial compressive strength of the intact rock material, s1 and s3 are the axial and confining principal stresses, respectively. When s3 ¼ 0, the compressive strength of good quality rock mass (a ¼ 0.5) can be written as scm ¼
pffiffi s sc .
(12)
When s1 ¼ 0, the tensile strength of good quality rock mass is stm ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sc ðm þ m2 þ 4sÞ. 2
(13)
A series of compatible (s1,s3) can be obtained from Eq. (11). Correspondingly, a series of (P,s) can be derived, where P is the hydrostatic pressure and s is the effective stress defined in Eq. (4). The constants A and N are then obtained by curve-fitting using Eq. (1). It should be mentioned that the uniaxial compressive strength scm is used to normalize the stress and pressure variables and the dynamic increase factor (DIF) will be considered separately. The material constants A and N in the strength model are regressed to be 3.32 and 0.82, respectively. Due to lack of the experimental data for the fully fractured rock mass, it is assumed that the residual strength of fractured material is 30% of that of the intact rock mass, thus, B ¼ 1.0, M ¼ 0.6 are derived. For the damage parameters, D1 ¼ 0.04, D2 ¼ 1.0 are assumed with reference to the values for ceramic, which is also a brittle material, similar to rock.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
4. Dynamic loading produced by borehole blasting
5. Numerical simulations
To simulate borehole blasting, a radial pressure wave may be applied to the elements on the borehole surface. The amplitude and the duration of the pressure wave can be determined by the internal energy of the explosive and the size and geometry of the borehole. To represent a large range of the borehole pressure, the following general form of a pulse function is used [25],
5.1. Stress loading rate effect
P ¼ P0 ½eat ebt ,
(14)
where P is the pressure at time t, P0 is the peak pressure, and a and b are constants. For convenient representation of the rising and decaying phases, two constants, i.e. x ¼ 1=ðeat0 ebt0 Þ, t0 ¼ (1/(b/a))ln(b/a), are defined. Eq. (14) is then rewritten as P ¼ P0 x½eat ebt
(15)
which was employed by Cho and Kaneko [25] to study the influence of the rising and decaying phases on the dynamic fracture process of rock mass.
Fig. 2. Pressure–time curves for applied waveforms with b/a ¼ 1.5.
969
Mechanical properties of rock materials are affected by strain rate, which may influence the rock fracture pattern. Earlier researches [2,26] showed that higher strain rate loading tends to produce a larger crushed zone followed by larger numbers of short fractures while lower strain rate results in a relatively smaller crushed zone and fewer and longer radial fractures. To verify these observations of the strain rate effect, different waveforms of the blast loading are numerically simulated. Fig. 2 shows the pressure–time histories used for the numerical simulation. The rising time t0 is varied between 5, 10, and 100 ms with b/a ¼ 1.5. The peak pressure is kept at constant of 100 MPa. The corresponding loading rates for the three cases are calculated to be 20, 10, and 1.0 MPa/ms, respectively. The model consists of a borehole in rock mass with non-reflecting boundaries. Instantaneous detonation is assumed for the explosive column along its entire length. The numerical model is a slice of rock mass perpendicular to the borehole. Stress variation along the borehole direction is not considered due to the constraint of the computer capacity. Threedimensional solid elements are used to simulate the rock mass. A square rock mass with 6 m of the sides and a borehole with 0.05 m radius are adopted in the simulations. The simulated fracture patterns with different loading rates are compared in Fig. 3. When the loading rate of the borehole internal pressure is high (20 MPa/ms), only a crushed zone is created (see Fig. 3a). The blast energy is primarily dissipated in creating the crushed zone. When the loading rate decreases to a lower level (10 MPa/ms), there is a crushed zone followed by many short radial fractures (see Fig. 3b). When the loading rate decreases further to 1.0 MPa/ms, the crushed zone disappears and there are only comparatively longer radial fractures as shown in Fig. 3c. The results agree well with the findings reported in [2,26]. From these results, it can be seen that the J–H model can simulate both compressive crushed zone and tensile radial fractures. In the immediate vicinity of the explosion cavity,
Fig. 3. Fracture patterns for different stress loading rates, P0 ¼ 100 MPa, b/a ¼ 1.5: (a) 20 MPa/ms, (b) 10 MPa/ms, and (c) 1.0 MPa/ms.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
970
when explosion peak pressure exceeds the dynamic compressive strength of the media, a crushed zone forms. In the crushed zone, material yielding occurs which is associated with a very high rate of energy dissipation. Considering the efficiency and damage control in the smooth blasting, the crushing zone around the borehole should be avoided. Following the formation of the crushed zone, tensile radial fractures may also occur. In order to improve the efficiency of a blasting operation, which is indicated by long radial fractures with a minimal crushed zone, the loading rate should be controlled as low as possible. It has been well understood in the mining industry, which is fulfilled by the decoupled sources [2].
3m
2.5m
Rock mass
0.1m Borehole 1m
5.2. Effect of a nearby free face In practice, rock blasting is often operated near free boundaries to split the rock mass. When the compressive stress wave strikes a free face, it will be reflected backward to produce a tensile stress wave. If the reflected tensile wave is sufficiently strong (higher than the dynamic tensile strength of the media), ‘spalling’ occurs progressively from the free face back towards the borehole [27]. Fig. 4 shows a numerical model with three non-reflecting faces and one free face. The parameters for the pressure wave at the borehole surface are P0 ¼ 30 MPa, b/a ¼ 1.5, t0 ¼ 100 ms. The fracture propagation process is shown in Fig. 5. Radial fractures are first initiated and propagated from the borehole surface without directional preference (Fig. 5a). With further propagation of the stress wave, some radial fractures strike the free face and the stress wave reflects back. The tensile stress wave reflected from the free surface induces ‘spalling’ failure of the rock mass (Fig. 5b and c). It should be noted that some fractures (marked with white arrow in Fig. 5d) are extended to a considerable length during the ‘spalling’, which conforms the observation of other researchers [1,27]. 5.3. Effect of pre-compressive stress
Non-reflecting boundary Fig. 4. Analytical model with one free face.
Free face
In deep underground operations, blasting can be affected by earth stresses [27]. High earth stress may induce
Fig. 5. Crack propagation process in one-free-face model.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
fractures around the borehole. These fractures may be extended or suppressed by the blasting-induced stress field. On the other hand, the high earth stress may induce stress
971
concentrations around the wall of the cavity and further lead to non-uniformity in the fracture pattern. Stress fields are usually biaxial. It is only the difference between two principal stresses in the plane normal to the borehole axis that will cause an asymmetrical fracture pattern. Therefore, only a uniaxial lateral stress is considered in this study. In the present numerical simulation, two opposite edges of the numerical model are pre-loaded with four different lateral stresses of 2, 10, 30, and 50 MPa, respectively, while
Fig. 6. Fracture pattern for the model subjected to different pre-existing compressive stresses.
Borehole
Rock mass
Borehole
Joint Plane
Rock mass
Joint Plane Non-reflecting boundary Fig. 7. Analytical models with a joint plane: (a) joint plane parallel to the free face and (b) joint plane normal to the free face.
Fig. 8. Fracture propagation processes for the models with one joint plane: (a) joint plane parallel to the free face and (b) joint plane normal to the free face.
ARTICLE IN PRESS 972
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
the other two edges are non-reflecting boundaries and free of stress. When the rock mass system becomes stable under the lateral stress, a pressure wave (P0 ¼ 30 MPa, b/a ¼ 1.5, t0 ¼ 100 ms) is then applied on the borehole wall. The simulated fracture patterns of the model with different lateral compressive stresses are shown in Fig. 6, which clearly demonstrates that the fractures are mainly aligned in the direction of the pre-loading axis. It is because the tensile stress perpendicular to the pre-loading axis is suppressed. This result conforms to the findings by other researchers [1,2]. 5.4. Effects of joint planes During the stress wave traveling, joint planes are frequently encountered in the rock mass. And these discontinuities may have strong effect on the fracture patterns. Zhu et al. [28] have investigated the effect of joint spacing, joint width, and joint filling material on the fracture patterns. In the current study, the joint orientation effect will be studied. Two cases are presented: one with the joint plane parallel to the free face (see Fig. 7a); and the other one with the joint plane normal to the free face (see Fig. 7b). The pressure wave and the borehole diameter remain the same. The joint is simplified as a gap with a width of 2 mm. A surface to surface contact algorithm which is available in LS-DYNA is applied to simulate the interaction of rock masses at the two sides of the joint. The simulated fractures for the above two cases are compared in Fig. 8. At the beginning, the radial fractures starts randomly from the borehole surface and extends toward the edges. When the stress wave strikes the joint, a tensile wave is reflected from the joint plane because the joint plane has an initial gap that acts as a free face. It is observed that the reflected stress wave significantly affects the fracture pattern. When the reflected tensile stress is higher than the tensile strength of the rock, ‘spalling’ occurs progressively from the joint plane towards the borehole (see Fig. 8). The reflected stress wave also interacts with some of the radial fractures and creates a
spalling zone near the joint plane. When the joint plane is parallel to the free face (Fig. 8a), the block between the joint and the free face are fragmented, which indicates that the stress wave is transferred to the block and causes tensile fractures at the free face. On the other hand, when the joint is normal to the free face (Fig. 8b), there is no damage incurred in the block ahead of the joint. The reason is that the stress wave is transferred through the non-reflecting boundary. This finding justifies the existing pre-split technique, which is designed to enhance free face quality and to prevent overbreak of blast operations in rock mass. However, it is worth noting that different angles between the joint plane and the free face may result in different form of damage. 5.5. Fracture control techniques Control of fracture initiation and propagation has wide applications in rock blasting operation. For example, in underground excavation, smooth blasting technique is applied to achieve the required smoothness and the specific cavity extents. There are two major blast control methods to fulfill the smooth blasting technique, namely redial notches and pre-split charge holder. The present study also investigates the efficiencies of these two techniques numerically. (1) Radial notches Radial notches are drilled at the borehole surface, where the fracture is intended to start. The fracture propagation is controlled by the stress concentrations at the tips of the radial notches along the borehole, which lead to the extension of the notch. In this section, two radial notches at the borehole surface are simulated. After detonation of the charge, the explosive gases will permeate into the pre-existing notches and exert a gas pressure on the notch walls. The radius of the borehole and the length of notches are 0.025 and 0.008 m, respectively. A 2 m square rock mass is
Fig. 9. Comparison of fracture patterns for the models with and without notches: (a) without notches, (b) with two diametrical notches, and (c) with two perpendicular notches.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
considered in the simulation. The simulation result is shown in Fig. 9. It is seen that the fracture starts and propagates mainly along the notches, while the fracture extents in other directions are much smaller. The effect of the notch wall pressure is also evaluated by introducing the ratio of notch wall pressure over the borehole wall pressure l. A larger pressure ratio indicates a better stemmed system in the blasting operation. In this study, the pressure ratio l is given various values of 0, 0.05, 0.1, 0.2, 0.3, and 0.4, respectively. When l is 0.4, the final length of the fracture from the notch is about 1 m. For convenient presentation of the simulation results, a parameter, length increase factor (LIF) is defined as the ratio of the final fracture length with the notch wall pressure to the final fracture length without notch wall pressure. The simulated variation of LIF with the pressure ratio l is shown in Fig. 10. It can be observed that the notch wall pressure has a positive effect on fracture propagation along the radial notches, while it reduces the number of fractures in other directions. A stemmed borehole favors the fracture control ability of the radial notch method. This result agrees well with Mohanty’s report [29], in which it was shown that the stress intensity ahead of the notches increases with the pressure ratio, l.
Fig. 10. Relationship between of length increase factor (LIF) and pressure ratio (l).
973
(2) Pre-split charge holder The pre-split charge holder is a thin-walled steel tube which contains longitudinal slits to control the location of fracture initiation and propagation direction. Fourney et al. [30] carried out experiments on the split charge-holding device and testified its fracture control ability. In our study, we simulate a model with a charge holder to evaluate its effectiveness on fracture control. A steel charge holder is considered and an elastic, isotropic plastic hardening material model is used to simulate the charge holder deformation. The J–H material model is again adopted for the rock mass. The simulated fracture patterns of the model with two-slit charge holder is shown in Fig. 11. The results demonstrate that the radial fractures initiate and propagate along the direction of the slits because of stress concentration. The two fractures in the direction of the diametrical slits form a single fracture plane. It is also interesting to note that there are no fractures in other directions, which is probably due to the reinforcement of the steel charge holder. Based on the numerical results, the two-slit charge holder can be used in presplitting operation effectively when a single fracture plane is desired. Charge holders with three slits are also simulated and the result is shown in Fig. 12. The directions of the three control planes are marked by white dashed lines. It can be found that fractures initiate and propagate along the three control planes with very little deviation, which again indicate the effectiveness of fracture control ability of the charge holder. In simple summary of the simulation results, it can be said that the charge holder is effective in controlling the initiation and propagation of fractures. By employing appropriate number of slits in a charge holder, desired fracture pattern can be fulfilled.
Fig. 11. (a)–(c) Two-slit charge holder with different orientations and the corresponding fracture patterns.
ARTICLE IN PRESS 974
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
(6) Charge holder is effective in controlling the initiation and propagation of fractures. Two-slit charge holder can be used in pre-splitting operation, while three-slit charge holder can be applied in smooth blasting where more fragmentation is desired at one side of the fractured plane.
References
Fig. 12. (a) and (b) Three-slit charge holder of different orientations and the corresponding fracture patterns.
6. Conclusions In the present study, the J–H material model is implemented into LS-DYNA to simulate the blastinginduced rock fractures. The influences of key parameters on the rock fracture pattern are extensively investigated. From the numerical simulation results, we come to the following conclusions: (1) J–H model can be applied to simulate effectively both compressive crushed zone and tensile radial fractures. (2) The stress-loading rate has a significant effect on the fracture pattern. When the loading rate is very high, only a crushed zone is created. When the loading rate is relatively low, a crushed zone is formed followed by many radial tensile fractures. And if the loading rate is very low, there will only be radial fractures. (3) For the one-free-face model, radial fractures are firstly initiated and propagated without directional preference from the borehole surface. With further propagation of the stress wave, some radial fractures strike the free face and the stress wave reflects back. The tensile stress wave induces ‘spalling’ failure of rock mass. Besides ‘spalling’, some fractures are extended to a great length. (4) Fractures are aligned in the direction of the precompressive stress axis. The larger the pre-loading, the more obvious is this phenomenon. (5) When the pre-existing joint is parallel to the free face, the block behind the joint will be well fragmented due to the free face. When the pre-existing joint is normal to the free face, there is no damage in the block behind the joint. It justifies the pre-splitting technique, which is designed to prevent overbreak.
[1] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8:181–202. [2] Donze FV, Bouchez J, Magnier SA. Modeling fractures in rock blasting. Int J Rock Mech Min Sci 1997;34(8):1153–63. [3] Preece DS, Burchell SL. Variation of spherical element packing angle and its influence on computer simulations of blasting induced rock motion. In: Proceedings of the 1st international conference on discrete element methods, MIT, Cambridge, MA; 1993. p. 25–99. [4] Donze FV, Magnier SA. Formulation of a three-dimensional numerical model of brittle behavior. Geophys J Int 1995;122:790–802. [5] Potyondy DO, Cundall PA. A bounded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329–64. [6] Donze FV, Magnier SA, Bouchez J. Numerical modeling of a highly explosive source in an elastic–brittle rock mass. J Geophys Res 1996;101:3103–12. [7] Song J, Kim K. Micromechanical modeling of the dynamic fracture process during rock blasting. Int J Rock Mech Min Sci 1996;33(4): 387–94. [8] Grady DE, Kipp ME. Continuum modelling of explosive fracture in oil shale. Int J Rock Mech Min Sci 1980;17:147–57. [9] Taylor LM, Chen EP, Kuszmaul JS. Micro-crack induced damage accumulation in brittle rock under dynamic loading. Comp Meth Appl Mech Eng 1986;55:301–20. [10] Grady DE, Kipp ME. Dynamic rock fragmentation. In: Atkinson BK, editor. Fracture mechanics of rock. London: Academic Press; 1987. p. 429–75. [11] Yang R, Bawden WF, Katsabanis PD. A new constitutive model for blast damage. Int J Rock Mech Min Sci 1996;33(3):245–54. [12] Liu L, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci 1997;34:217–31. [13] Hao H, Ma GW, Zhou YX. Numerical simulation of underground explosions. Fragblast—Int J Blast Fragment 1998;2:383–95. [14] Ma GW, Hao H, Zhou YX. Modeling of wave propagation induced by underground explosion. Comp Geotech 1998;22(3/4):283–303. [15] Holmquist TJ, Templeton DW, Bishnoi KD. Constitutive modeling of aluminum nitride for large strain, high-strain rate, and highpressure applications. Int J Impact Eng 2001;25:211–31. [16] Johnson GR, Holmquist TJ, Beissel SR. Response of aluminum nitride (including a phase change) to large strains, high strain rates, and high pressures. J Appl Phys 2003;94(3):1639–46. [17] Zhang YQ, Lu Y, Hao H. Analysis of fragment size and ejection velocity at high strain rate. Int J Mech Sci 2004;46:27–34. [18] Lu Y, Xu K. Modeling of dynamic behavior of concrete materials under blast loading. Int J Solids Struct 2004;41:131–43. [19] Munjiza A, Owen RJ, Bicanic N. A combined finite–discrete element method in transient dynamics of fracturing solids. Eng Comput 1995; 12:145–74. [20] Johnson GR, Holmquist TJ. An improved computational constitutive model for brittle materials. In: Schmidt SC, Shaner JW, Samara GA, Ross M, editors. High pressure science and technology—1993. Woodbury, NY: AIP Press; 1994. p. 981–4. [21] LSTC. LS-DYNA theoretical manual, Livermore Software Technology Corporation, Livermore, CA, 1998. [22] LSTC. LS-DYNA keyword user’s manual, Version 970, Livermore Software Technology Corporation, Livermore, CA, 2003.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975 [23] Hoek E, Brown ET. Underground excavations in rock. London: Institute of Mining and Metallurgy; 1980. [24] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci 1997;34(8):1165–86. [25] Cho SH, Kaneko K. Influence of the applied pressure waveform on the dynamic fracture processes in rock. Int J Rock Mech Min Sci 2004;41:771–84. [26] Cho SH, Ogata Y, Kaneko K. Strain rate dependency of the dynamic tensile strength of rock. Int J Rock Mech Min Sci 2003;40:763–77.
975
[27] Hagan TN. Rock breakage by explosives. Acta Astron 1979;6: 329–40. [28] Zhu ZM, Mohanty B, Xie HP. Numerical investigation of blastinginduced crack initiation and propagation in rocks. Int J Rock Mech Min Sci 2007;44:412–24. [29] Mohanty B. Explosion generated fractures in rock and rock-like materials. Eng Fract Mech 1990;35(4/5):889–98. [30] Fourney WL, Dally JW, Holloway DC. Controlled blasting with ligamented charge holders. Int J Rock Mech Min Sci 1978;15:121–9.
International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975 www.elsevier.com/locate/ijrmms
Numerical simulation of blasting-induced rock fractures G.W. Ma, X.M. An School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore Received 3 January 2007; received in revised form 3 December 2007; accepted 4 December 2007 Available online 20 February 2008
Abstract In the present study, the Johnson–Holmquist (J–H) material model is implemented into the commercial software LS-DYNA through user-subroutines to simulate the blasting-induced rock fractures. The J–H model consists of strength models for both intact and fully fractured materials, a polynomial equation of state, and a damage model that represents the material from an intact state to a fully fractured state. Influences of the key parameters in smooth blasting, viz., loading rate, distance from a free face, earth stress, and preexisting joint planes, etc., on fracture patterns are explored. According to the simulation results, the rock fracture pattern is significantly influenced by the loading rate. Fracture control techniques, namely, notched borehole and charge holder with slits are also simulated. Effectiveness of the fracture control techniques is demonstrated. The numerical simulation in the present study reproduces some of the well-known phenomena observed by other researchers. It has the potential to be applied in practical blast control and gas and hydraulic fracturing engineering. r 2007 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; LS-DYNA; Johnson–Holmquist material model; Rock fracture; Blast control
1. Introduction Drill and blast method has been widely used in underground excavation and construction for many years and it is still a popular method of rock breakage. However, a clear description of the fracture process investigated in the rock is still lacking. When an explosive is detonated, a chemical reaction occurs very rapidly and a relatively small quantity of explosive is converted into gas of very high temperature. This reaction results in two types of loadings applied on the borehole wall, namely a stress wave and a gas pressure with longer duration [1,2]. Rock fractures will be subsequently initiated and propagated in the surrounding rock mass. In a coupled borehole blasting, a crushing zone may be generated at a specific depth from the borehole surface. Recent studies have revealed that the stress wave is responsible for initiation of the crushing zone and the surrounding radial fractures, while the gas pressure further extends the fractures [2].
Corresponding author. Tel.: +65 6790 4984; fax: +65 6790 6841.
E-mail address: [email protected] (G.W. Ma). 1365-1609/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2007.12.002
Investigation of rock blasting and its effects by scaled or full size experiment is very expensive and time-consuming. On the other hand, numerical method, derived from sound mechanical principles and validated against experimental data, indicates a promising approach to reveal the blastinginduced rock fracture process. Considerable efforts have been directed towards developing discontinuum [2,3–7], continuum [8–18], and continuum–discontinuum [19] constitutive models of the rock materials. For a discontinuum model, the rock mass is represented by an assembly of jointed elements (e.g. particles and springs) and a fracture is initiated when the stress in the joint between two elements exceeds a critical value. However, it requires large amount of data about the spatial distribution, geometrical and mechanical properties of the rock joints, which are normally not readily available. Continuum damage mechanics is also considered appropriate to describe the blast damage process from a stress wave. The continuum damage mechanics is based on the fact that the brittle rock has pre-existing defects. The material failure is a process of growth and nucleation of these defects. Material strength depends on the strength and distribution of the critical defects. A single variable, or
ARTICLE IN PRESS 967
Intact (D=0)
Damaged (0
HEL
Fully Fractured (D=1.0)
T*
Normalized pressure, P*
D=1
Pressure, P
a set of variables, can be used to describe the changing state of a material. For a continuum damage model, it is critical to define an appropriate damage evolutionary law. Johnson and Holmquist [15,16,20] suggested a constitutive mode for brittle materials (J–H model) that has been applied to high-purity ceramics. The J–H model consists of strength models for both intact and fully fractured materials, a polynomial equation of state, and a damage model that represents the material from an intact state to a fully fractured state. It is pressure- and strain-ratedependent and suitable for large strain, high strain rate and high pressure problems. In the present study, the J–H model is implemented into a commercial software LS-DYNA [21,22], which is an explicit, three-dimensional, dynamic, nonlinear finite element program, through its user-subroutines. Blastinginduced rock fractures are simulated for different loading and boundary conditions of a borehole. It is found that the numerical simulation in the present study reproduces some of the well-known phenomena observed by other researchers. Fracture control techniques, namely, notched borehole and charge holder with slits are also simulated. Effectiveness of the fracture control techniques is clearly demonstrated.
Normalized Equivalent Stress, σ∗
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
0
The J–H material model is schematically illustrated in Fig. 1. The model gives the intact and the totally fractured strengths (Fig. 1a), a polynomial equation of state (Fig. 1b), and a damage model that represents the material from an intact state to a fractured state (Fig. 1c). The intact strength (upper solid line in Fig. 1a) is given as
Relative Volume, μ
si ¼ AðP þ T ÞN ð1 þ C ln _ Þ.
Equivalent Plastic Fracture Strain, εpf
2. J–H material model
(1)
The fully fractured strength (lower solid line in Fig. 1a) is given as sf ¼ BðP ÞM ð1 þ C ln _ Þ.
(2)
The strength of damaged material (dashed line in Fig. 1a) is expressed as s ¼ si Dðsi sf Þ,
(3)
where A, N, B, M, and C are material constants; si , sf , s, P, and T are effective stress of the intact material, effective stress of the fully fractured material, current effective stress, pressure, and tensile strength, respectively, which are normalized by the following expressions: si ¼ si =sHEL , sf ¼ sf =sHEL , s ¼ s=sHEL , P ¼ P=PHEL , and T ¼ T=PHEL , in which P is the current hydrostatic pressure, T is the maximum tensile hydrostatic pressure the material can withstand, sHEL is the effective stress at the Hugoniot elastic limit (HEL) which is defined in the original J–H model to characterize material yielding in a shock wave condition, PHEL is the associated hydrostatic pressure at the HEL, D is the damage variable, s is the
T∗
Normalized pressure P∗
Fig. 1. Schematic illustration of J–H material model: (a) strength model, (b) pressure–volume relationship, and (c) damage model.
effective stress which has the general form of s¼ (
ðsx sy Þ2 þ ðsx sz Þ2 þ ðsy sz Þ2 þ 6ðt2xy þ t2xz þ t2yz Þ 2
)1=2 ,
ð4Þ where sx, sy, and sz are the three normal stresses, and txy, tyz, and tzx are the three shear stresses. It should be mentioned that rate dependence of the J–H model is represented by a dynamic strength increase factor
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
968
(DIF) as follows,
DIF ¼ 1 þ C ln _ ,
(5)
where C is a constant, _ ¼ _=_0 , where _ is the actual equivalent strain rate and _ 0 ¼ 1:0 s1 is a reference strain rate used to normalize the strain rate. The equivalent strain rate is expressed as _ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x _y Þ2 þ ð_x _z Þ2 þ ð_y _z Þ2 þ 32ð_g2xy þ g_ 2xz þ g_ 2yz Þ, 9½ð_ ð6Þ where _x , _y , and _z are the three normal strain rates, and g_ xy , g_ xz , and g_ yz are the three shear strain rates. A polynomial equation of state is employed in the J–H model, the pressure variable is expressed as P ¼ k1 m þ k2 m2 þ k3 m3 þ Dp,
(7)
where k1, k2, k3 are constants (k1 is the bulk modulus), m is compression variable, m ¼ r/r01, r is the current density, r0 is the initial density, and Dp is an additional pressure increment determined from energy considerations to include the dilatation effect after compression failure. For tension with mo0, Eq. (7) is rewritten as P ¼ k1 m.
Table 1 Properties of granite rock mass with good quality Properties
Value
Density, r (kg/m3) Unconfined compressive strength of rock, sci (MPa) Critical tensile strength of intact rock, scri (MPa) Intact rock Young’s modulus, Ei (GPa) Intact rock shear modulus, Gi (GPa) Rock mass deformation modulus, Em (GPa) Poisson’s ratio, v Hoek–Brown constant for intact rock, mi s a m
2650 157 14 70 29 47 0.2 33 0.062 0.5 13.5
criterion [23,24] are listed in Table 1 [23]. The material constants A and N in the J–H model for the intact rock mass are derived from the result of the Hoek and Brown [24]. The Hoek–Brown criterion is an empirical criterion originally developed for applications in underground excavation design, which is expressed as s 1 ¼ s 3 þ sc
(8)
The damage variable is expressed as X Dp D¼ , fp
(9)
where Dep is the effective plastic strain during a cycle of integration and fp ¼ f ðPÞ is the plastic strain to fracture under a constant pressure P and expressed as D2
fp ¼ D1 ðP þ T Þ ,
(10)
where D1 and D2 are two damage constants. The material cannot undertake any plastic strain at P ¼ T, but fp increases with the increase of P. Details of the J–H material model can be found in [15,16,20]. The above J–H model was originally proposed for highpurity ceramics in a shock wave condition. In the present study, it is extended to simulate the fracture process of rock mass subjected to blasting-induced stress wave. Instead of the Hugoniot elastic limit for a shock wave condition, the uniaxial compressive strength of a rock material is used to normalize the stress and pressure variables in the strength models. The uniaxial compressive strength can be easily obtained and in the present loading condition, the loading intensity to the rock mass is not as high as in the shock wave condition produced by a hypervelocity impact. 3. Determination of J–H material constants For application of the J–H material model, a group of material constants should be given. In the present study, rock mass of granite is considered. Typical material constants for granite rock mass based on the Hoek–Brown
a ms3 þs , sc
(11)
where m, s, and a are material parameters, sc is the uniaxial compressive strength of the intact rock material, s1 and s3 are the axial and confining principal stresses, respectively. When s3 ¼ 0, the compressive strength of good quality rock mass (a ¼ 0.5) can be written as scm ¼
pffiffi s sc .
(12)
When s1 ¼ 0, the tensile strength of good quality rock mass is stm ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sc ðm þ m2 þ 4sÞ. 2
(13)
A series of compatible (s1,s3) can be obtained from Eq. (11). Correspondingly, a series of (P,s) can be derived, where P is the hydrostatic pressure and s is the effective stress defined in Eq. (4). The constants A and N are then obtained by curve-fitting using Eq. (1). It should be mentioned that the uniaxial compressive strength scm is used to normalize the stress and pressure variables and the dynamic increase factor (DIF) will be considered separately. The material constants A and N in the strength model are regressed to be 3.32 and 0.82, respectively. Due to lack of the experimental data for the fully fractured rock mass, it is assumed that the residual strength of fractured material is 30% of that of the intact rock mass, thus, B ¼ 1.0, M ¼ 0.6 are derived. For the damage parameters, D1 ¼ 0.04, D2 ¼ 1.0 are assumed with reference to the values for ceramic, which is also a brittle material, similar to rock.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
4. Dynamic loading produced by borehole blasting
5. Numerical simulations
To simulate borehole blasting, a radial pressure wave may be applied to the elements on the borehole surface. The amplitude and the duration of the pressure wave can be determined by the internal energy of the explosive and the size and geometry of the borehole. To represent a large range of the borehole pressure, the following general form of a pulse function is used [25],
5.1. Stress loading rate effect
P ¼ P0 ½eat ebt ,
(14)
where P is the pressure at time t, P0 is the peak pressure, and a and b are constants. For convenient representation of the rising and decaying phases, two constants, i.e. x ¼ 1=ðeat0 ebt0 Þ, t0 ¼ (1/(b/a))ln(b/a), are defined. Eq. (14) is then rewritten as P ¼ P0 x½eat ebt
(15)
which was employed by Cho and Kaneko [25] to study the influence of the rising and decaying phases on the dynamic fracture process of rock mass.
Fig. 2. Pressure–time curves for applied waveforms with b/a ¼ 1.5.
969
Mechanical properties of rock materials are affected by strain rate, which may influence the rock fracture pattern. Earlier researches [2,26] showed that higher strain rate loading tends to produce a larger crushed zone followed by larger numbers of short fractures while lower strain rate results in a relatively smaller crushed zone and fewer and longer radial fractures. To verify these observations of the strain rate effect, different waveforms of the blast loading are numerically simulated. Fig. 2 shows the pressure–time histories used for the numerical simulation. The rising time t0 is varied between 5, 10, and 100 ms with b/a ¼ 1.5. The peak pressure is kept at constant of 100 MPa. The corresponding loading rates for the three cases are calculated to be 20, 10, and 1.0 MPa/ms, respectively. The model consists of a borehole in rock mass with non-reflecting boundaries. Instantaneous detonation is assumed for the explosive column along its entire length. The numerical model is a slice of rock mass perpendicular to the borehole. Stress variation along the borehole direction is not considered due to the constraint of the computer capacity. Threedimensional solid elements are used to simulate the rock mass. A square rock mass with 6 m of the sides and a borehole with 0.05 m radius are adopted in the simulations. The simulated fracture patterns with different loading rates are compared in Fig. 3. When the loading rate of the borehole internal pressure is high (20 MPa/ms), only a crushed zone is created (see Fig. 3a). The blast energy is primarily dissipated in creating the crushed zone. When the loading rate decreases to a lower level (10 MPa/ms), there is a crushed zone followed by many short radial fractures (see Fig. 3b). When the loading rate decreases further to 1.0 MPa/ms, the crushed zone disappears and there are only comparatively longer radial fractures as shown in Fig. 3c. The results agree well with the findings reported in [2,26]. From these results, it can be seen that the J–H model can simulate both compressive crushed zone and tensile radial fractures. In the immediate vicinity of the explosion cavity,
Fig. 3. Fracture patterns for different stress loading rates, P0 ¼ 100 MPa, b/a ¼ 1.5: (a) 20 MPa/ms, (b) 10 MPa/ms, and (c) 1.0 MPa/ms.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
970
when explosion peak pressure exceeds the dynamic compressive strength of the media, a crushed zone forms. In the crushed zone, material yielding occurs which is associated with a very high rate of energy dissipation. Considering the efficiency and damage control in the smooth blasting, the crushing zone around the borehole should be avoided. Following the formation of the crushed zone, tensile radial fractures may also occur. In order to improve the efficiency of a blasting operation, which is indicated by long radial fractures with a minimal crushed zone, the loading rate should be controlled as low as possible. It has been well understood in the mining industry, which is fulfilled by the decoupled sources [2].
3m
2.5m
Rock mass
0.1m Borehole 1m
5.2. Effect of a nearby free face In practice, rock blasting is often operated near free boundaries to split the rock mass. When the compressive stress wave strikes a free face, it will be reflected backward to produce a tensile stress wave. If the reflected tensile wave is sufficiently strong (higher than the dynamic tensile strength of the media), ‘spalling’ occurs progressively from the free face back towards the borehole [27]. Fig. 4 shows a numerical model with three non-reflecting faces and one free face. The parameters for the pressure wave at the borehole surface are P0 ¼ 30 MPa, b/a ¼ 1.5, t0 ¼ 100 ms. The fracture propagation process is shown in Fig. 5. Radial fractures are first initiated and propagated from the borehole surface without directional preference (Fig. 5a). With further propagation of the stress wave, some radial fractures strike the free face and the stress wave reflects back. The tensile stress wave reflected from the free surface induces ‘spalling’ failure of the rock mass (Fig. 5b and c). It should be noted that some fractures (marked with white arrow in Fig. 5d) are extended to a considerable length during the ‘spalling’, which conforms the observation of other researchers [1,27]. 5.3. Effect of pre-compressive stress
Non-reflecting boundary Fig. 4. Analytical model with one free face.
Free face
In deep underground operations, blasting can be affected by earth stresses [27]. High earth stress may induce
Fig. 5. Crack propagation process in one-free-face model.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
fractures around the borehole. These fractures may be extended or suppressed by the blasting-induced stress field. On the other hand, the high earth stress may induce stress
971
concentrations around the wall of the cavity and further lead to non-uniformity in the fracture pattern. Stress fields are usually biaxial. It is only the difference between two principal stresses in the plane normal to the borehole axis that will cause an asymmetrical fracture pattern. Therefore, only a uniaxial lateral stress is considered in this study. In the present numerical simulation, two opposite edges of the numerical model are pre-loaded with four different lateral stresses of 2, 10, 30, and 50 MPa, respectively, while
Fig. 6. Fracture pattern for the model subjected to different pre-existing compressive stresses.
Borehole
Rock mass
Borehole
Joint Plane
Rock mass
Joint Plane Non-reflecting boundary Fig. 7. Analytical models with a joint plane: (a) joint plane parallel to the free face and (b) joint plane normal to the free face.
Fig. 8. Fracture propagation processes for the models with one joint plane: (a) joint plane parallel to the free face and (b) joint plane normal to the free face.
ARTICLE IN PRESS 972
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
the other two edges are non-reflecting boundaries and free of stress. When the rock mass system becomes stable under the lateral stress, a pressure wave (P0 ¼ 30 MPa, b/a ¼ 1.5, t0 ¼ 100 ms) is then applied on the borehole wall. The simulated fracture patterns of the model with different lateral compressive stresses are shown in Fig. 6, which clearly demonstrates that the fractures are mainly aligned in the direction of the pre-loading axis. It is because the tensile stress perpendicular to the pre-loading axis is suppressed. This result conforms to the findings by other researchers [1,2]. 5.4. Effects of joint planes During the stress wave traveling, joint planes are frequently encountered in the rock mass. And these discontinuities may have strong effect on the fracture patterns. Zhu et al. [28] have investigated the effect of joint spacing, joint width, and joint filling material on the fracture patterns. In the current study, the joint orientation effect will be studied. Two cases are presented: one with the joint plane parallel to the free face (see Fig. 7a); and the other one with the joint plane normal to the free face (see Fig. 7b). The pressure wave and the borehole diameter remain the same. The joint is simplified as a gap with a width of 2 mm. A surface to surface contact algorithm which is available in LS-DYNA is applied to simulate the interaction of rock masses at the two sides of the joint. The simulated fractures for the above two cases are compared in Fig. 8. At the beginning, the radial fractures starts randomly from the borehole surface and extends toward the edges. When the stress wave strikes the joint, a tensile wave is reflected from the joint plane because the joint plane has an initial gap that acts as a free face. It is observed that the reflected stress wave significantly affects the fracture pattern. When the reflected tensile stress is higher than the tensile strength of the rock, ‘spalling’ occurs progressively from the joint plane towards the borehole (see Fig. 8). The reflected stress wave also interacts with some of the radial fractures and creates a
spalling zone near the joint plane. When the joint plane is parallel to the free face (Fig. 8a), the block between the joint and the free face are fragmented, which indicates that the stress wave is transferred to the block and causes tensile fractures at the free face. On the other hand, when the joint is normal to the free face (Fig. 8b), there is no damage incurred in the block ahead of the joint. The reason is that the stress wave is transferred through the non-reflecting boundary. This finding justifies the existing pre-split technique, which is designed to enhance free face quality and to prevent overbreak of blast operations in rock mass. However, it is worth noting that different angles between the joint plane and the free face may result in different form of damage. 5.5. Fracture control techniques Control of fracture initiation and propagation has wide applications in rock blasting operation. For example, in underground excavation, smooth blasting technique is applied to achieve the required smoothness and the specific cavity extents. There are two major blast control methods to fulfill the smooth blasting technique, namely redial notches and pre-split charge holder. The present study also investigates the efficiencies of these two techniques numerically. (1) Radial notches Radial notches are drilled at the borehole surface, where the fracture is intended to start. The fracture propagation is controlled by the stress concentrations at the tips of the radial notches along the borehole, which lead to the extension of the notch. In this section, two radial notches at the borehole surface are simulated. After detonation of the charge, the explosive gases will permeate into the pre-existing notches and exert a gas pressure on the notch walls. The radius of the borehole and the length of notches are 0.025 and 0.008 m, respectively. A 2 m square rock mass is
Fig. 9. Comparison of fracture patterns for the models with and without notches: (a) without notches, (b) with two diametrical notches, and (c) with two perpendicular notches.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
considered in the simulation. The simulation result is shown in Fig. 9. It is seen that the fracture starts and propagates mainly along the notches, while the fracture extents in other directions are much smaller. The effect of the notch wall pressure is also evaluated by introducing the ratio of notch wall pressure over the borehole wall pressure l. A larger pressure ratio indicates a better stemmed system in the blasting operation. In this study, the pressure ratio l is given various values of 0, 0.05, 0.1, 0.2, 0.3, and 0.4, respectively. When l is 0.4, the final length of the fracture from the notch is about 1 m. For convenient presentation of the simulation results, a parameter, length increase factor (LIF) is defined as the ratio of the final fracture length with the notch wall pressure to the final fracture length without notch wall pressure. The simulated variation of LIF with the pressure ratio l is shown in Fig. 10. It can be observed that the notch wall pressure has a positive effect on fracture propagation along the radial notches, while it reduces the number of fractures in other directions. A stemmed borehole favors the fracture control ability of the radial notch method. This result agrees well with Mohanty’s report [29], in which it was shown that the stress intensity ahead of the notches increases with the pressure ratio, l.
Fig. 10. Relationship between of length increase factor (LIF) and pressure ratio (l).
973
(2) Pre-split charge holder The pre-split charge holder is a thin-walled steel tube which contains longitudinal slits to control the location of fracture initiation and propagation direction. Fourney et al. [30] carried out experiments on the split charge-holding device and testified its fracture control ability. In our study, we simulate a model with a charge holder to evaluate its effectiveness on fracture control. A steel charge holder is considered and an elastic, isotropic plastic hardening material model is used to simulate the charge holder deformation. The J–H material model is again adopted for the rock mass. The simulated fracture patterns of the model with two-slit charge holder is shown in Fig. 11. The results demonstrate that the radial fractures initiate and propagate along the direction of the slits because of stress concentration. The two fractures in the direction of the diametrical slits form a single fracture plane. It is also interesting to note that there are no fractures in other directions, which is probably due to the reinforcement of the steel charge holder. Based on the numerical results, the two-slit charge holder can be used in presplitting operation effectively when a single fracture plane is desired. Charge holders with three slits are also simulated and the result is shown in Fig. 12. The directions of the three control planes are marked by white dashed lines. It can be found that fractures initiate and propagate along the three control planes with very little deviation, which again indicate the effectiveness of fracture control ability of the charge holder. In simple summary of the simulation results, it can be said that the charge holder is effective in controlling the initiation and propagation of fractures. By employing appropriate number of slits in a charge holder, desired fracture pattern can be fulfilled.
Fig. 11. (a)–(c) Two-slit charge holder with different orientations and the corresponding fracture patterns.
ARTICLE IN PRESS 974
G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975
(6) Charge holder is effective in controlling the initiation and propagation of fractures. Two-slit charge holder can be used in pre-splitting operation, while three-slit charge holder can be applied in smooth blasting where more fragmentation is desired at one side of the fractured plane.
References
Fig. 12. (a) and (b) Three-slit charge holder of different orientations and the corresponding fracture patterns.
6. Conclusions In the present study, the J–H material model is implemented into LS-DYNA to simulate the blastinginduced rock fractures. The influences of key parameters on the rock fracture pattern are extensively investigated. From the numerical simulation results, we come to the following conclusions: (1) J–H model can be applied to simulate effectively both compressive crushed zone and tensile radial fractures. (2) The stress-loading rate has a significant effect on the fracture pattern. When the loading rate is very high, only a crushed zone is created. When the loading rate is relatively low, a crushed zone is formed followed by many radial tensile fractures. And if the loading rate is very low, there will only be radial fractures. (3) For the one-free-face model, radial fractures are firstly initiated and propagated without directional preference from the borehole surface. With further propagation of the stress wave, some radial fractures strike the free face and the stress wave reflects back. The tensile stress wave induces ‘spalling’ failure of rock mass. Besides ‘spalling’, some fractures are extended to a great length. (4) Fractures are aligned in the direction of the precompressive stress axis. The larger the pre-loading, the more obvious is this phenomenon. (5) When the pre-existing joint is parallel to the free face, the block behind the joint will be well fragmented due to the free face. When the pre-existing joint is normal to the free face, there is no damage in the block behind the joint. It justifies the pre-splitting technique, which is designed to prevent overbreak.
[1] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8:181–202. [2] Donze FV, Bouchez J, Magnier SA. Modeling fractures in rock blasting. Int J Rock Mech Min Sci 1997;34(8):1153–63. [3] Preece DS, Burchell SL. Variation of spherical element packing angle and its influence on computer simulations of blasting induced rock motion. In: Proceedings of the 1st international conference on discrete element methods, MIT, Cambridge, MA; 1993. p. 25–99. [4] Donze FV, Magnier SA. Formulation of a three-dimensional numerical model of brittle behavior. Geophys J Int 1995;122:790–802. [5] Potyondy DO, Cundall PA. A bounded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329–64. [6] Donze FV, Magnier SA, Bouchez J. Numerical modeling of a highly explosive source in an elastic–brittle rock mass. J Geophys Res 1996;101:3103–12. [7] Song J, Kim K. Micromechanical modeling of the dynamic fracture process during rock blasting. Int J Rock Mech Min Sci 1996;33(4): 387–94. [8] Grady DE, Kipp ME. Continuum modelling of explosive fracture in oil shale. Int J Rock Mech Min Sci 1980;17:147–57. [9] Taylor LM, Chen EP, Kuszmaul JS. Micro-crack induced damage accumulation in brittle rock under dynamic loading. Comp Meth Appl Mech Eng 1986;55:301–20. [10] Grady DE, Kipp ME. Dynamic rock fragmentation. In: Atkinson BK, editor. Fracture mechanics of rock. London: Academic Press; 1987. p. 429–75. [11] Yang R, Bawden WF, Katsabanis PD. A new constitutive model for blast damage. Int J Rock Mech Min Sci 1996;33(3):245–54. [12] Liu L, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci 1997;34:217–31. [13] Hao H, Ma GW, Zhou YX. Numerical simulation of underground explosions. Fragblast—Int J Blast Fragment 1998;2:383–95. [14] Ma GW, Hao H, Zhou YX. Modeling of wave propagation induced by underground explosion. Comp Geotech 1998;22(3/4):283–303. [15] Holmquist TJ, Templeton DW, Bishnoi KD. Constitutive modeling of aluminum nitride for large strain, high-strain rate, and highpressure applications. Int J Impact Eng 2001;25:211–31. [16] Johnson GR, Holmquist TJ, Beissel SR. Response of aluminum nitride (including a phase change) to large strains, high strain rates, and high pressures. J Appl Phys 2003;94(3):1639–46. [17] Zhang YQ, Lu Y, Hao H. Analysis of fragment size and ejection velocity at high strain rate. Int J Mech Sci 2004;46:27–34. [18] Lu Y, Xu K. Modeling of dynamic behavior of concrete materials under blast loading. Int J Solids Struct 2004;41:131–43. [19] Munjiza A, Owen RJ, Bicanic N. A combined finite–discrete element method in transient dynamics of fracturing solids. Eng Comput 1995; 12:145–74. [20] Johnson GR, Holmquist TJ. An improved computational constitutive model for brittle materials. In: Schmidt SC, Shaner JW, Samara GA, Ross M, editors. High pressure science and technology—1993. Woodbury, NY: AIP Press; 1994. p. 981–4. [21] LSTC. LS-DYNA theoretical manual, Livermore Software Technology Corporation, Livermore, CA, 1998. [22] LSTC. LS-DYNA keyword user’s manual, Version 970, Livermore Software Technology Corporation, Livermore, CA, 2003.
ARTICLE IN PRESS G.W. Ma, X.M. An / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 966–975 [23] Hoek E, Brown ET. Underground excavations in rock. London: Institute of Mining and Metallurgy; 1980. [24] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci 1997;34(8):1165–86. [25] Cho SH, Kaneko K. Influence of the applied pressure waveform on the dynamic fracture processes in rock. Int J Rock Mech Min Sci 2004;41:771–84. [26] Cho SH, Ogata Y, Kaneko K. Strain rate dependency of the dynamic tensile strength of rock. Int J Rock Mech Min Sci 2003;40:763–77.
975
[27] Hagan TN. Rock breakage by explosives. Acta Astron 1979;6: 329–40. [28] Zhu ZM, Mohanty B, Xie HP. Numerical investigation of blastinginduced crack initiation and propagation in rocks. Int J Rock Mech Min Sci 2007;44:412–24. [29] Mohanty B. Explosion generated fractures in rock and rock-like materials. Eng Fract Mech 1990;35(4/5):889–98. [30] Fourney WL, Dally JW, Holloway DC. Controlled blasting with ligamented charge holders. Int J Rock Mech Min Sci 1978;15:121–9.