- Email: [email protected]
Modelling of blast-induced fractures in jointed rock masses
Engineering Fracture Mechanics 76 (2009) 1945–1955
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Modelling of blast-induced fractures in jointed rock masses Z.L. Wang a,b,*, H. Konietzky b a b
Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China Institut für Geotechnik, Technische Universität Bergakademie Freiberg, Freiberg 09596, Germany
a r t i c l e
i n f o
Article history: Received 25 August 2008 Received in revised form 16 March 2009 Accepted 9 May 2009 Available online 18 May 2009 Keywords: Jointed rock masses Blasting Coupled method Dynamic fracture Numerical modelling
a b s t r a c t This paper focuses on the dynamic fracturing process of jointed rock masses due to blast wave loading. A coupled numerical method using both LS-DYNA (a transient dynamic finite element program) and UDEC (an universal discrete element code) is outlined first. The blast-induced fracture extension in two representative jointed rock masses is studied subsequently. Besides, the effects of loading density of explosive charge on fracture evolution and the effects of the preexisting earth stress and free face on the fracturing process of rock masses are also explored in detail. The numerical results capture some of the well-known phenomena observed in the field and expected theoretically, and thus can offer useful guidelines to blasting design in rock masses. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction The destruction of hard rocks by means of blasting usually involves the drilling of a borehole, placement of an explosive charge and stemming prior to detonation. When the explosive is detonated, an extremely high pressure pulse is generated which is transmitted into rock mass adjacent to the borehole, producing a dilatational wave that propagates away from the charge. This may cause damage to rock and, furthermore, when the compressive stress wave reaches a free face or fissure (non-transmission), it will be reflected and converted into tensile wave, which may produce tensile cracking or cause spalling of surficial slabs if the tensile strength of the rock is exceeded [1–3]. Some researchers [4,5] believe that cracking is mainly caused by the incident dilatational wave and any reflected waves, while other investigators [6] consider the action of the compressed gases forcing its way through the cracks from the borehole more important. Until recently, it is generally agreed that both stress wave and gas pressure loadings play an important role in the process of rock fracture and fragmentation. In fact, our understanding of the blasting process is far from thorough, as both the explosive and the rock are complex materials. The release characteristics of energy stored in the explosive charge are highly variable, depending on the predominant factors like borehole parameters and charge loading density. Similarly, the response of the circumambient rock mass to high strain-rate dynamic loadings, which may last only for a few milliseconds, remains largely unknown. Under this scenario, the best approach to the study may be to first generate an extensive experimental database on these properties. At the same time, it is necessary to explore the fracture and fragmentation processes through numerical tools so as to obtain a better understanding of the underlying mechanism [3,7]. It is well recognized that the properties of a rock mass are determined by the properties of the intact rock and the discontinuities such as joint structures and faults [8,9]. The presence of the discontinuities has significant influence on the responses of the rock mass to either static or dynamic loading, and renders the numerical simulations more complicated [9,10]. * Corresponding author. Address: Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China. Fax: +86 551 3606459. E-mail address: [email protected] (Z.L. Wang). 0013-7944/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2009.05.004
1946
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Nomenclature
Drn Dun kn ks
normal stress increment normal displacement increment normal stiffness shear stiffness rn normal stress ss shear stress rt limiting tensile strength smax maximum shear stress C joint cohesion joint friction angle / w joint dilation angle Dues elastic incremental shear displacement Dus total incremental shear displacement critical shear displacement ucs ql loading density of charge q0 mass density of rock E young’s modulus of rock v poisson’s ratio of rock C1, r1, C2, r2, x constants of explosive P pressure V specific volume e specific energy initial specific energy e0 k wavelength Dl spatial element size t time horizontal effective stress Sh vertical effective stress Sv
There are some numerical codes available for simulating discontinuous rock mass at present, the popular ones being the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM). For the joints, they can be modelled as a kind of special ‘‘interface elements’’ in FEM [11], or simulated using ‘‘slide lines” in FDM, etc [12]. However, these numerical treatments are incapable of dealing with many intersecting interfaces and have no an automatic scheme for recognizing new contacts. In addition, they are limited to small displacements and/or rotation, etc. UDEC, a two-dimensional discrete element code, is proven as an ideal numerical technique for modelling the jointor fault-related problem in discontinuous rock system [13,14]. The rock mass is represented as an assemblage of discrete blocks, and the internal discontinuities are treated as boundary conditions between blocks. Large displacements along discontinuities and rotations of blocks are allowed. Individual blocks behave as either rigid or deformable materials. Deformable blocks are subdivided into constant-strain elements or zones, and each zone follows a prescribed stress-strain law. In addition, an appropriate material model must also be assigned to all discontinuities in calculations. Generally, the area contact Coulomb slip model is the most commonly used model for engineering studies [14]. This study concentrates on the role of stress-wave induced fracture in jointed rock masses because the rock fracturing process under stress wave loading is considered the crucial stage as the subsequent rock fragmentation and large-scale movement due to gas pressure loading [7]. After the coupled method combining LS-DYNA and UDEC is introduced in detail, the effects of joint structure, loading density, earth stress as well as free face on the fracture creation and evolution are explored by assuming that the blocks behave elastically and the joints follow the Coulomb slip behavior. 2. Numerical simulation tools In the present study, the numerical analyses are carried out using a coupled method which combines the LS-DYNA and the UDEC together. The former generates the blast loading which is used by the latter to simulate the blast-wave propagation in jointed rock mass. LS-DYNA is a general-purpose commercial code, which can be applied to solve a wide variety of non-linear problems in solid, fluid and gas dynamics [15]. This code has been proved by numerous studies as a suitable tool for large deformation analysis including the rock fracturing [3,16]. It performs dynamic analysis by seeking a solution to the momentum equation, which satisfies all boundary conditions, while integrating the energy equation to be used as a balance for global energy. Kinematic boundary conditions are applied when the equations of motion and constitutive equations are solved. Time integration is conducted via the central difference method when solving for displacements, velocities and accelerations. The timestep is
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1947
determined by the smallest element in the entire model. Databases are then generated which record history variables such as stress, strain and kinetic energy. Finally, LS-DYNA updates the current time and checks it against the simulation termination time. It has been well known that handling discontinuities in rock masses is a challenging task. UDEC is specially developed to model discontinuous problems [14]. It can accommodate a large number of discontinuities and permits the modelling system undergoing larger geometrical change through the use of contact updating scheme. In UDEC, the deformation of a fractured rock mass consists of the elastic/plastic deformation of blocks of intact rock, together with the displacements along and across fractures. The motion of block is characterized by the Newton’s second law of motion, which is expressed in the central finite difference form with respect to time. Calculations are performed over one timestep in an explicit time-marching algorithm. For deformable blocks, numerical integration of the differential equation of motion is used to determine the incremental displacements at the gridpoints of the triangular constant strain element within the blocks. The incremental displacements are then used to calculate the new stresses within the element through an appropriate constitutive equation. The amount of normal and tangential displacement between two adjacent blocks can be determined directly from the block geometry and the translation and rotation of the centroid of each block, as will be described in the upcoming sections.
3. Coupled method and joint model In this study, the coupled method combining LS-DYNA and UDEC is employed to simulate the blast-wave propagation and associated dynamic fracturing in jointed rock masses. An illustration of the above-mentioned coupled method (refer to Fig. 1) is elaborated as follows: Firstly, LS-DYNA is used to model the explosion process of the explosive whereby the rock mass is temporarily assumed as a continuous medium. The explosion history is measured on the wall of the borehole (or chamber) in terms of particle velocity or pressure (see Fig. 1b). Secondly, after the explosion history is converted into the formatted data compatible with UDEC, it is fed into the UDEC modelling as a radial velocity history. Thirdly, the blast-wave propagation and fracture evolution in the jointed rock mass, induced by detonating the explosive charge, is simulated using the UDEC modelling (see Fig. 1c). The intact rock blocks are assumed to be perfectly elastic in the present study. As for the joint behavior model, many types of constitutive laws may be contemplated. It has been proven that the Coulomb slip model with residual strength (area contact) is capable of capturing several of the features which are representative of the physical response of joints [14] and is thus adopted in the present study. In the normal direction, the stress-displacement relation is assumed to be linear and governed by the stiffness kn (expressed in terms of incremental stress per unit of displacement) as
Drn ¼ kn Dun
ð1Þ
where Drn is the normal stress increment of the interface between blocks, Dun is the normal displacement increment. Also, there is a limiting tensile strength, rt, for the joint. If the tensile strength is exceeded, i.e.,
rn < rt
ð2Þ
then let rn = 0. Similarly, in the shear direction, the response is controlled by a constant interface shear stiffness, ks. The shear stress ss, is limited by a combination of joint cohesion (C) and friction angle (/). Thus, if
ss 6 smax ¼ C þ rn tan /
ð3Þ
Dss ¼ ks Dues
ð4Þ
then,
or else, if
ss P smax
ð5Þ
ss ¼ signðDus Þsmax
ð6Þ
then,
where Dues is the elastic component of the incremental shear displacement, Dus is the total incremental shear displacement. In addition, joint dilation may occur at the onset of slip of the joint, the dilation is restricted such that [14]:
w¼0
if if
ss 6 smax ss ¼ smax and jus j P ucs
where w is the dilation angle, ucs denotes the critical shear displacement.
ð7Þ
1948
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Fig. 1. Illustration of the coupled method.
Based on this joint model, a fracture is considered to be initiated when the joint between two elements exceeds its tensile or shear strength, and the fracture begins to propagate when new discontinuities are created at its tip.
1949
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
4. Dynamic fracturing processes and discussions Drilling and blasting has been widely used in underground excavation and construction over the years and it is still a popular method of rock breakage. In this regard, a clear description of the fracturing process in a rock mass is highly desirable [17]. In this section, the influences of joint pattern, loading density, stress field and free face on initiation and propagation of the explosion-generated fractures are explored and discussed in detail. The problem involves an explosion of a cylindrical explosive column in a borehole with a radius of 30 cm (see Fig. 2a). The radius of the explosive column can be greater than zero and less than or equal to 30 cm, thus corresponding to different loading densities. Fig. 2b shows the UDEC modelling by cutting a horizontal cross-section at the charge level in Fig. 2a. The crosssection was set as 5.0 m10.0 m in calculations which can be simulated as a plane-strain problem. The non-reflecting boundaries were given at the top, bottom and right sides, and the left side was treated as symmetric boundary. The explosion input to the UDEC modelling, in terms of the velocity time history at the wall of borehole, was computed via the LS-DYNA modelling. The rock blocks were assumed to behave elastically with the following typical properties: mass density q0 = 2650.0 kg/ m3, Young’s modulus E = 90.0GPa, Poisson’s ration v = 0.16. Coulomb slip behavior was assumed for the joints. Table 1 lists the typical textbook values for joint properties. The Jones-Wilkens-Lee (JWL) equation of state was used to model the pressure generated by the expansion of the detonation product of the chemical explosive. It has been widely used in engineering calculations [18] and can be written as
x r1 V x r2 V xe e e P ¼ C1 1 þ C2 1 þ V r1 V r2 V
ð8Þ
where C1, r1, C2, r2 and x are material constants, P is the pressure, V is the specific volume, and e specific energy with an initial value of e0. The JWL parameters for the explosives (TNT) used in this work are presented in Table 2.
Fig. 2. Computation model for underground blasting.
Table 1 Material properties of joint. kn (GPa/m)
ks (GPa/m)
C (MPa)
/ (°)
w (°)
rt (MPa)
/r (°)
300.0
300.0
100.0
30.0
15.0
20.0
19.0
/r, the residual friction angle.
1950
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Table 2 JWL Parameters for TNT explosive.
q0 (kg/m3)
C1 (GPa)
C2 (GPa)
DC–J (m/s)
r1
r2
x
e0 (GPa)
1640.0
373.77
3.23
6930.0
4.15
0.95
0.30
7.0
DC–J, the C–J detonation velocity.
Fig. 3 plots the blast loading (x-velocity at the borehole wall) computed by LS-DYNA. Due to this measured point (x = 0.3 m, y = 0 m) as marked in Fig. 2a on the x-axis, the obtained explosion history can be regarded as the radial velocity. On the other hand, since the explosion history from the LS-DYNA is obtained randomly with time, it must be converted into a compatible form (see, for example, equal time-spacing) before feeding to the gridpoints on the borehole boundary in the UDEC modelling (see Fig. 2b). The damping in the numerical simulations should attempt to reproduce the energy loss in the system when subjected to dynamic loading [14]. Rayleigh damping is commonly used to provide damping that is approximately frequency-independent over a restricted range of frequencies. In order to find the principal frequency of the input velocity history (refer to Fig. 3), a spectral analysis should be made before simulations. As can be seen from Fig. 4, the principal frequency of the blast wave (about 1657.0 Hz) is obtained via fast Fourier transform algorithm.
20
Velocity (m/s)
15
10
5
0
-5
0
2
4 6 Time (ms)
8
10
Fig. 3. Explosion history calculated by LS-DYNA (ql = 410.0 kg/m3).
0.2
Fourier amplitude
0.15
0.1
0.05 Principal frequency
0
0
500
1000
1500 2000 2500 Frequency (Hz)
3000
Fig. 4. Power spectrum of input velocity–time history.
3500
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1951
Damping parameters are very important for UDEC dynamic analyses. It also relates to the maximum edge length of triangle-shaped finite difference zones of UDEC modelling. Some researches show that for accurate representation of wave transmission through a model, the spatial element size Dl must be smaller than approximately on-tenth to one-eighth of the wavelength k associated with the highest frequency component that contains appreciable energy [14], i.e.,
Dl 6
1 1 k 10 8
ð9Þ
Again, in order to satisfy the comparability, each block in the UDEC modelling was subdivided into constant-strain finitedifference triangular zones of the equal size as the elements in the LS-DYNA modelling. Preliminary analyses show that the adoption of a zone size of 0.15 m for the present problem is small enough to allow waves to propagate accurately at the input frequency without distortion. Another important item in a dynamic analysis is the reproduction of frequency-independent damping of materials at the correct level. For geomedia (soil and rock), the natural damping commonly falls in the range of 2.0–5.0% of critical value [14,17]. Herein, the fraction of critical damping of 0.02 is adopted for the present study. 4.1. Effect of joint structure Two joint patterns are shown in Fig. 5. One is randomly polygonal joints (namely Voronoi joints) while the other is the two sets of orthogonally aligned joints dipping at 15° and 75°, respectively. These two representative joints are used to examine the effect of blast-induced fracture creation and evolution. When the converted dynamic input is applied at the semi-circular boundary of radius 0.3 m (see Fig. 2b), the dynamic fracturing process of rock can be simulated using UDEC. As shown in Fig. 6, the radial fractures are first initiated and then propagate from the borehole surface to exterior boundaries. With time lapse, the fractures continuously evolve. At t = 5.0 ms, the fracture extension comes to an end and the maximum length approximates to 2.5 m. The radius of the fractured zone was found to approach 7.0–9.0 times the borehole radius, which also conforms to some established results for a cylindrical charge situation [19]. The similar phenomenon can be observed in the case of oriented joints (see Fig. 7). The fracture creation and propagation in this second case tends to be of lower consistency or higher anisotropy compared to the case of the previous Voronoi joint pattern. It can be seen that Voronoi joint generator is more illustrative in demonstrating the fracture propagation. Fracturing spontaneously occurs when the joint strength between neighboring blocks is exceeded. As such, the Voronoi jointed rock mass is focused in the following further numerical explorations. 4.2. Effect of loading density The dynamic responses of a rock mass to blasting are much affected by loading density of the explosive charge, which may influence the rock fracture pattern. Three loading densities ql of 182.2 kg/m3, 410.0 kg/m3 and 728.9 kg/m3 (the corresponding radii of explosive column are 10 cm, 15 cm and 20 cm, respectively) were adopted in the numerical simulations. Figs. 3 and 8 show the velocity time histories for the three different cases. Clearly, the larger the loading density, the higher the peak value of compressive wave is. Also, it can be found that the peak values of the tensile waves (negative) has the same changing tendency.
(a) Voronoi joints
(b) Orthogonal joints
Fig. 5. Joint patterns in calculation domain.
1952
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
(a) t=0.5ms
(b) t=1.0ms
(c) t=5.0ms
Fig. 6. Fracture evolution in Voronoi jointed rock mass.
(a) t=0.5ms
(b) t=1.0ms
(c) t=5.0ms
Fig. 7. Fracture evolution in orthogonally jointed rock mass.
Fig. 9 shows the effect of loading density on the fracture pattern of rock mass at the same time (t = 5 ms). A higher loading density of charge (728.9 kg/m3) increases the number of fractures and causes the intense stress release around the running fractures. The stress release induced by adjacent fractures disperses the concentrated stresses around the tips of the following fractures and thus affects fracture propagation, resulting in many shorter factures [20]. When the loading density decreases to lower level (410.0 kg/m3), the fracturing due to tension and shear loads markedly reduces and there are comparatively longer and continuous fracture extension (see Fig. 9b). When the loading density of charge decreases further to 182.2 kg/m3, only fewer fractures as shown in Fig. 9a can be observed. 4.3. Effect of preexisting stress field In deep underground operations, blasting can be affected by earth stress. High earth stress field may influence dynamic fracture process and orientation around the borehole. For example, the rock fractures may be extended or suppressed by preexisting stress field. In general, the high earth stress field can induce stress concentrations around the wall of the borehole and further lead to non-uniformity in the fracture pattern. In the present numerical simulations, two opposite edges of the numerical modelling (ql = 728.9 kg/m3) are pre-loaded with three different stress fields (see Fig. 10). The minor horizontal effective stress Sh is set equal to 0.3 MPa, acting parallel to the x-axis. The major vertical effective stress Sv is set equal to 10MP, 30 MPa or 50 MPa for different models acting parallel
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
8
Velocity (m/s)
6
4
2
0
-2
0
2
4 6 Time (ms)
8
10
8
10
(a) ρ l =182.2kg/m3 25 20
Velocity (m/s)
15 10 5 0 -5 -10 -15
0
2
4 6 Time (ms)
(b) ρ l =728.9kg/m3 Fig. 8. Explosion histories of different loading densities.
(a)
ρ l =182.2 kg/m3
(b)
ρ l =410.0 kg/m3
(c)
ρ l =728.9 kg/m3
Fig. 9. Effect of loading density on rock fracturing (t = 5.0 ms).
1953
1954
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
(a) Sh =0.3Mpa, Sv =10MPa
(b) Sh =0.3Mpa, Sv =30MPa
(c) Sh =0.3Mpa, Sv =50MPa
Fig. 10. Effect of stress field on rock fracturing (t = 1.0 ms).
to the y-axis. Fig. 10 clearly demonstrates that the blast-induced fractures are aligned in the direction of the principal stress axis (vertical), while the crack extension in the horizontal direction is intensively suppressed under a higher Sv. The results are in good agreement with some reported numerical and experimental observations in the literature [16,21]. 4.4. Effect of nearby free face In practice, rock breakage events caused by applied borehole pressures are usually performed near free boundaries to increase the fracturing effect on rock masses. It is well known that when a compressive stress-wave reaches the free face of a rock mass, it will be reflected backward to produce a tensile wave. If the reflected tensile wave is sufficiently strong (higher than the dynamic tensile strength of rock), it may cause cracking and even spalling of surficial slabs [22]. Fig. 11 shows the effect of free face on rock fracturing at t = 2.0 ms. Three different numerical models with one free face or two free faces are presented. It can be noticed that fractures are first initiated and extend from the borehole surface without directional preference. With further propagation of the stress waves, the fractures caused by the tensile wave starts to appear near the non-transmission or reflecting boundary and are approximately parallel to the boundary. Some fractures interlace together and are extended to a considerable length during the spalling.
(a) Upper free-face
(b) Right free-face
(c) Right and lower free-face
Fig. 11. Effect of free face on rock fracturing (t = 2.0 ms).
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1955
5. Conclusions The dynamic fracturing process in rock blasting has received considerable attention over the last few decades. This paper uses the hybrid method of UDEC and LS-DYNA to model the underground explosion in jointed rock masses and to investigate the influences of the dominant factors, including joint pattern, loading density, stress field and free face, on the rock fracturing process. The following conclusions can be drawn from the present study: (1) The coupled method is a good numerical tool for simulating blasting in jointed rock masses. LS-DYNA is adopted to model the explosive detonation and to provide the blast loading to be used in UDEC, which in turn simulates the stress-wave propagation and rock fracturing process. (2) Blasting in rock masses with different joint structures will behave differently. In a rock mass containing orthogonal joints, the crack evolution exhibits strong anisotropy. The Voronoi joint generator is found to be very conducive for numerical simulation of rock crack creation and propagation due to blasting. The fracturing happens when the joint strength between Voronoi blocks is exceeded. (3) The loading density of charge has a significant effect on the rock fracture pattern. A higher ql implies a higher peak of dilatational wave, which induces a larger fractured zone with many shorter cracks. However, there are only fewer fractures when the loading density is lower. (4) The existence of stress field considerably affects the fracture extension and causes the non-uniformity of rock fracture. The fractures tend to evolve along the major stress axis. Besides, near the free face, the stress wave is reflected and converted into tensile wave, resulting in the well-known phenomenon of surficial spalling. (5) Because the discrete blocks of rock is considered as elastic media in our calculations, no plastic crushed zone can be observed in the immediate vicinity of borehole. Therefore, there is still room for further refinement in future research work.
Acknowledgements This study was supported by the Program for New Century Excellent Talents in University (No. NCET-08-0525) the Alexander von Humboldt Foundation, the Civil Aeronautics Joint Research Foundation (No. 60776821) and the Specialized Research Fund for the Doctoral Program of Higher Education (No. 20070358073). The first author would like to thank Dr. R.F. Shen for his embellishment of this paper. References [1] Langefors U, Kihlstrom B. The modern technique of rock blasting. 3rd ed. Stockholm, Sweden: Almqvist & Wiksell Forlag AB; 1978. [2] Dally JW, Fourney WL, Holloway DC. Influence of containment of the borehole pressure on explosive induced fracture. Int J Rock Mech Min Sci 1975;12:5–12. [3] Wang ZL, Li YC, Shen RF. Numerical simulation of tensile damage and blast crater in brittle rock due to underground explosion. Int J Rock Mech Min Sci 2007;44:730–8. [4] Liu LQ, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci Geomech Abstr 1997;34:217–31. [5] Atkinson BK. Fracture mechanics of rock. New York: Academic Press Geology Series; 1987. [6] Paine AS, Please CP. An improved model of fracture propagation by gas during rock blasting—some analytical results. Int J Rock Mech Min Sci 1994;31:699–706. [7] Zhu ZM, Mohanty B, Xie HP. Numerical investigation of blasting-induced crack initiation and propagation in rocks. Int J Rock Mech Min Sci 2007;44:412–24. [8] Goodman RE. Methods of geological engineering in discontinuous rocks. St. Paul: West Publishing; 1976. [9] Cundall PA. Numerical modelling of jointed and faulted rock. In: Rossmanith HP, editor. Mechanics of jointed and faulted rock. Rotterdam: AA Balkema; 1990. [10] Chen SG, Zhao J. A study of UDEC modelling for blast wave propagation in jointed rock masses. Int J Rock Mech Min Sci 1998;35:93–9. [11] Desai CS, Zaman MM, Lightner JG, Siriwardane HJ. Thin-layer element for interfaces and joints. Int J Numer Anal Methods Geomech 1984;8:19–43. [12] Coates RT, Schoenberg M. Finite difference modelling of faults and fractures. Geophysics 1995;60:1514–26. [13] Hart RD. An introduction to distinct element modelling for rock engineering. In: Hudson JA, editor. Comprehensive rock engineering. Oxford: Pergamon Press; 1993. [14] Itasca Consulting Group, Inc. UDEC-Universal Distinct Element Code. Version 3.0, User Manual. Minnesota, USA; 2004. [15] Hallquist JO. LS-DYNA theoretical manual. Livermore, CA: Livermore Software Technology Corporation; 2003. [16] Ma GW, An XM. Numerical simulation of blasting-induced rock fractures. Int J Rock Mech Min Sci 2008;45:966–75. [17] Wang ZL, Konietzky H, Shen RF. Coupled finite element and discrete element method for underground blast in faulted rock masses. Soil Dynam Earthq Eng 2009;29:939–45. [18] Livermore Software Technology Corporation (LSTC). LS-DYNA Keyword User’s manual. Livermore, CA; 2003. [19] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8:181–202. [20] 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. [21] Donze FV, Bouchez J, Magnier SA. Modelling fractures in rock blasting. Int J Rock Mech Mining Sci 1997;34:1153–63. [22] Wang ZL, Li YC, Wang JG. Numerical analysis of blast-induced wave propagation and spalling damage in a rock plate. Int J Rock Mech Min Sci 2008;45:600–8.
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Modelling of blast-induced fractures in jointed rock masses Z.L. Wang a,b,*, H. Konietzky b a b
Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China Institut für Geotechnik, Technische Universität Bergakademie Freiberg, Freiberg 09596, Germany
a r t i c l e
i n f o
Article history: Received 25 August 2008 Received in revised form 16 March 2009 Accepted 9 May 2009 Available online 18 May 2009 Keywords: Jointed rock masses Blasting Coupled method Dynamic fracture Numerical modelling
a b s t r a c t This paper focuses on the dynamic fracturing process of jointed rock masses due to blast wave loading. A coupled numerical method using both LS-DYNA (a transient dynamic finite element program) and UDEC (an universal discrete element code) is outlined first. The blast-induced fracture extension in two representative jointed rock masses is studied subsequently. Besides, the effects of loading density of explosive charge on fracture evolution and the effects of the preexisting earth stress and free face on the fracturing process of rock masses are also explored in detail. The numerical results capture some of the well-known phenomena observed in the field and expected theoretically, and thus can offer useful guidelines to blasting design in rock masses. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction The destruction of hard rocks by means of blasting usually involves the drilling of a borehole, placement of an explosive charge and stemming prior to detonation. When the explosive is detonated, an extremely high pressure pulse is generated which is transmitted into rock mass adjacent to the borehole, producing a dilatational wave that propagates away from the charge. This may cause damage to rock and, furthermore, when the compressive stress wave reaches a free face or fissure (non-transmission), it will be reflected and converted into tensile wave, which may produce tensile cracking or cause spalling of surficial slabs if the tensile strength of the rock is exceeded [1–3]. Some researchers [4,5] believe that cracking is mainly caused by the incident dilatational wave and any reflected waves, while other investigators [6] consider the action of the compressed gases forcing its way through the cracks from the borehole more important. Until recently, it is generally agreed that both stress wave and gas pressure loadings play an important role in the process of rock fracture and fragmentation. In fact, our understanding of the blasting process is far from thorough, as both the explosive and the rock are complex materials. The release characteristics of energy stored in the explosive charge are highly variable, depending on the predominant factors like borehole parameters and charge loading density. Similarly, the response of the circumambient rock mass to high strain-rate dynamic loadings, which may last only for a few milliseconds, remains largely unknown. Under this scenario, the best approach to the study may be to first generate an extensive experimental database on these properties. At the same time, it is necessary to explore the fracture and fragmentation processes through numerical tools so as to obtain a better understanding of the underlying mechanism [3,7]. It is well recognized that the properties of a rock mass are determined by the properties of the intact rock and the discontinuities such as joint structures and faults [8,9]. The presence of the discontinuities has significant influence on the responses of the rock mass to either static or dynamic loading, and renders the numerical simulations more complicated [9,10]. * Corresponding author. Address: Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China. Fax: +86 551 3606459. E-mail address: [email protected] (Z.L. Wang). 0013-7944/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2009.05.004
1946
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Nomenclature
Drn Dun kn ks
normal stress increment normal displacement increment normal stiffness shear stiffness rn normal stress ss shear stress rt limiting tensile strength smax maximum shear stress C joint cohesion joint friction angle / w joint dilation angle Dues elastic incremental shear displacement Dus total incremental shear displacement critical shear displacement ucs ql loading density of charge q0 mass density of rock E young’s modulus of rock v poisson’s ratio of rock C1, r1, C2, r2, x constants of explosive P pressure V specific volume e specific energy initial specific energy e0 k wavelength Dl spatial element size t time horizontal effective stress Sh vertical effective stress Sv
There are some numerical codes available for simulating discontinuous rock mass at present, the popular ones being the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM). For the joints, they can be modelled as a kind of special ‘‘interface elements’’ in FEM [11], or simulated using ‘‘slide lines” in FDM, etc [12]. However, these numerical treatments are incapable of dealing with many intersecting interfaces and have no an automatic scheme for recognizing new contacts. In addition, they are limited to small displacements and/or rotation, etc. UDEC, a two-dimensional discrete element code, is proven as an ideal numerical technique for modelling the jointor fault-related problem in discontinuous rock system [13,14]. The rock mass is represented as an assemblage of discrete blocks, and the internal discontinuities are treated as boundary conditions between blocks. Large displacements along discontinuities and rotations of blocks are allowed. Individual blocks behave as either rigid or deformable materials. Deformable blocks are subdivided into constant-strain elements or zones, and each zone follows a prescribed stress-strain law. In addition, an appropriate material model must also be assigned to all discontinuities in calculations. Generally, the area contact Coulomb slip model is the most commonly used model for engineering studies [14]. This study concentrates on the role of stress-wave induced fracture in jointed rock masses because the rock fracturing process under stress wave loading is considered the crucial stage as the subsequent rock fragmentation and large-scale movement due to gas pressure loading [7]. After the coupled method combining LS-DYNA and UDEC is introduced in detail, the effects of joint structure, loading density, earth stress as well as free face on the fracture creation and evolution are explored by assuming that the blocks behave elastically and the joints follow the Coulomb slip behavior. 2. Numerical simulation tools In the present study, the numerical analyses are carried out using a coupled method which combines the LS-DYNA and the UDEC together. The former generates the blast loading which is used by the latter to simulate the blast-wave propagation in jointed rock mass. LS-DYNA is a general-purpose commercial code, which can be applied to solve a wide variety of non-linear problems in solid, fluid and gas dynamics [15]. This code has been proved by numerous studies as a suitable tool for large deformation analysis including the rock fracturing [3,16]. It performs dynamic analysis by seeking a solution to the momentum equation, which satisfies all boundary conditions, while integrating the energy equation to be used as a balance for global energy. Kinematic boundary conditions are applied when the equations of motion and constitutive equations are solved. Time integration is conducted via the central difference method when solving for displacements, velocities and accelerations. The timestep is
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1947
determined by the smallest element in the entire model. Databases are then generated which record history variables such as stress, strain and kinetic energy. Finally, LS-DYNA updates the current time and checks it against the simulation termination time. It has been well known that handling discontinuities in rock masses is a challenging task. UDEC is specially developed to model discontinuous problems [14]. It can accommodate a large number of discontinuities and permits the modelling system undergoing larger geometrical change through the use of contact updating scheme. In UDEC, the deformation of a fractured rock mass consists of the elastic/plastic deformation of blocks of intact rock, together with the displacements along and across fractures. The motion of block is characterized by the Newton’s second law of motion, which is expressed in the central finite difference form with respect to time. Calculations are performed over one timestep in an explicit time-marching algorithm. For deformable blocks, numerical integration of the differential equation of motion is used to determine the incremental displacements at the gridpoints of the triangular constant strain element within the blocks. The incremental displacements are then used to calculate the new stresses within the element through an appropriate constitutive equation. The amount of normal and tangential displacement between two adjacent blocks can be determined directly from the block geometry and the translation and rotation of the centroid of each block, as will be described in the upcoming sections.
3. Coupled method and joint model In this study, the coupled method combining LS-DYNA and UDEC is employed to simulate the blast-wave propagation and associated dynamic fracturing in jointed rock masses. An illustration of the above-mentioned coupled method (refer to Fig. 1) is elaborated as follows: Firstly, LS-DYNA is used to model the explosion process of the explosive whereby the rock mass is temporarily assumed as a continuous medium. The explosion history is measured on the wall of the borehole (or chamber) in terms of particle velocity or pressure (see Fig. 1b). Secondly, after the explosion history is converted into the formatted data compatible with UDEC, it is fed into the UDEC modelling as a radial velocity history. Thirdly, the blast-wave propagation and fracture evolution in the jointed rock mass, induced by detonating the explosive charge, is simulated using the UDEC modelling (see Fig. 1c). The intact rock blocks are assumed to be perfectly elastic in the present study. As for the joint behavior model, many types of constitutive laws may be contemplated. It has been proven that the Coulomb slip model with residual strength (area contact) is capable of capturing several of the features which are representative of the physical response of joints [14] and is thus adopted in the present study. In the normal direction, the stress-displacement relation is assumed to be linear and governed by the stiffness kn (expressed in terms of incremental stress per unit of displacement) as
Drn ¼ kn Dun
ð1Þ
where Drn is the normal stress increment of the interface between blocks, Dun is the normal displacement increment. Also, there is a limiting tensile strength, rt, for the joint. If the tensile strength is exceeded, i.e.,
rn < rt
ð2Þ
then let rn = 0. Similarly, in the shear direction, the response is controlled by a constant interface shear stiffness, ks. The shear stress ss, is limited by a combination of joint cohesion (C) and friction angle (/). Thus, if
ss 6 smax ¼ C þ rn tan /
ð3Þ
Dss ¼ ks Dues
ð4Þ
then,
or else, if
ss P smax
ð5Þ
ss ¼ signðDus Þsmax
ð6Þ
then,
where Dues is the elastic component of the incremental shear displacement, Dus is the total incremental shear displacement. In addition, joint dilation may occur at the onset of slip of the joint, the dilation is restricted such that [14]:
w¼0
if if
ss 6 smax ss ¼ smax and jus j P ucs
where w is the dilation angle, ucs denotes the critical shear displacement.
ð7Þ
1948
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Fig. 1. Illustration of the coupled method.
Based on this joint model, a fracture is considered to be initiated when the joint between two elements exceeds its tensile or shear strength, and the fracture begins to propagate when new discontinuities are created at its tip.
1949
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
4. Dynamic fracturing processes and discussions Drilling and blasting has been widely used in underground excavation and construction over the years and it is still a popular method of rock breakage. In this regard, a clear description of the fracturing process in a rock mass is highly desirable [17]. In this section, the influences of joint pattern, loading density, stress field and free face on initiation and propagation of the explosion-generated fractures are explored and discussed in detail. The problem involves an explosion of a cylindrical explosive column in a borehole with a radius of 30 cm (see Fig. 2a). The radius of the explosive column can be greater than zero and less than or equal to 30 cm, thus corresponding to different loading densities. Fig. 2b shows the UDEC modelling by cutting a horizontal cross-section at the charge level in Fig. 2a. The crosssection was set as 5.0 m10.0 m in calculations which can be simulated as a plane-strain problem. The non-reflecting boundaries were given at the top, bottom and right sides, and the left side was treated as symmetric boundary. The explosion input to the UDEC modelling, in terms of the velocity time history at the wall of borehole, was computed via the LS-DYNA modelling. The rock blocks were assumed to behave elastically with the following typical properties: mass density q0 = 2650.0 kg/ m3, Young’s modulus E = 90.0GPa, Poisson’s ration v = 0.16. Coulomb slip behavior was assumed for the joints. Table 1 lists the typical textbook values for joint properties. The Jones-Wilkens-Lee (JWL) equation of state was used to model the pressure generated by the expansion of the detonation product of the chemical explosive. It has been widely used in engineering calculations [18] and can be written as
x r1 V x r2 V xe e e P ¼ C1 1 þ C2 1 þ V r1 V r2 V
ð8Þ
where C1, r1, C2, r2 and x are material constants, P is the pressure, V is the specific volume, and e specific energy with an initial value of e0. The JWL parameters for the explosives (TNT) used in this work are presented in Table 2.
Fig. 2. Computation model for underground blasting.
Table 1 Material properties of joint. kn (GPa/m)
ks (GPa/m)
C (MPa)
/ (°)
w (°)
rt (MPa)
/r (°)
300.0
300.0
100.0
30.0
15.0
20.0
19.0
/r, the residual friction angle.
1950
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
Table 2 JWL Parameters for TNT explosive.
q0 (kg/m3)
C1 (GPa)
C2 (GPa)
DC–J (m/s)
r1
r2
x
e0 (GPa)
1640.0
373.77
3.23
6930.0
4.15
0.95
0.30
7.0
DC–J, the C–J detonation velocity.
Fig. 3 plots the blast loading (x-velocity at the borehole wall) computed by LS-DYNA. Due to this measured point (x = 0.3 m, y = 0 m) as marked in Fig. 2a on the x-axis, the obtained explosion history can be regarded as the radial velocity. On the other hand, since the explosion history from the LS-DYNA is obtained randomly with time, it must be converted into a compatible form (see, for example, equal time-spacing) before feeding to the gridpoints on the borehole boundary in the UDEC modelling (see Fig. 2b). The damping in the numerical simulations should attempt to reproduce the energy loss in the system when subjected to dynamic loading [14]. Rayleigh damping is commonly used to provide damping that is approximately frequency-independent over a restricted range of frequencies. In order to find the principal frequency of the input velocity history (refer to Fig. 3), a spectral analysis should be made before simulations. As can be seen from Fig. 4, the principal frequency of the blast wave (about 1657.0 Hz) is obtained via fast Fourier transform algorithm.
20
Velocity (m/s)
15
10
5
0
-5
0
2
4 6 Time (ms)
8
10
Fig. 3. Explosion history calculated by LS-DYNA (ql = 410.0 kg/m3).
0.2
Fourier amplitude
0.15
0.1
0.05 Principal frequency
0
0
500
1000
1500 2000 2500 Frequency (Hz)
3000
Fig. 4. Power spectrum of input velocity–time history.
3500
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1951
Damping parameters are very important for UDEC dynamic analyses. It also relates to the maximum edge length of triangle-shaped finite difference zones of UDEC modelling. Some researches show that for accurate representation of wave transmission through a model, the spatial element size Dl must be smaller than approximately on-tenth to one-eighth of the wavelength k associated with the highest frequency component that contains appreciable energy [14], i.e.,
Dl 6
1 1 k 10 8
ð9Þ
Again, in order to satisfy the comparability, each block in the UDEC modelling was subdivided into constant-strain finitedifference triangular zones of the equal size as the elements in the LS-DYNA modelling. Preliminary analyses show that the adoption of a zone size of 0.15 m for the present problem is small enough to allow waves to propagate accurately at the input frequency without distortion. Another important item in a dynamic analysis is the reproduction of frequency-independent damping of materials at the correct level. For geomedia (soil and rock), the natural damping commonly falls in the range of 2.0–5.0% of critical value [14,17]. Herein, the fraction of critical damping of 0.02 is adopted for the present study. 4.1. Effect of joint structure Two joint patterns are shown in Fig. 5. One is randomly polygonal joints (namely Voronoi joints) while the other is the two sets of orthogonally aligned joints dipping at 15° and 75°, respectively. These two representative joints are used to examine the effect of blast-induced fracture creation and evolution. When the converted dynamic input is applied at the semi-circular boundary of radius 0.3 m (see Fig. 2b), the dynamic fracturing process of rock can be simulated using UDEC. As shown in Fig. 6, the radial fractures are first initiated and then propagate from the borehole surface to exterior boundaries. With time lapse, the fractures continuously evolve. At t = 5.0 ms, the fracture extension comes to an end and the maximum length approximates to 2.5 m. The radius of the fractured zone was found to approach 7.0–9.0 times the borehole radius, which also conforms to some established results for a cylindrical charge situation [19]. The similar phenomenon can be observed in the case of oriented joints (see Fig. 7). The fracture creation and propagation in this second case tends to be of lower consistency or higher anisotropy compared to the case of the previous Voronoi joint pattern. It can be seen that Voronoi joint generator is more illustrative in demonstrating the fracture propagation. Fracturing spontaneously occurs when the joint strength between neighboring blocks is exceeded. As such, the Voronoi jointed rock mass is focused in the following further numerical explorations. 4.2. Effect of loading density The dynamic responses of a rock mass to blasting are much affected by loading density of the explosive charge, which may influence the rock fracture pattern. Three loading densities ql of 182.2 kg/m3, 410.0 kg/m3 and 728.9 kg/m3 (the corresponding radii of explosive column are 10 cm, 15 cm and 20 cm, respectively) were adopted in the numerical simulations. Figs. 3 and 8 show the velocity time histories for the three different cases. Clearly, the larger the loading density, the higher the peak value of compressive wave is. Also, it can be found that the peak values of the tensile waves (negative) has the same changing tendency.
(a) Voronoi joints
(b) Orthogonal joints
Fig. 5. Joint patterns in calculation domain.
1952
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
(a) t=0.5ms
(b) t=1.0ms
(c) t=5.0ms
Fig. 6. Fracture evolution in Voronoi jointed rock mass.
(a) t=0.5ms
(b) t=1.0ms
(c) t=5.0ms
Fig. 7. Fracture evolution in orthogonally jointed rock mass.
Fig. 9 shows the effect of loading density on the fracture pattern of rock mass at the same time (t = 5 ms). A higher loading density of charge (728.9 kg/m3) increases the number of fractures and causes the intense stress release around the running fractures. The stress release induced by adjacent fractures disperses the concentrated stresses around the tips of the following fractures and thus affects fracture propagation, resulting in many shorter factures [20]. When the loading density decreases to lower level (410.0 kg/m3), the fracturing due to tension and shear loads markedly reduces and there are comparatively longer and continuous fracture extension (see Fig. 9b). When the loading density of charge decreases further to 182.2 kg/m3, only fewer fractures as shown in Fig. 9a can be observed. 4.3. Effect of preexisting stress field In deep underground operations, blasting can be affected by earth stress. High earth stress field may influence dynamic fracture process and orientation around the borehole. For example, the rock fractures may be extended or suppressed by preexisting stress field. In general, the high earth stress field can induce stress concentrations around the wall of the borehole and further lead to non-uniformity in the fracture pattern. In the present numerical simulations, two opposite edges of the numerical modelling (ql = 728.9 kg/m3) are pre-loaded with three different stress fields (see Fig. 10). The minor horizontal effective stress Sh is set equal to 0.3 MPa, acting parallel to the x-axis. The major vertical effective stress Sv is set equal to 10MP, 30 MPa or 50 MPa for different models acting parallel
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
8
Velocity (m/s)
6
4
2
0
-2
0
2
4 6 Time (ms)
8
10
8
10
(a) ρ l =182.2kg/m3 25 20
Velocity (m/s)
15 10 5 0 -5 -10 -15
0
2
4 6 Time (ms)
(b) ρ l =728.9kg/m3 Fig. 8. Explosion histories of different loading densities.
(a)
ρ l =182.2 kg/m3
(b)
ρ l =410.0 kg/m3
(c)
ρ l =728.9 kg/m3
Fig. 9. Effect of loading density on rock fracturing (t = 5.0 ms).
1953
1954
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
(a) Sh =0.3Mpa, Sv =10MPa
(b) Sh =0.3Mpa, Sv =30MPa
(c) Sh =0.3Mpa, Sv =50MPa
Fig. 10. Effect of stress field on rock fracturing (t = 1.0 ms).
to the y-axis. Fig. 10 clearly demonstrates that the blast-induced fractures are aligned in the direction of the principal stress axis (vertical), while the crack extension in the horizontal direction is intensively suppressed under a higher Sv. The results are in good agreement with some reported numerical and experimental observations in the literature [16,21]. 4.4. Effect of nearby free face In practice, rock breakage events caused by applied borehole pressures are usually performed near free boundaries to increase the fracturing effect on rock masses. It is well known that when a compressive stress-wave reaches the free face of a rock mass, it will be reflected backward to produce a tensile wave. If the reflected tensile wave is sufficiently strong (higher than the dynamic tensile strength of rock), it may cause cracking and even spalling of surficial slabs [22]. Fig. 11 shows the effect of free face on rock fracturing at t = 2.0 ms. Three different numerical models with one free face or two free faces are presented. It can be noticed that fractures are first initiated and extend from the borehole surface without directional preference. With further propagation of the stress waves, the fractures caused by the tensile wave starts to appear near the non-transmission or reflecting boundary and are approximately parallel to the boundary. Some fractures interlace together and are extended to a considerable length during the spalling.
(a) Upper free-face
(b) Right free-face
(c) Right and lower free-face
Fig. 11. Effect of free face on rock fracturing (t = 2.0 ms).
Z.L. Wang, H. Konietzky / Engineering Fracture Mechanics 76 (2009) 1945–1955
1955
5. Conclusions The dynamic fracturing process in rock blasting has received considerable attention over the last few decades. This paper uses the hybrid method of UDEC and LS-DYNA to model the underground explosion in jointed rock masses and to investigate the influences of the dominant factors, including joint pattern, loading density, stress field and free face, on the rock fracturing process. The following conclusions can be drawn from the present study: (1) The coupled method is a good numerical tool for simulating blasting in jointed rock masses. LS-DYNA is adopted to model the explosive detonation and to provide the blast loading to be used in UDEC, which in turn simulates the stress-wave propagation and rock fracturing process. (2) Blasting in rock masses with different joint structures will behave differently. In a rock mass containing orthogonal joints, the crack evolution exhibits strong anisotropy. The Voronoi joint generator is found to be very conducive for numerical simulation of rock crack creation and propagation due to blasting. The fracturing happens when the joint strength between Voronoi blocks is exceeded. (3) The loading density of charge has a significant effect on the rock fracture pattern. A higher ql implies a higher peak of dilatational wave, which induces a larger fractured zone with many shorter cracks. However, there are only fewer fractures when the loading density is lower. (4) The existence of stress field considerably affects the fracture extension and causes the non-uniformity of rock fracture. The fractures tend to evolve along the major stress axis. Besides, near the free face, the stress wave is reflected and converted into tensile wave, resulting in the well-known phenomenon of surficial spalling. (5) Because the discrete blocks of rock is considered as elastic media in our calculations, no plastic crushed zone can be observed in the immediate vicinity of borehole. Therefore, there is still room for further refinement in future research work.
Acknowledgements This study was supported by the Program for New Century Excellent Talents in University (No. NCET-08-0525) the Alexander von Humboldt Foundation, the Civil Aeronautics Joint Research Foundation (No. 60776821) and the Specialized Research Fund for the Doctoral Program of Higher Education (No. 20070358073). The first author would like to thank Dr. R.F. Shen for his embellishment of this paper. References [1] Langefors U, Kihlstrom B. The modern technique of rock blasting. 3rd ed. Stockholm, Sweden: Almqvist & Wiksell Forlag AB; 1978. [2] Dally JW, Fourney WL, Holloway DC. Influence of containment of the borehole pressure on explosive induced fracture. Int J Rock Mech Min Sci 1975;12:5–12. [3] Wang ZL, Li YC, Shen RF. Numerical simulation of tensile damage and blast crater in brittle rock due to underground explosion. Int J Rock Mech Min Sci 2007;44:730–8. [4] Liu LQ, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci Geomech Abstr 1997;34:217–31. [5] Atkinson BK. Fracture mechanics of rock. New York: Academic Press Geology Series; 1987. [6] Paine AS, Please CP. An improved model of fracture propagation by gas during rock blasting—some analytical results. Int J Rock Mech Min Sci 1994;31:699–706. [7] Zhu ZM, Mohanty B, Xie HP. Numerical investigation of blasting-induced crack initiation and propagation in rocks. Int J Rock Mech Min Sci 2007;44:412–24. [8] Goodman RE. Methods of geological engineering in discontinuous rocks. St. Paul: West Publishing; 1976. [9] Cundall PA. Numerical modelling of jointed and faulted rock. In: Rossmanith HP, editor. Mechanics of jointed and faulted rock. Rotterdam: AA Balkema; 1990. [10] Chen SG, Zhao J. A study of UDEC modelling for blast wave propagation in jointed rock masses. Int J Rock Mech Min Sci 1998;35:93–9. [11] Desai CS, Zaman MM, Lightner JG, Siriwardane HJ. Thin-layer element for interfaces and joints. Int J Numer Anal Methods Geomech 1984;8:19–43. [12] Coates RT, Schoenberg M. Finite difference modelling of faults and fractures. Geophysics 1995;60:1514–26. [13] Hart RD. An introduction to distinct element modelling for rock engineering. In: Hudson JA, editor. Comprehensive rock engineering. Oxford: Pergamon Press; 1993. [14] Itasca Consulting Group, Inc. UDEC-Universal Distinct Element Code. Version 3.0, User Manual. Minnesota, USA; 2004. [15] Hallquist JO. LS-DYNA theoretical manual. Livermore, CA: Livermore Software Technology Corporation; 2003. [16] Ma GW, An XM. Numerical simulation of blasting-induced rock fractures. Int J Rock Mech Min Sci 2008;45:966–75. [17] Wang ZL, Konietzky H, Shen RF. Coupled finite element and discrete element method for underground blast in faulted rock masses. Soil Dynam Earthq Eng 2009;29:939–45. [18] Livermore Software Technology Corporation (LSTC). LS-DYNA Keyword User’s manual. Livermore, CA; 2003. [19] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8:181–202. [20] 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. [21] Donze FV, Bouchez J, Magnier SA. Modelling fractures in rock blasting. Int J Rock Mech Mining Sci 1997;34:1153–63. [22] Wang ZL, Li YC, Wang JG. Numerical analysis of blast-induced wave propagation and spalling damage in a rock plate. Int J Rock Mech Min Sci 2008;45:600–8.