Hybrid finite-discrete element modelling of dynamic fracture and resultant fragment casting and muck-piling by rock blast

Hybrid finite-discrete element modelling of dynamic fracture and resultant fragment casting and muck-piling by rock blast

Computers and Geotechnics 81 (2017) 322–345 Contents lists available at ScienceDirect Computers and Geotechnics journal homepage: www.elsevier.com/l...

5MB Sizes 0 Downloads 25 Views

Computers and Geotechnics 81 (2017) 322–345

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Hybrid finite-discrete element modelling of dynamic fracture and resultant fragment casting and muck-piling by rock blast H.M. An a,b, H.Y. Liu a,⇑, Haoyu Han a, Xin Zheng a,c, X.G. Wang b a

School of Engineering and ICT, University of Tasmania, TAS 7001, Australia School of Civil & Environmental Engineering, University of Science and Technology Beijing, China c School of Resources and Civil Engineering, Northeastern University, Liaoning, China b

a r t i c l e

i n f o

Article history: Received 5 July 2016 Received in revised form 13 September 2016 Accepted 20 September 2016

Keywords: Hybrid FEM-DEM Dynamic fracture Rock blast Numerical method Failure mechanism Fragment muck-piling

a b s t r a c t A state-of-the-art review is conducted to highlight the fracture mechanism in rock blast and advantages and limitations of various methods in modelling it. A hybrid finite-discrete element method (FEM-DEM) is implemented to simulate rock fracture and resultant fragment muck-piling in various blasting scenarios. The modelled crushed, cracked and long radial crack zones are compared with those in literatures to calibrate the hybrid FEM-DEM. Moreover, the hybrid modelling reproduces the rock fragmentation process during blasting. It is concluded that the hybrid FEM-DEM is superior to continuous and discontinuous methods in terms of modelling dynamic fracture of rock under blast-induced impact load. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Blasting has been widely employed in the mining industry for many centuries and it remains a popular method of rock fragmentation, hard rock tunnelling and structure demolition in modern mining and civil engineerings. In mining engineering, the rock fragmentation by blasting is the first stage of the comminution process in mines and is the activity that may have the most leverage in the efficiency of a mining operation, as the output from a blast impacts every downstream operation. An improved fragmentation associated with blasting can result in efficient rock breakage, reduced costs for both secondary fragmentation and transportation of the blasted rock, effective destressing of rockburst prone area, improved environmental aspects, and reductions in energy consumption during crushing and grinding of the ore, as well as improved metal recovery. In civil engineering, controlling both fragmentation and the degree of blast induced damage in the hard rock tunnelling and structure demolition are important aspects of the project design process. Poor blasting practices are typified by excessive damage and over-break, oversize fragmentation, restricted access, increased local reinforcement requirements and increased project cycle times and costs. Therefore, it is essential ⇑ Corresponding author. E-mail address: [email protected] (H.Y. Liu). http://dx.doi.org/10.1016/j.compgeo.2016.09.007 0266-352X/Ó 2016 Elsevier Ltd. All rights reserved.

to study the rock fragmentation by blasting with an aim of whether inducing considerable fracture and fragmentation of rock or preventing failure of rock engineering structure under blast loads. 2. Review of previous studies on rock fragmentation process by blasting Considerable efforts were made in the last four decades to understand the rock fracturing mechanism in rock blasting and many studies concluded that the current situation in rock blasting is far from the theoretical optimum fragmentation [1–5]. The potential for improving blasting in mining, civil, petroleum and defence engineerings is huge and the economic potential is enormous. Blast fragmentation modelling is an important step to achieve the optimum fragmentation, which allows the estimation of blast fragmentation distributions for a number of different rock masses, blast geometries and explosive parameters. In the past, empirical models were put forward by some researchers [6–8]. These models were mainly developed based on the attributes observed after rock blasting, especially the size distribution of rock fragments. In 1930s, Rossin and Rammler [9] proposed the RosinRammler equation to characterize the particle-size distribution of material. Later, Kuznetsov [7] developed a semi-empirical equation for estimating the size distribution of rock fragments. The Kuz-Ram model [8,10] is probably the most widely used approach and the

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

modified Kuz-Ram model was used by Gheibie et al. [11] at Sungun Mine to predict the fines produced during blasting. The other important empirical models include the JKMRC blast fragmentation models [3], e.g. the Two-Component Model [12] and the Crushed Zone Model [13], and the Swebrec functions [14]. Since the empirical models require only a few of input parameters for engineering applications, they can be easily applied in routine blast design layout spreadsheets [8,10]. However, limited input parameters may lead to inaccurate prediction. In addition, in regards to the practice, several blastings must be tested frequently in advance, which makes the implementation of the empirical models expensive and time consuming. Moreover, the empirical model may not be able to satisfy the requirements of the modern rock blasting engineering since rock fracture and fragmentation progresses are ignored in the empirical model. As a matter of fact, the rock blasting is an extremely complex process and, generally, involves explosive detonation, gas expansion, stress wave propagation, rock fracturing and resultant rock fragment throwing and muck-piling. The lack of understanding of the complex process of the rock blasting has limited engineers to optimize rock blast design. Thus, this paper is intended to study the rock blasting process from the mechanics point of view. 2.1. Review of dynamic rock fracture mechanism during rock blasting Many researchers have conducted experimental and theoretical studies of the rock blasting processes in order to understand the rock fragmentation mechanism and then improve rock blasting efficiency. Latham et al. [15] charted the researches on the understanding of rock blasting processes drawing upon the work of researchers worldwide in the 1980’s and 1990’s. More recently, Saharan et al. [16] presented a detailed state-of-the-art review on the study of the dynamic rock fracture initiation and propagation due to explosive energy. In their review, explosive energy dissipation in crushing and fracturing was examined, the various means to enhance the explosive energy utilization for dynamic rock fracturing were reviewed, and the need for a better understanding of the dynamic fracturing process was particularly highlighted. According to these reviews, the different rock breaking mechanisms in rock blasting, such as tensile reflected waves, compressive stress wave, gas pressure, flexural rupture and nuclide stress flow, were propounded by various researchers in literatures but there is no agreement between the various researchers. However, it is generally accepted that the tensile reflected wave, and the coalescence of stress wave and gas pressure are the mainly reasons contributing to the rock fracture and fragmentation in rock blasting. Correspondingly, they are reviewed here in detail. At one stage, the tensile reflected waves were considered as the predominant means of rock fracture and fragmentation in rock blasting, which is true especially for the situations when there are one or more free surfaces. Hino [17] proposed that the major rock failure was caused by the reflection in tension of the primary shock wave. According to his theory, as depicted in Fig. 1(a–c) and (e), while the detonation of an explosive charge produces a crushed zone, a shockwave with high peak pressure but short of duration propagates outwards as a compressive wave. As the shock wave attenuates, it does not produce more breakage besides the cracks in the crushed zone. After it reaches free surface, the compressive wave reflects as a tensile wave. Since the tensile strength of rock is smaller than its compressive strength, the rock fracture occurs at the areas with the intensive tensile wave, i.e. the free surface area and its vicinities. If the compressive wave remains after the first shock wave reaches free surfaces and is reflected as tensile wave, the processes are repeated at newly produced free surface [17]. Theoretically, the processes will also occur at joints which naturally exist in the rock. As well known, in mining and civil engineer-

323

ing blasting practices, a free surface is normally made before blasting and the fragmentation is produced mainly from the free surface areas. Thus, this model seems reasonable. Nevertheless, if there is no any free surface, the tensile reflected wave may not play such a significant role. Bhandari [18] conducted the single blasthole tests in both cement-mortar and granite and concluded that in large fragments (e.g. boulders) the stress wave reflection and scabbing actions were weak. The gas pressure was once considered playing a significant role in rock blast [19]. Dally et al. [20] focused their study on the effect of the gas pressure from the combustion products of an explosive charge on the rock fracture process. Two plane models with a centrally located circular charge were constructed by them. In one of these two models, the charge was vented, so that the cracks produced around borehole could be caused only by the stress wave. In the other model, the charge was contained with a special sealing device in order to observe the action of the gas pressure in extending the cracks. They concluded that the containment of the charge (i.e. the containment of gas) in the hole caused more extensive fracturing. Besides Hino [17] and Dally et al. [20], there were other researchers concentrating on just limited aspects of the blasting process, such as either the initial explosive strain pulse [21–24] or the gas pressure [25–28], which were thought as the main cause of the rock fracture and fragmentation, while other factors were neglected [29]. However, more sophisticated theory should consider all the aspects controlling the rock fracture and fragmentation process during the rock blasting. The coalescence of the compressive stress wave and the gas pressure are widely accepted to result in the rock fracture and fragmentation in rock blasting in the literature [30]. Kisslinger [6] divided the region around the detonation into three zones: a strong shock zone (which was mainly produced by the shock wave), a transitional and nonlinear zone (in which both the shock wave and the gas pressure played a significant role), and the elastic region (which was produced under the gas pressure). Generally speaking, the three zones are produced by the forces exerted by the gas pressure and the stress wave simultaneously and it is almost impossible to separate the two principal blast forces. However, Kutter and Fairhurst [29] studied the roles of the stress wave and the gas pressure respectively in producing the rock fragmentation by underground blasts. In their studies, a pulse generated by an underwater spark discharge was used to simulate the explosive wave, and the pressurized oil was used to simulate expanding combustion i.e. the gas pressure. Since this theory has been widely accepted, we introduced it here in detail in order to provide a comparison with the numerical results to be presented later in this paper. Fig. 1(a)–(d) depicted the blast-induced rock fracture process proposed by Kutter and Fairhurst [29]. As shown in Fig. 1(a), the strong-shock zone lies in the region immediately around the borehole. When the chemical charge is detonated, a high temperature and density gas coupled with an extremely high pressure pulse, i.e. explosive wave, is generated. The high pressure pulse transmits in the rock adjacent to the borehole producing a dilatational wave, which propagates away from the charge at the sonic velocity in the rock [20]. The corresponding high pressure, which may exceed several-fold of the compressive strength of the rock, is exerted on the rock immediately around the borehole resulting in the vicinity around the borehole being intensively crushed and shattered. The first zone, i.e. the crush zone, is formed as shown in Fig. 1(a). In the crushed zone, the elastic rigidity of the rock is completely insignificant [29] and the crush is caused by both the compressive stress and the tangential stress but the compressive stress plays a more significant role. Moreover, the crushed zone is characterized by the shattered, smallest and relatively uniform particle size. The second zone, i.e. the transitional and non-linear zone or the cracked zone shown in Fig. 1(b), is characterized by

324

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

Transitional, non-linear zone

Fragment formation zone

Crushed zone

a) Formation of crushed zone

b) Formation of cracked zone

c) Long radial crack initiation and propagation Reflected tensile wave

Crack or Joint

Shock wave d) Fracture pattern of a fully contained explosion

Crushed zone

e) Stress wave propagation and reflection at the free surface or joint

Fig. 1. Rock fracture process by blasts in an infinite rock mass (a–d) (after Kutter and Fairhurst [29]) or a rock mass with one free surface (a–c, and e) (after Hino [17]).

the fracture phenomena from severe crushing through plastic deformation to partial fracturing and the rapidly increasing particle size as the radial distance increases [29]. As the explosive wave attenuates, the proportion of the fracture initiated by the compressive stress decreases, and, instead, that initiated by the tensile stress increases. In this zone, the fractures propagate predominately along the radial direction. The third zone, i.e. the fragment formation zone according to Saharan et al. [16] or the elastic zone according to Kutter and Fairhurst [29], appears as the explosive wave continues to attenuate, where the compressive stress rarely meets the critical strength of the rock. Therefore, there are no more fractures caused by the compressive stress. However, the tangential stress is still large enough to cause radial fractures because the tensile strength of the rock is considerably smaller than its compressive strength. As shown in Fig. 1(c), the fractures mainly propagate along the radial direction. Besides the explosive wave, a high-pressure gas is generated simultaneously when the chemical charge is detonated. The high pressure gas promotes and accelerates the propagation of the fractures. As shown in Fig. 1(d), the crushed zone and the non-liner zone are expanded due to the penetration of the high pressure gas into the fractures in these zones. At the same time, the gas pressure gradually decreases during its penetration process. It should be noted that, although the explosive wave and high-pressure gas are produced simultaneously and they are exerted on the rock almost at the same time, the expansion of the high pressure gas lags behind the propagation of the explosive wave. Controversies exist for the respective roles of the gas pressure and the shock wave energy in dynamic fracture initiation and propagation [15,16]. It is uncertain what respective amount of

the energy is liberated and imparted in different situations of engineering blasting operations in different rock types. Moreover, the explanations put forth so far also ignore the strain-rate dependant material response. The controversies probably stem from the observation difficulties associated with prevalent experimental techniques, which are due to the extremely short duration of the explosive shock wave load on rocks (in the order of microseconds), the very fast rock fracturing process (the cracking speed may be up to one third of the wave propagation speed in rock) occurring under the explosive high-pressure gaseous products, and the rock debris casting apart from the heterogeneity of rocks. The current experimental techniques do not have full control on the blasting experiments and the development of other experimental tools is much needed in this regard. Moreover, the coupling of the dynamic and quasi-static loading mechanisms and the inherently dynamic nature of the fracture process inhibits a robust analytical treatment of the damage process induced in rock by blasting. With the rapid developments in computing power, interactive computer graphics and topological data structure, numerical models derived from sound mechanical principle have provided a promising approach to untangle the blasting problem and to predict the rock fragmentation process induced by blasting under certain conditions. 2.2. State-of-the-art review of numerical modelling of rock fracture process by rock blasting Since Preece et al. [31] and Latham et al. [15] presented excellent reviews on the computer simulation of rock blasting conducted till 1990s, only the numerical modellings completed

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

thereafter are reviewed here (Table 1). As can be seen from Table 1, the numerical modelling of rock blasting available in literatures can be classified according to the types of numerical method, software code, and material constitutive model used in the numerical simulation, as well as the types of rock/rockmass, fracture and fragmentation modelled in the numerical study, which are briefly explained below. As can be seen from Table 1, the continuous and discontinuous methods are mainly implemented to investigate the fundamental processes taking place in blasting although the continuous method varies from the finite element software such as ABAQUS, LS-DYNA, AUTODYN and RFPA-RT to the finite different software such as FLAC and the discontinuous method varies from the commercial software such as PFC2D/3D to the free and opensource code such as DDA. The continuous method usually focuses on modelling the blast-induced damage while the discontinuous method concentrates on displacement of fragments resultant from detonation-induced crack propagation. Thorne et al. [32] described a finite element model based on the constitutive damage concept developed by Grady and Kipp [33], Taylor et al. [34], and Kuszmaul [35], where the fragmentation is achieved by the continual accumulation of damage in an isotropic material that includes microcracks. An important aspect of their simulation was the recognition of the gradual degradation of the material strength with time and the incorporation of microcracks which could interact and grow into macrocracks. Fourney et al. [36] implemented a finite element code to determine the displacement of several points in a crater region and compared these displacements with those measured in 3D laboratory tests, which were conducted in the models made from a very high-strength and low-porosity cement. Cho and Kaneko [37] proposed a finite element method to verify the dynamic fracture mechanism related to the blast-induced borehole breakdown. The stress wave propagation was modelled while the flow of the explosive gas was ignored. Their results reproduced the formation of the crushed zone and the long radial crack propagation. Liu et al. [38] investigated the crack propagation induced by the blast in the TASQ tunnel of Äspö Hard Rock Laboratory using the rock and tool interaction code (RFPA-RT) on the bases of the finite element method coupled with a gas expansion model. It is worth noting that the gas expansion model considers the gas immigration during the crack propagation process, where the explosive gas entered the formed cracks and further drove the cracks to propagate forming long radial cracks. Rouabhi et al. [39] developed a phenomenological constitutive model and implemented it into a finite element code to model the dynamic behaviour and fragmentation in rock blasting although the fragmentation was modelled analytically and there was no evidence of cracks or fragments shown in results. Liu et al. [40] used LS-DYNA to model the stress wave propagation in mining production blasting with the detonators in various locations of the boreholes in mining sublevel caving and compared the overstressed regions with the observed phenomena in drifts after the mining production blasting in LKAB’s Malmberget mine. Later, a number of researchers [41–45] implemented LS-DYNA with more complicated material models such as Taylor-Chen-Kuszmaul (TCK) continuum damage model, Johnson-Holmquist (JH) material model and Riedel-Thoma-Hiermaier (RHT) material model to simulate the stress wave propagation, the formation of the crushed zone and the crack zone, and the initiation and propagation of the long radial cracks induced by rock blasting and their behaviours under different loading and boundary conditions. The fracture process was modelled using the so-called erosion algorithm, which had the capability of treating the excessive element distortion problem. The element was immediately deleted when the material response in an element reached a certain erosion criterion. For example, in the study conducted by Wang et al. [41], the critical tensile damage was used as the erosion criterion. When the tensile damage

325

reached the critical value of 0.6, the associated element was immediately deleted. The deletion process was irreversible, which meant that, when the applied load was reversed, the deleted material would not be able to offer further resistance. This technique could be employed to capture the physical fracture process if no significant reverse loading occurred. These research projects have made LS-DYNA, or the explicit dynamic finite element method, one of the most popular and powerful methods for the numerical simulation of rock blasting due to its small time step and explicit time integration scheme. The small time steps ensure numerical accuracy and make the explicit method be especially suitable for rock blasting modelling since the duration of rock blasting is very short, i.e. in the order of milliseconds. The explicit time integration scheme does not require a factorization of the stiffness matrix in the step-by-step solution and the solution can essentially be carried out on the element level and relatively little high-speed storage is required, which greatly improves the calculation efficiency and saves the computation time. The other explicit dynamic finite element methods widely used in modelling rock blasting are ABAQUS [46,47] and AUTODYN [48,49]. For example, Saharan and Mitri [46] implemented ABAQUS with a brittle, Rankine failure-type material model and the element elimination technique to simulate the initiation and growth of fractures in the rock under the effect of blast-induced dynamic pressure pulse. Liu et al. [47] employed ABAQUS with a continuum damage model and the element erosion algorithm to model a series of rock blasts including the blasting in a single borehole in an unbounded rock mass, that in a rock mass with one free surface, the simultaneous blasting and the consecutive blasting with various delay times. AUTODYN was implemented by Zhu et al. [48] to simulate the blasting-induced crack initiation and propagation in rocks with the focus on the shock or stress wave induced cracks in the rock. In their model, coupling materials were placed in the borehole between the explosive charge and the borehole wall to simulate the effect of the subsequent high-pressure explosive gases on further development of the cracks induced by the stress waves. Liang et al. [49] implemented AUTODYN to investigate the effects of decoupled charge blasting on rock fragmentation efficiency. On the basis of the comparison between the numerical modelling and laboratory experiments, it was concluded that the gap between the explosive and the borehole wall in the decoupled charge blasting decreased the peak of the shock wave, extended the action time, and weakened the loading strain rate. An application of the discrete element modelling to rock blasting includes the work by Preece [50] who developed a discrete element method to simulate the motion of circular particles caused by the explosive gases resulting from a blast. However, that study did not model the rock fragmentation induced by the explosive, since the rock was assumed to have already been fragmented at the start of the calculations. Unlike the work of Preece [50], Song and Kim [51] developed a discrete element model based on Newton’s second law of motion to study the rock fracture process induced by blasting. In their model, the rock was represented using point masses or particles which were interconnected through massless springs and the inhomogeneities, which were the characteristics of geological materials such as rock, were incorporated. Their modelling results indicated that the discrete model was capable of simulating various dynamic fracture phenomena including crack branching and tortuous propagation associated with the rock blasting. Similar to Song and Kim [51], Potyondy et al. [52] implemented a parallel bond model into the commercial discrete element code PFC2D to model the detonation-induced crack initiation and propagation. In their parallel bond model, a material strength capability was incorporated by adding links or bonds between the discrete elements. These links were used to hold the particles together, maintain the geometrical shape of the particles,

326

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

Table 1 Literature review on numerical modelling of rock fragmentation process by rock blast. Continuous (C), discontinuous (D) or hybrid

Numerical method

Numerical code

Constitutive model

Type of rocks modelled

Modelled results including stress wave, gas flow, crack initiation and propagation, formation of crushed and cracked zones, fragment separation, casting, and muck-piling

References

C

FEM(2D)

Elastic

Concrete, homogeneous

FEM(2D)

FPZ, fracture mechanics

Granite, inhomogeneity, statistical distribution

C

FEM(2D)

RFPA-RT

Elastic-brittle

Diorite, heterogeneous

C

FEM(2D) FEM(2D)

Quasi-brittle material model Elastic damage model

Limestone, homogeneous

C

In-house code LS-DYNA

C

FEM(2D)

LS-DYNA

Oil shale, homogeneous

C

FEM(2D)

AUTODYN

C

FEM(2D)

ABAQUS

C

FEM(2D)

LS-DYNA

TCK continuum damage model with an erosion algorithm Modified principal stress failure criterion Brittle, Rankine failure-type material model JH material model

C

FEM(3D)

ANSYSLSDYNA

Damage model

Granite, homogeneous

C

FEM(2D)

ABAQUS

Damage model

Generic, homogeneous

C

FEM(2D)

AUTODYN

Concrete, homogeneous

C

FEM(3D)

LS-DYNA

Modified principal stress failure criterion RHT material model

C

FDM(3D)

FLAC3D

Mohr-Coulomb failure criterion

Granite and sandstone, homogeneous

Determine displacements of several points in the crater region. No cracks or fragments are modelled Crushed zone and long radial crack are modelled. No fragment or fragmentation separation, No gas pressurization phase Radial crack propagation and the crack zone is modelled. No fragment separation. Stress wave is ignored and only gas pressurization phase is modelled Fragmentation has been modelled analytically. There is no evidence of cracks or fragments Stress wave propagation is modelled and fragment is predicted. However, there is no crack propagation and fragment separation Stress wave propagation is modelled while the gas pressurization phase is ignored Tensile failure and shear failure ate modelled. No fragment or fragment separation Crack formed by smeared elements. No fragment or fragment separation Crushed zone and long radial crack are modelled. No fragment or fragmentation separation, No gas pressurization phase Damage zone is modelled using peak particle velocity (PPV) damage criterion and plastic strain criterion. No fragment or fragmentation separation, No gas pressurization phase Crack formed by smeared elements. No fragment or fragment separation Tensile failure and shear failure ate modelled. No fragment or fragment separation Crushed zone and long radial crack are modelled. No fragment or fragmentation separation, No gas pressurization phase Crushed zone is modelled. No crack initiation or propagation. No fragment separation

[36]

C

In-house code In-house code

D

DEM(2D)

Elastic

Generic, homogeneous

DEM(2D)

Lattice network model: elastic and brittle

Charcoal gray granite, heterogeneous

Gas flow is modelled. No fracture or fragmentation are modelled Stress wave propagation, quasistatic pressure acting against the walls of the borehole and induced fractures. Long radial crack propagation and crushed zone is modelled. No fragment separation Radial crack is formed. No crushed zone, no fragments. shock wave and fluid coupling Crushed zone formation and long radial crack propagation are modelled. No fragment separation Ignore the stress wave propagation and only consider the gas pressurization process. No crack propagation. Fragment throw and muckpile is modelled Stress wave propagation across fractured rock masses No formation of the crushed and cracked zone and long radial crack initiation, fragment cast is modelled No gas pressurization phase and throw were modelled

[50]

D

In-house code In-house code

Magnetite, homogeneous

Diorite, homogeneous Generic, homogeneous Granite, homogeneous

Magnetite

D

DEM(3D)

PFC3D

Parallel bond model

Generic, homogeneous

D

DEM(2D) DEM(2D)

Cohesive elastic-brittle model Elastic

Generic, homogeneous

D

In-house code DDA_BLAST

D

DEM(2D)

UDEC

Elastic

Generic, jointed rock mass

D

DEM(2D)

DDA

Artificial joint model

Granite, concrete, homogeneous

D

DEM(2D)

In-house code

Mohr-Coulomb failure criterion

Generic, homogeneous

Hybrid

FEM-DEM (2D)

Y2D

Strain-softening-based model

Generic, homogeneous

Hybrid

FDM-DEM (2D)

FLAC-PFC

Granite, homogeneous

Hybrid

FDM-DEM (3D)

HSBM

Mohr-Coulomb constitutive model for FLAC and parallel bond model for PFC Elastic plastic and lattice model

Hybrid

FEM-DEM (2D)

Y-2D/3D IDE

Combined single and smear fracture model

Generic, heterogeneity based on Weibull distribution

Hybrid

DEM-SPH (2D)

In-house code

Bonded particle model

Generic, homogeneous

Generic, homogeneous

Synthetic rock, homogeneous, isotropic

[37]

[38]

[39] [40]

[41,42] [48] [46] [43]

[44]

[47] [49] [45]

[63]

[51]

[52] [53] [55]

[56] [58,59]

[64]

Gas penetration and opening of fracture. No fragment separation, no effect of geological discontinuities Blast-induced radial cracks are modelled

[54]

Both stress wave and gas flow are modelled. No clear crack formed. Fragment separated, but mostly particles dispersed in space Both stress wave and gas pressurization are modelled. No clear crack initiation and propagation, no formation of the crushed zone and the cracked zone. Fragment separation and muckpiling Crack initiation and propagation, crushed zone and cracked zone are modelled. Fragment is separated but no casted up. Solid-fluid coupling is modelled. No stress wave propagation

[60]

[57]

[62]

[61]

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

and add the strength to the various geological formations. These links could break due to the tension, shear or compression stress satisfying the strength of the material. A discrete element method was developed by Donze et al. [53] to examine the fracturing process during the stress wave loading in a cohesive elastic-brittle material, which took into account the initiation and propagation of narrow discontinuities around the explosion cavity during a blasting event. Munjiza et al. [54] described such an effort of modelling the blast-induced fragmentation using a combination of finite element and discrete element methods with a strainsoftening-based fracture model, in which each discrete element was broken into two separate discrete elements whenever the stress satisfied the material strength. Munjiza et al.’s [54] model could handle the gas penetration and the opening of fractures due to the gas pressurization induced by rock blasting although it couldn’t consider the effect of the geological discontinuities on the rock fracture process. Mortazavi and Katsabanis [55] employed DDA_BLAST to simulate typical blasting problems in jointed rock mass and delve into the macro scale mechanisms involved in the gas pressurization phase of the blasting process. On the basis of the numerical results, they discussed the effects of the discontinuities on the burden breakage process. However, their DDA model only considered the gas pressurization phase of the blasting process and assumed that the rock was already fragmented in-situ due to the intersection of pre-existing discontinuities and the passage of the stress wave. Zhao et al. [56] studied the stress wave attenuation across jointed rock masses using UDEC but no damage or fracture was modelled in their study. Saiang [57] implemented a coupled continuum-discontinuum model using FLAC and PFC2D to simulate the blast-induced damaged zone around excavations, where the inner segment of the model was simulated using PFC2D while the outer segment was simulated using FLAC. However, in Saiang’s [57] study, the artificial blast-induced radial cracks were individually implemented into the continuum-discontinuum model instead of being induced by blasting. Ning et al. [58,59] implemented DDA to simulate the rock mass failures by the blast-induced high pressure expansion, in which the instant explosion gas pressure is loaded on the main blast chamber walls and the artificially-placed surrounding connected fracture surfaces. Their DDA simulation modelled the rock mass casting process and the formation of the final blasting pile but ignored the stress wave and then the formation of the crushed and cracked zone and the long crack initiation and propagation. More recently, Onederra et al. [60] developed a hybrid stress blasting model (HSBM) to model the blasting process from the explosive detonation to the rock breakage and fragmentation. The HSBM code used a combination of the continuum numerical technique and the discrete element method to model the explosive and rock interaction. The near-borehole area was represented using FLAC while DEM was used to represent the rock mass. The DEM representation used a lattice-scheme where the rock matrix was created as a collection of randomly distributed points connected by spring. However, although the HSBM code had the ability of combining the effect of the stress wave and the gas flow into the calculations and modelling the fragment separation, there were not the clear crack initiation and propagation, stress wave propagation, formation of the crushed and cracked zones, and fragments except separated particles mostly dispersed in space. Fakhimi and Lanari [61] proposed a hybrid DEM-SPH (smooth particle hydrodynamics) method to simulate the rock blasting in which DEM with a bonded particle model was utilized to mimic the behaviour of rock and SPH was implemented to simulate gas flow. Their model successfully simulated the crack propagation in rock, in particular, the crushed zone, radial crack and surface spalling. However, no fragment casting and muck-piling were modelled in their study.

327

In sum, the finite element method, discrete element method and their hybrid are the main methods used in the literatures to simulate the blast-induced rock fracture and fragmentation although other methods such as the finite difference method [57,60] and SPH [61] are also reported in the literatures. The finite element method is good at modelling the blast-induced damage and crack propagation. Most of the finite element modellings presented above have successfully modelled the formation of the crushed zone and the cracked zone, and the radial crack initiation and propagation although they may focus on either the stress wave propagation or the gas pressurization instead of both. However, none of the finite element modellings is able to simulate the fragment formation, separation, casting, and muck-piling processes due to its inherent continuous hypothesis. Compared with the studies on the finite element modellings, most of the discrete element modellings reviewed above successfully model the fragment formation and separation. However, most of them have not well simulated the crack initiation and propagation, the shape of the fragments except separated particles mostly dispersed in space, and the fragment casting and muck-piling processes. Correspondingly, the hybrid finite-discrete element method is a promising approach to model the rock blasting as it takes advantage of the strengths of each method. The literature review presented above has already indicated the emerging hybrid finite-discrete element method is indeed the most suitable method to model both the stress wave propagation and the explosive gas pressurization in spite of its shortcomings reviewed above. However, until this moment, very few studies have attempted to model the complete blasting process from the detonation of the explosive to the formation of the muckpiling. This paper aims to model the complete rock blasting process using a hybrid finite-discrete element method implemented by Liu et al. [62] in order to understand the rock blasting mechanisms and then improve the rock blasting efficiency. 3. Hybrid finite-discrete element method for modelling rock blasting The hybrid finite-discrete element method to be used for the blast modelling in the following sections is implemented by Liu et al. [62] on the basis of their previous enriched finite element codes RFPA-RT2D [65,66] and TunGeo3D [67], and the opensource combined finite-discrete element libraries Y2D and Y3D originally developed by Munjiza [68] and Xiang et al. [69], respectively. Similar implementation of the hybrid finite-discrete element method has been conducted by Tatone and Grasselli [70] and Mahabadi et al. [71] to model the rock failure processes using the Brazilian tensile strength tests and uniaxial compressive strength tests in 2D and 3D, respectively. The hybrid finitediscrete element method consists of the following components: contact detection, contact interaction between individual bodies, deformability and transition from continuum to discontinuum, temporal integration scheme, and computational fluid dynamics. The blast-induced rock fracture and fragmentation is a coupled problem involving the expansion of the detonation gas with high pressure and temperature, and the rock fracture and fragmentation. In other words, the following two aspects are essential for successful rock blasting modelling: (1) modelling the deformability, fracture and fragmentation of rock, and (2) modelling the gas pressure as surface loads on the rock around the borehole and its flow into subsequently generated fractures. Among the components of the hybrid finite-discrete element method to be used in the following sections, the transition from continuum to discontinuum through fracture and fragmentation is the key feature for the first aspect of the rock blasting modelling, which makes the hybrid

328

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

finite-discrete element method distinguish from the continuumbased method and the discontinuum-based method reviewed in Section 2.2 in terms of modelling the rock blasting. The computational fluid dynamics is the key component for the second aspect of the rock blasting modelling. Thus, the algorithms of the transition from continuum to discontinuum through fracture and fragmentation and the computational fluid dynamics implemented in the hybrid finite-discrete element method are introduced in the following sections. Compared with the previous studies conducted by Munjiza et al. [54,68], the transition from continuum to discontinuum in the hybrid finite-discrete element method occurs through three fracture modes, i.e. model I (i.e. tensile failure mode), model II (i.e. shear failure mode) and mixed-mode I–II while the transition only occurs through the mode I fracture in the original open-source combined finite-discrete element library Y2D. Correspondingly, the transition from continuum to discontinuum through fracture and fragmentation is introduced in Section 3.1 in detail by focusing on the transition through the mode II and mix-mode I–II fractures. As for the computational fluid dynamics, no improvements are made except an equation of state is used to obtain the relationship between the pressure and time in the hybrid finite-discrete element method compared with the previous studies [54,68]. Thus, Section 3.2 briefly introduces the gas expansion and flow by stressing the equation of state. Moreover, since the behaviour of rock under impact load is loading rate dependent and the original open-source combined finite-discrete element library Y2D has not taken it into account, an empirical relationship between the physical-mechanical properties of rock and the loading rate is implemented into the hybrid finitediscrete element method. Besides, a non-reflection boundary, which is unavailable in the original open-source combined finitediscrete element library Y2D, too, is implemented into the hybrid finite-discrete element method in order to model the infinite rock mass involved in the rock blasting modelling. 3.1. Transition from continuum to discontinuum through fracture and fragmentation The hybrid finite-discrete element method regards the rock blasting problem to be modelled as a number of interactive discrete elements with general shape and size. The interaction of the discrete elements is governed by the so-called contact law: a stiffness (i.e. spring) in the normal direction and a stiffness and friction angle (i.e. spring-slip) in the tangential directions. The interaction forces developed at contact points are determined using linear functions of the deformations of the spring and/or spring-slip surfaces and resolved into normal and tangential components. Each discrete element is of a general shape and size and is then discretised into a series of finite elements to analyse its deformability. An explicit and large strain Lagrangian formulation for the 10-node tetrahedral elements in 3D or constant strain elements in 2D is used to represent the element deformations. Based on the Gauss’ theorem of converting volume integrals into surface integrals in 3D or area integrals into line integrals in 2D, the increments of element strain can be obtained for each time step. The stress increments can then be obtained by invoking the constitutive equations of the modelled materials. In the hybrid finite-discrete element method, the discrete elements may fracture and fragment depending on the calculated stress and strain of the discretised finite elements in the discrete elements. The fracture and fragmentation of the discrete elements result in the transition from continuum to discontinuum and irregular-shaped and further breakable fragments although the shape and size of the fragments depend on the mesh discretisation. The finite elements are bonded together using a four-node joint element before fracturing and the cracks are assumed to coincide

with the surfaces of the joint element during fracturing. Separation of these surfaces induces a bonding stress, which is taken to be a function of the size of separation. At any point on the surfaces of a crack, the separation d can be divided into two components in Eq. (1)

d ¼ dn n þ ds t

ð1Þ

where n and t are the unit vectors in the normal and tangential directions, respectively, of the surface at such a point, dn and ds are the magnitudes of the components of d in the normal and tangential directions, respectively. Depending on the relative separations in the normal and tangential directions and the local stress, the joint element may damage and fracture in the mode I (i.e. tensile failure mode), the mode II (i.e. shear failure mode) or the mixed mode I– II. As shown in Fig. 2(a), if the opening of the joint element, i.e. the separation in the normal direction dn , reaches a critical displacement, dnp , prescribed by the tensile strength of the element, rt , the tensile failure (i.e. mode I) occurs. As the opening of the joint element is beyond dnp , i.e. dn > dnp , the bonding stress in the normal direction, i.e. rn , is assumed to gradually decrease till an ultimate stress according to a mechanical damage model. Finally, when the opening of the joint element exceeds the ultimate opening displacement, dnu the joint element is completely broken and the adjacent two finite elements are completely separated to become two discrete elements. The ultimate opening displacement dnu is determined by the mode I fracture energy release rate Gf I , which includes the contributions of the surface energy 2c of the newly formed crack surface along with the energy consumed in the damage process. The mode I fracture energy release rate Gf I is equal to the area under the curve of the bonding stress and the opening displacement depicted in Fig. 2(a), i.e.

Z Gf I ¼

dnu

dnp

rn ðdn Þddn

ð2Þ

In summary, if the mode I fracture initiation and propagation occur, the bonding stress in the normal direction can be calculated using Eq. (3)

rn ¼

8  2  > dn dn > 2   rt  > < dnp dnp > > > :

if

0 6 dn 6 dnp

f ðDÞ  rt

if

dnp 6 dn 6 dnu

0

if

dn P dnu

ð3Þ

where D is a damage variable between 0 and 1, f ðDÞ is a damage function described in the mechanical damage model [65,66], and all other parameters have the same meanings as those introduced above. It should be noted that, in the combined single and smeared crack model originally implemented into the combined finitediscrete element method by Munjiza et al. [72], f ðDÞ is a heuristic function that approximates the shape of the experimental postpeak stress-displacement curves measured in direct tension tests of concrete [73]. If the separation of the joint element in the tangential direction, i.e. the sliding displacement ds , reaches a critical slip, dsp , defined by the shear strength of the joint element, rc , the shear failure (i.e. mode II) occurs. As the sliding displacement of the joint element is beyond dsp , i.e. ds P dsp , the bonding stress in the tangential direction, i.e. s, is assumed to gradually decrease till a residual stress according to a mechanical damage model. Finally, when the sliding displacement of the joint element exceeds the residual opening displacement dsr , a physical separation is formed along the tangential direction and the bonding stress in the tangential direction becomes a purely frictional resistance defined by the Coloumb’s model. The residual sliding displacement dsr is prescribed by the mode II fracture energy release rate Gf II , which includes the contributions of the surface energy 2c of the newly formed crack surface

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

329

where all the parameters have the same meanings as those introduced above. When the joint element is yielding and damaging, the bonding stresses rn and s are applied to the nodes of the adjacent finite elements as the equivalent nodal forces. When the ultimate opening or the residual sliding are reached, the joint element is effectively deleted from the model and the initially bonded finite elements separate to become the discrete elements. Fig. 3 illustrates the flowchart of the fracture and fragmentation algorithm introduced above. After fracture and fragmentation, an explicit central difference scheme is applied in the hybrid finitediscrete element method to integrate the equations of motion of either the initially discrete elements or the discrete elements formed by the fracture and fragmentation algorithm. The explicit integration scheme is much faster than the implicit one and requires few high-speed storages since no factorization of the stiffness matrix is involved although it requires the time step must be smaller than a critical value related to material properties and meshes, which makes the time step becomes very small when the mesh size is small. The unknown variables, i.e. contact forces on the discrete elements’ boundary or stresses in the internal elements, are determined locally at each time step from the known variables on the boundaries and in the elements and their immediate neighbours.

a) Under tension conditions

3.2. Detonation-induced gas expansion and flow through fracturing rock

b) Under shear conditions Fig. 2. Relationship between the bonding stress and opening/sliding displacement under tension and shear conditions.

along with the energy consumed not only in the damage process but also by the friction. The mode II fracture energy release rate Gf II is equal to the area under the curve of the bonding stress and the sliding displacement depicted in Fig. 2(b), i.e.

Z Gf II ¼

dsr

dsp

½sðds Þ  sr dds

ð4Þ

In summary, if the mode II fracture initiation and propagation occur, the bonding stress in the tangential direction can be calculated using Eq. (5)



8  2  > ds ds > 2   rc  > dsp dsp < > > > :

if

0 6 ds 6 dsp

gðDÞ  rc

if

dsp 6 ds 6 dsr

rn  tanð;f Þ

if

ds P dsr

ð5Þ

where gðDÞ is a damage function described in the mechanical damage model [65,66], ;f is the joint friction angle, and all other parameters have the same meanings as those introduced above. Under the blasting loads, the joint element may be subjected to a combination of the opening and sliding displacements instead of either pure mode I (i.e. opening) or pure mode II (i.e. sliding) separation. In this case, the resultant joint displacement may be large even if the opening or sliding displacements have not reached the ultimate opening or the residual sliding, respectively. Correspondingly, the mixed-mode I–II yield and damage occur if Eq. (6) is satisfied



dn  dnp dnu  dnp

2

 þ

ds  dsp dsr  dsp

For the rock blasting analysis, the modelling of the interaction behaviour between the detonation products and the surrounding rock mass is a major requirement. In most of existing methods for rock blasting modelling, the pressure-time histories generated from empirical equations are used to model the detonationinduced pressure [37,46,53]. The approach has many limitations and involves in various crude approximations. This study is not intended to model the complicated chemical reaction and the associated early-time rock crushing behaviour in blast using the hybrid finite-discrete element method. Instead, the chemical reaction in blast is modelled using the commercial explicit finite element code AUTODYN. A relationship between the blast-induced instantaneous gas pressure and the time is obtained through the AUTODYN modelling and then implemented into the hybrid finite-discrete element method. As the detonation gas expands, the rock mass around the borehole fails to form the crushed zone and cracks, and then the detonation gas flows through the crushed zone and the cracks propagating out of the crushed zone. Meanwhile the detonation gas exerts pressures on the cracks and the crushed zone to result in further fracture and fragmentation. In order to model this process, a relatively simple model for gas flow through the fractured solid is implemented into the hybrid finite-discrete element method, which is actually based on the detonation wave speed and the distance of the observation point from the detonation point. The explosive gas is presented only within the gas zone, which is best defined as a circle around the borehole. The spatial and temporal distribution of the gas pressure in the circle are determined on the basis of the distance of the circle from the detonation point and the instant time following a iteration procedure described by Munjiza et al. [54,68]. From the procedure briefly introduced above, the spatial and temporal pressure distribution of the detonation gas is finally acquired. 3.3. Effect of loading rate on dynamic behaviour of rock

2 P1

ð6Þ

The effect of the loading rate on the deformation and fracture behaviour of rock is often ignored, which is usually acceptable

330

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

Calculate interaction forces between discrete elements (bodies)

Calculate strain and stress for finite elements inside discrete bodies

Compute separations and of joint elements between finite elements

No

No

Eq. 6 is satisfied

No

Yes: Mix mode damage

Tensile damage Yes

Yes Shear damage

No

No

Yes Shear fracture

Tensile fracture Yes Compute normal bonding stress according to Eq. 3

Compute shear bonding stress according to Eq. 5

Update forces applied on original and newly formed discrete bodies

Calculate acceleration and angular rotation of discrete bodies (elements)

Fig. 3. Flowchart of fracture algorithm implemented into the hybrid finite-discrete element method.

since most rock engineering applications can be regarded as a static or quasi-static problem and the loading rate won’t vary in a wide range. However, for the rock engineering applications involving in strong dynamic loadings, such as the rock blasting considered in this study, it is essential to take into account the effect of loading rate on the dynamic deformation and fracture behaviour of rock. Previous studies have indicated that the strength [74] and fracture toughness [75] of rock material generally increase with the loading rate increasing and several mechanisms have proposed to explain the change of the strength and fracture toughness due to the effect of the loading rate. A semi-log formula in Eq. (7) was proposed by Zhao [74] to describe the relation between the uniaxial compressive strength and the loading rate:



rcd ¼ A  log



r_ cd þ rc r_ c

ð7Þ

where rcd is the dynamic uniaxial compressive strength (MPa), r_ cd is the dynamic loading rate (MPa/s), r_ c is the quasi-static loading

rate (approximately 5  102 MPa/s), rc is the uniaxial compressive strength at the quasi-static loading rate (MPa) and A is a material parameter depending on rock type. Zhao et al. [75] proposed a formula similar to Eq. (7) to fit the static and dynamic fracture toughness data obtained from their experiments. In the hybrid finitediscrete element method, the fracture toughness has not been directly used but the fracture energy release rate is instead implemented to regulate the damage process and the ultimate opening displacement or the residual sliding displacement. For the plane strain simplification of the rock blast modelling considered in this study, Gf I ¼ K 2I ð1  v 2 Þ=E, where Gf I is the fracture energy release rate, K I is the fracture toughness, v is the Poisson’s ratio and E is the elastic modulus. Correspondingly, the dynamic fracture energy release rate is assumed to increase with the loading rate increasing

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

according to Eq. (7). Therefore, in order to consider the effect of the loading rate on the dynamic fracture behaviour of rock in rock blasting, Eq. (7) is implemented into the hybrid finite-discrete element method to model the increases of the tensile and shear strengths and the mode I and mode II fracture energy release rates, which has been calibrated against the dynamic Brazilian tensile strength experiments. 3.4. Non-reflection boundary for modelling infinite rock In the hybrid finite-discrete element method, all boundaries are naturally reflection boundaries. However, the rock blasting modelling may involve in infinite rock and the need for nonreflection boundary naturally arises. Of course, the simplest way is to extend the model so that the boundaries are far enough and the modelling is not affected by the boundaries. However, due to the high speed of the stress waves induced by rock blasting, the extension of the model is not computationally practical. Lisjak and Grasselli [76] reviewed the different absorbing boundary conditions available in literatures and found the local solution proposed by Lysmer and Kuhlemeyer [77] was most suitable to be implemented into the hybrid finite-discrete element method. Following their approach, the viscous boundary tractions are used to numerically absorb kinetic energy of incident waves. For a wave with a particle velocity of VðV x ; V y Þ approaching a boundary, the surface tractions in the horizontal and vertical directions can be calculated using Eqs. (8) and (9), respectively:

rx ¼ a  q  V P 

@V x @t

ð8Þ

ry ¼ b  q  V S 

@V y @t

ð9Þ

where q is the material density, V P and V S are the velocities of longitudinal (P) and shear (S) waves, respectively, and a and b are constants, which are a ¼ b ¼ 1 for this study. The calculated forces are then applied to the boundaries as equivalent nodal forces to cancel the tractions at the boundaries caused by the blast-induced stress waves. Correspondingly, the non-reflection boundaries allows the necessary energy radiation, which are essential for modelling the blasting in infinite rock. 4. Hybrid modelling of rock fracture process and resultant fragment casting by rock blasting This section will use the hybrid finite-discrete element method introduced in Section 3 to model the stress wave propagation and the fracture initiation and propagation caused by various blasts, i.e.  a single-borehole blast in a rock mass with a free surface,  two-borehole simultaneous blasts, and  two-borehole consecutive blasts. The obtained results will then be compared with those documented in the literature. After that, a conceptual numerical model will be built for the mining production blast in mining sublevel caving to demonstrate the potential application of the hybrid finite-discrete element method in mining fragmentation by blasts. 4.1. Rock fracture process and resultant fragment casting by a single borehole blast One of the basic principles of designing rock blasting is the presence of a free face parallel or sub-parallel to the blast holes as detonation occurs. Thus a rock mass with one free surface and a borehole parallel to the free surface is first modelled using the

331

hybrid finite-discrete element method. As the denotation velocity of the explosive along the length of the borehole is very fast, it can be supposed that the explosive along the length of the borehole is detonated at the same time and the rock fracture mechanisms and the resultant fragmentation are similar in each section along the length of the borehole. Thus, the problem can be simplified as a two-dimensional plane strain problem although the simplification may limit the practical usefulness of the model. As a matter of fact, the two-dimensional model has been widely employed to study the basting mechanism under both experimental condition and numerical environment [29,36,46]. Furthermore, the results from the plane strain modelling can be directly compared with those documented in literatures to qualitatively calibrate the hybrid finite-discrete element method. As shown in Fig. 4(a), the numerical model comprises a borehole with a radius of 50 mm, which is 2 m far from the free surface, i.e. the top boundary. The other boundaries of the rock mass are put far away from the borehole and are treated as the nonreflection boundary according to Section 3.4 in order to eliminate any influence caused by the reflected waves from the boundary except the free surface. It should be noted that the size of the borehole in the numerical model is much bigger than that in the blasting practice since it corresponds to the size after the explosive chemical reaction in the borehole and the associated early-time rock crushing behaviour. As stated in Section 3.2, the complicated chemical reaction in the blast is first modelled using the commercial finite element code AUTODYN and the resultant instantaneous pressure is then applied in the borehole. Triangular elements are employed to discretize the model, and dense elements are employed in the interested area, i.e. the vicinity of the borehole, as shown in Fig. 4(b). Points A, B and C are marked in Fig. 4(a) for the convenience of presenting the simulation results. The static material properties of the rock mass and other input parameters are listed in Table 2. The dynamic strength and fracture energy of the rock mass are calculated according to Eq. (7) introduced in Section 3.3. Figs. 5 and 6 depict the temporal distributions of the blastinduced stress wave propagation and fracture process modelled using the hybrid finite-discrete element method. It should be noted that, following the sign convention in solid mechanics, the compressive stresses are taken as negative, while the tensile stresses are regarded as positive. Moreover, the colour in Fig. 5 represents the size of minor (compressive) principal stress and the magnitudes can be referred to the legend shown in Fig. 5(a). In Fig. 6, the blue colour represents the tensile failure while the red colour represents the compressive and mix-mode failures. 4.1.1. Stress and fracture initiation and propagation in the crushed zone While the explosive is detonated, the contained chemical charge in the borehole reacts rapidly, which results in the creation of high-pressure gases [29]. As the rapidly expanding highpressure gas interacts with the rock material surrounding the borehole, an intense stress wave is emitted into the rock and propagates radially out of the borehole, as shown in Fig. 5(a)–(c). At the time of 0.45 ms, the stress wave reaches the free surface, as shown in Fig. 5(d). During this period, in the immediate vicinity of the borehole, the wave energy is highly concentrated and produces extensive compressive breakage in that region, as shown by the cracks marked by the red colour in Fig. 6(a)–(d). This significantly modifies the borehole boundary on which the gas pressure subsequently acts. The high pressure, coupled with the stress wave emitted into the rock, cause the rock mass around the borehole to crush and shatter forming the so-called strongshock zone or crushed zone according to Kutter and Fairhurst [29] and Kisslinger [6].

332

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

Free surface Borehole: D = 100 mm 2m Point C Point B

0.55 m 1.65 m

Point A

3.35 m

12 m Non-reflection boundary

20 m

a) Geometrical model

b) Numerical model

Fig. 4. Geometrical and numerical models for rock blast in a rock mass with one free surface (the monitoring points A, B and C are placed in the elastic, cracked and crushed zones, respectively, to be induced by the rock blast).

In this study, the static tensile and compressive strengths of the rock mass are 20 MPa and 200 MPa, respectively, while the dynamic tensile and compressive strengths are determined automatically by the hybrid finite-discrete element method using Eq. (7). Hence the rock mass can more easily fail in the tensile mode than the compressive mode in both static and dynamic conditions. However, in the crushed zone, the rock mass is compressed into failure since the compressive stresses there are several-fold larger than the static compressive strength of the rock mass due to the blast-induced high gas pressure and stress wave propagation. As shown by the cracks marked using the red colour in Fig. 6(a)–(d), all of the blast-induced failures in the crushed zone till 0.45 ms are shear failures and none of them are tensile failures. Kutter and Fairhurst [29] and Kisslinger [6] conducted physical experiments of rock blasting and found that in the region immediately around the borehole, the explosive-induced stresses exceeded several-fold of the static compressive strength of the rock mass, and concluded that the rock mass in the crushed zone was intensely crushed and shattered by the compressive stress and the elastic rigidity of the rock mass in the crushed zone was completely insignificant. Thus, the hybrid finite-discrete element method has well modelled the formation of the crushed zone induced by rock blast. 4.1.2. Stress wave reflection and resultant tensile failure at the free surface As the compressive stress wave propagates radially out of the detonation borehole, the stress wave then reaches the free surface and is reflected back to form tensile stress wave propagating toward the borehole, which is quite clear from the colour changes near the free surface in Fig. 5(d)–(f), i.e. from the green1 colour (representing the compressive stress) to the red colour (representing the tensile stress). Correspondingly, the rock mass around the free surface starts to crack predominately in the tensile mode instead of the compressive mode when the tensile stress induced by the reflection of the compressive stress wave at the free surface is strong enough to satisfy the dynamic tensile strength of the rock mass. These tensile failures are marked using the blue colour in Fig. 6 (e)–(f). The modelled tensile failures around the free surface conform with the rock blast theory proposed by Hino [17], which concluded that the reflective tensile stress wave was the main reason of the rock fracture and fragmentation in rock blast, as depicted in Fig. 1 (e), since the tensile strength of the rock was much smaller than its compressive strength and the tensile failure was the most efficient way to fragment the rock in blast. 1 For interpretation of color in Fig. 5, the reader is referred to the web version of this article.

Table 2 Input parameters for the hybrid finite-discrete element model. Properties

Symbols

Values

Units

Young’s modulus Poisson’s ratio Density Tensile strength Compressive strength Internal friction angle Surface friction coefficient Model-I fracture energy release Mode-II fracture energy release

E

60 0.26 2600 20 200 30 0.1 50 250

GPa N/A Kg m3 MPa MPa  (degrees) N/A N m1 N m1

v q rt rc

; u GfI GfII

4.1.3. Long radial crack initiation and propagation With the exception of the stress wave reflections at the top free surface, the compressive stress wave in other three directions continues to propagate radially out of the borehole and results in the compressive failure causing the crushed zone to expand, as shown in Figs. 5 and 6(e)–(f). However, the compressive stress wave attenuates as it propagates. During this process, more fractures are initiated, which propagate from the edge of the crushed zone in the radial and circumferential direction, as shown in Fig. 6(g)– (i). The newly formed failure area is called the non-linear zone or the cracked zone [6,29]. The cracked zone is characterized by increasing importance of the shear resistance, although the compressive stress are still substantially above the dynamic compressive strength of the rock mass causing the compressive failures in some area of the cracked zone. However, as the stress wave attenuates, the proportion of the fractures induced by the tensile stresses is increasing while the proportion of the compressive failure is decreasing, as indicated by the increasing cracks marked by the blue colour (representing tensile failure) Fig. 6(g)–(i). As the stress wave continues to attenuate, the compressive stress at the leading edge of the stress wave can barely meet its critical value and no more failures can be induced by the compressive stress. Nevertheless, in the tailing part of the outgoing wave, there are tensile radial and tangential stresses and the tensile stresses are large enough to cause the tensile failures in the radial direction since the tensile strength of the rock is much smaller than its compressive strength. Correspondingly, the tensile cracks are further generated in the cracked zone and initiated at the edge of the cracked zone to form the fragment formation zone [16] or the elastic zone [29], as shown in Fig. 6(j)–(l). Meanwhile, gas plays an essential role from the very beginning of the explosive detonation. As introduced in Section 3.2, the detonation gas flow is implemented into the hybrid finite-discrete element method on the basis of the detonation wave speed and the distance of the observation point from the detonation point

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

333

2.19803e-5 Pa

a) t = 0.15 ms

-7.38696e+9 Pa

b) t = 0.25 ms

c) t = 0.35 ms

d) t = 0.45 ms

e) t = 0.55 ms

f) t = 0.75 ms

g) t = 1 ms

h) t = 1.5 ms

i) t = 2 ms

j) t = 2.5 ms

k) t = 3.5 ms

l) t = 4.5 ms

Fig. 5. Stress wave propagation induced by rock blast in a rock mass with one free surface.

according to an iteration procedure described by Munjiza et al. [54,68]. The detonation gas then penetrates into the cracks in all the three zones introduced above, i.e. the crushed zone, the cracked zone and the elastic zone, forces them to further propagate, and even drives some of them to propagate out of the cracked zone forming the long radial cracks. It should be noted that if the borehole is detonated in an infinite rock mass without free surface, more similar outline of the three zones and the long radial cracks to that caused by the fully contained explosion [29] can be obtained using the hybrid finite-discrete element method. However, since there is a free surface at the top boundary of the numerical model, the tensile cracks induced by the reflected stress wave, the long radial cracks driven by the stress wave and gas flow, and the shear cracks caused by the stress wave interact and coalesce with each other in the area near the free surface resulting in the rock fragments casting out there. As shown in Figs. 5 and 6(i)–(l), the hybrid finite-discrete element method can model well not only the transition from continuum to discontinuum of the rock mass under the coupled action of the stress wave propagation and the detonation gas flow in the rock blast but also the casting process of the resultant fragments, which is superior to both traditional finite element method and discrete element method.

4.1.4. Quantitative analysis on the crushed zone sizes Since the crushed zone plays an important role in controlling the effectiveness of the rock fragmentation by blast, the sizes of the crushed zone are obtained here using both the empirical equations and the hybrid finite-discrete element method and are compared with each other to quantitatively calibrate the hybrid finitediscrete element method. According to Esen et al. [78], the radius of the crushed zone can be calculated using Eq. (10):

rc ¼ 0:812r 0 ðCZIÞ0:219

ð10Þ

where r c is the crushing zone radius (mm), r0 is the borehole radius (mm), and CZI is the crushing zone index, which is a dimensionless index that identifies the crushing potential of a charged borehole. The CZI takes into account both explosive (borehole pressure) and rock (uniaxial compressive strength and stiffness) properties as well as borehole radius, which can be estimated using Eq. (11):

CZI ¼

ðPb Þ3 ðKÞ  r2c

ð11Þ

where Pb is the borehole pressure (Pa), K is the rock stiffness (Pa) and rc is the uniaxial compressive strength (Pa). If the borehole pressure is Pb ¼ 10 GPa, the crush zone index is 492.80 according

334

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) t = 0.15 ms

b) t = 0.25 ms

c) t = 0.35 ms

d) t = 0.45 ms

e) t = 0.55 ms

f) t = 0.75 ms

g) t = 1 ms

h) t = 1.5 ms

i) t = 2 ms

j) t = 2.5 ms

k) t = 3.5 ms

l) t = 4.5 ms

Fig. 6. Fracture initiation and propagation induced by rock blast in a rock mass with one free surface (Tensile crack is marked using blue colour while shear crack using red colour). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

to Eq. (11) using the input parameters listed in Table 2. Therefore, the radius of the crushed zone is 0.1530 m according to Eq. (10). As can be seen from Fig. 6, the crushed zone is formed immediately after the denotation. At 0.25 ms, the modelled crushed zone is approximately 0.17 m, which is close to 0.1530 calculated from the above empirical formula. The difference between the analytical and numerical results is about 10%, which indicates that the hybrid finite-discrete element method can well model not only the formation process of the crushed zone but also its size. It should be noted that, due to the expansion of the detonation gas, the formed crushed zone is expanded, too, as the borehole expands. However, the size of the formed crushed zone from the profile of the expanded borehole to that of the expanded crushed zone has not increased much.

4.1.5. Evolution of velocity, force and stress at specific monitoring points Fig. 7 depicts the variations of the velocity, force and stress at the specific monitoring points A, B and C located in the elastic zone, the cracked zone, and the crushed zone, respectively, induced by the single borehole blast in a rock mass with one free surface. At the time of 0 ms, the detonation is ignited. Since the detonation is put in the centre of the borehole, the explosive inside the

borehole is detonated with a velocity of V D ¼ 5500 m=s. Thus, it takes the time t = 0.05 m/5500 m/s  0.009 ms for the detonation to reach the borehole boundary. As soon as the explosive is detonated, high compressive stress wave is generated and propagates outward with a velocity of V P in Eq. (12)

VP ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ed 73:85 ¼ 5519:3 m=s ¼ qð1  v 2d Þ 2:6  106  ð1  0:262 Þ ð12Þ

where Ed is the dynamic Young’s modulus of rock, v d is the dynamic Poisson’s ratio, and q is the rock density. At the time m0:05 mÞ t ¼ 0:009 ms þ ð0:55 ¼ 0:0996 ms, the stress wave propa5519:3 m=s

gates to the monitoring point C, which is why the monitored velocity, force and stress at the monitoring point C reach the peak values at that time, shown in Fig. 7. Subsequently, at the times of 0.3989 ms and 0.6069 ms, the stress wave reaches the monitoring points B and A, respectively. Correspondingly, the monitored velocity, force and stress at these monitoring points reach their peak values, too. However, as shown in Fig. 7, their peak velocity, force and stress are much smaller than those monitored at the point C since the stress wave attenuates during its propagation and the monitoring point C is closer to the detonation borehole than the monitoring

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

points B and A. The comparison between the peak velocities, forces and stresses monitored at the points B and A shows the same trend. As the time increases, the stress wave further propagates and attenuates. Finally, as shown in Fig. 6, the velocity in the X-direction, force, and stress monitored at all points approach to zero, which indicates that the implementation of the non-reflection boundary into the hybrid finite-discrete element method has been proven to be effective. However, considerable velocities in the Y-direction are observed at all monitoring points, which is probably caused by the free surface at the top boundary of the model. As can be seen from Figs. 5 and 6, the blast-induced fragments are actually casted out of the free surface at the later stage of the rock blast, which is why the velocities in the Y-direction have not become zero as the stress wave attenuates. In short, the evolution of the monitored velocity, force and stress with respect to time follows that caused by the typical shock wave propagation. These variables first increase rapidly to the peak as the stress wave propagates to the monitoring point and then drop sharply to zero as the stress wave passes and attenuates. 4.2. Rock fracture process and resultant fragment casting by simultaneous and consecutive blasts In this section, a two-dimensional plane-strain model with two boreholes and one free surface is built to model the simultaneous and consecutive blasts. As depicted in Fig. 8(a), the space between the two boreholes is 2 m and the two boreholes are 2 m far from the free surface. In order to eliminate the influence caused by the reflected waves from the boundaries other than the free surface, the boundaries except the free surface are put far away from the borehole vicinity and are implemented as the non-reflection ones according to Section 3.4. As can be seen from Fig. 8(b), triangular elements are again used to discretize the model and the dense elements are employed in the vicinity of the borehole. The material properties of the rock mass follow those in the model for the single borehole blast. 4.2.1. Stress wave and fracture initiation and propagation Fig. 9(i) shows the stress wave propagation induced by the simultaneous blast. Before the stress waves induced by the two simultaneously detonated boreholes propagate to each other, the stress wave propagation radiated from each individual borehole (Fig. 9(i-a) is the same as that in the single borehole blast. The stress waves meet each other at around 0.25 ms (Fig. 9(i-b) and then interact with each other. Due to the superposition of the stress waves induced by the two individual borehole detonations, strong intensive stresses are observed in the vicinity between the two boreholes, as shown in Fig. 9(i-c). After that, it seems that the stress wave is initiated from a single combined borehole due to the close position of the two boreholes and propagate outward to cause the rock fracture and resultant fragment casting (Fig. 9 (i-d). Fig. 9(ii) depicts the stress wave propagation process induced by the consecutive blast with 1 ms delay time. Before the detonation of the second (right) borehole, the stress wave propagation (Fig. 9(ii-a, b) is the same as that in previous. At the time of 1 ms, the charge in the right borehole is ignited. As can be seen from Fig. 9(ii-c), a strong stress wave initiating from the second (right) borehole superposes that propagating outward from the first (left) borehole. Finally, the superposed stress wave and the merged detonation gas flow work together to result in the rock fracture and the eruption of the resultant fragments out of the free surface, as shown in Fig. 9(ii-d). Fig. 10(i) depicts the corresponding rock fragmentation process induced by the simultaneous blast. As soon as the boreholes are detonated, two individual crushed zones (Fig. 10(i-a) are formed

335

simultaneously around the two denoted boreholes. As the stress waves initiated from the two borehole detonations propagate to each other, the initiated cracks in the two crushed zones begin to interact with each other, too. When the compressive stress waves propagate to the free surface, they are then reflected above to form the tensile stress wave resulting in tensile failures in the vicinity of the free surface, as shown in Fig. 10(i-c). Meanwhile, the superposed stress waves propagate into the rock mass in other directions other then the top free surface, which work together with the merged detonation gas flow to generate a combined single cracked zone and long radial cracks (Fig. 10(i-d). Finally, the radial cracks are driven by the merged detonation gas flow to propagate out of the cracked zone and those propagating along the top direction coalesce with the tensile cracks near the free surface caused by the reflected stress wave to result in separate rock fragments, which are then casted out by the detonation gas, as shown in Fig. 10(i-e) and (f). The rock fracture process and the resultant fragment casting process induced by the consecutive blast are shown in Fig. 10(ii). After the detonation of the first (left) borehole, a crushed zone is formed around the borehole (Fig. 11(ii-a) as discussed previously in detail. After the time of 0.35 ms, the stress wave propagates to the second (right) borehole, which has not been detonated till this moment. As can be seen from Fig. 10(ii-b) and (c), the second borehole acts as another free surface since the compressive stress wave reflects there and the reflected tensile stress wave induces extensive tensile failures around the borehole profile. At the same time, tensile cracks are initiated near the top free surface due to the reflected tensile stress wave, too. At the time of 1.0 ms, the second (right) borehole is detonated and then a crushed zone is formed around the borehole. However, the detonation of the first (left) borehole creates an expanded opening, which actually works as another free surface for the detonation of the second (right) borehole. Moreover, the cracks induced by the detonation of the first (left) borehole greatly weaken the burdens in both directions toward the first (left) borehole and the originally designed top free surface. Correspondingly, the formed cracks around the second (right) borehole is no longer symmetrical. As shown in Fig. 10(iid), it is clear that the cracks in the left side of the borehole are better developed than those in the right side. As the stress wave propagates and the detonation gas flows, the cracked zone and the long radial cracks are well developed in the right side of the second borehole, too. However, due to the cracks caused by the detonation of the first (left) borehole, which greatly weaken the rock mass, the second (right) borehole is extremely expanded, as shown in Fig. 10 (ii-d). Finally, the rock fragments resultant from the coalescences of the cracks are casted out of the top surface. The comparison between the fragments caused by the simultaneous and consecutive blasts indicates that the consecutive blast generates more fragments than the simultaneous blast since more free surfaces are implemented/created in the consecutive blast. 4.2.2. Evolution of velocity and stress at specific monitoring points Similar to the case of the single borehole blast, three points A, B and C in the crushed, cracked, and elastic zones, respectively, to be formed during the simultaneous and consecutive blasts are monitored, as shown in Fig. 8(a). The monitored results are illustrated in Fig. 11 in terms of the relationship between the velocity and time and that between the stress and time. As can be seen from Fig. 11, the amplitude of the monitored peak velocity/stress and the time when the peak is monitored depend on the distance from the monitoring points and the detonation boreholes in both simultaneous and consecutive blasts. Basically, the further the monitoring point from the detonation borehole, the smaller the peak is and the later the peak appears, which follows the same variational law detailed in Section 4.1.5 although it is much more complicated due to the

336

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

100

100

X-velocity (m/s)

Y-velocity (m/s)

50

50

0

0

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -50

Time (ms)

-50

Time (ms)

-100

-100 -150

Point A in the elastic zone

-150

Point A in the elastic zone

-200

Point B in the cracked zone

-200

Point B in the cracked zone

Point C in the crushed zone

Point C in the crushed zone -250

-250

a) Velocity in the X-direction 15

b) Velocity in the Y-direction 15

X-Force (MN)

Y-Force (MN)

10

10 5

Time (ms)

0

5

Time (ms)

0

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -5

Point A in the elastic zone -10

Point B in the cracked zone Point C in the crushed zone

-15

-10 -15

c) Force in the X-direction 2000

X-stress (MPa)

d) Force in the Y-direction 2000

Time (ms)

0

Point A in the elastic zone Point B in the cracked zone Point C in the crushed zone

Y-stress (MPa)

Time (ms)

0

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2000

-4000

Point A in the elastic zone Point B in the cracked zone

-6000

-4000

Point A in the elastic zone -6000

-8000

Point B in the cracked zone Point C in the crushed zone

Point C in the crushed zone -8000

e) Stress in the X-direction

f) Stress in the Y-direction

Fig. 7. Evolution of velocity, force and stress at monitoring points A, B and C in the elastic zone, the cracked zone and the crushed zone, respectively, induced by rock blasting in a rock mass with one free surface.

Free surface

1m1m 2m 0.6 m

Point A

Point B

4.5 m

Point C

12 m

Non-reflection boundary

20 m

a) Geometrical model

b) Numerical model

Fig. 8. Geometrical and numerical models for simultaneous and consecutive blasts in rock masses with one free surface (The monitoring points A, B and C are placed in the elastic, cracked and crushed zones, respectively, to be induced by the rock blast).

complex interaction between the stress waves caused by the neighbouring borehole detonations in the either simultaneous and consecutive blasts. Moreover, it is noted from Fig. 11 that there are two peaks for each variable monitored in the consecutive blast while there is

only one peak in the simultaneous blast. That is because, in the simultaneous blast, the stress waves induced by the two borehole detonation arrive at the monitoring points at the same time since the monitoring points locate in the symmetrical line. However, in the consecutive blast, the 1 ms delay time causes the stress waves

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) t = 0.15 ms

a) t = 0.35 ms

b) t = 0.25 ms

b) t = 0.45 ms

c) t = 0.45 ms

c) t = 1.5 ms

d) t = 4.5 ms

d) t = 4.5 ms

i) Simultaneous blast

337

ii) Consecutive blast (1 ms delay)

Fig. 9. Stress wave propagation induced by simultaneous blast and consecutive blast with 1 ms delay time in rock masses with one free surface.

generated by the two borehole detonations to arrive at the monitoring points in different times and then lead to the two peaks illustrated in Fig. 11(ii). As the stress waves pass by the monitoring points and attenuate, the amplitude of the variable drops rapidly and finally approaches zero except the velocity. In the simultaneous blast, the velocity in the X-direction approaches zero, too, while that in the Y-direction keeps a considerable value since the resultant fragments are casted out of the top free surface along the Y-direction. However, in the consecutive blast, the velocity in the X-direction at the monitoring points A and B in the crushed and cracked zone keeps a considerable negative value (moves along the left direction), too, since the resultant fragments at these monitoring points are casted out along the X-direction, too, by the second (right) borehole detonation. 4.3. Rock fracture process and resultant fragment muck-piling caused by mining production blast This section applied the hybrid finite-discrete element method to model the rock fracture and fragmentation process and the

resultant irregular-shaped, deformable and fracture-able fragment muck-piling process caused by mining production blast in underground mining sub-level caving. As illustrated in Fig. 12(a), during the sub-level caving, horizontal sub-levels are created in the rock mass and ore passes extend vertically down connecting each sub-level to the main level. Access tunnels are excavated to run over the length of the ore body within a sub-level and crosscuts are drilled perpendicular to the access tunnels. From the crosscuts, rings of boreholes are drilled upward in a fan-shape pattern (middle picture in Fig. 12). Explosives are charged into the boreholes and the rings are blasted in sequence to recover the ore on each sublevel starting with the overlying sublevels and proceeding downwards. Within each sublevel, the ore is removed from the hanging wall to the footwall. As the ore is recovered from a sublevel, the hanging wall collapses by mine design. The mining production blast in a single sub-level in Fig. 12(a) can be simplified as a plane strain problem, as shown in Fig. 12 (b) with five rings in a single sub-level. Using the Swedish LKAB Malmberget mine as an example for the numerical model Fig. 12 (b), the spacing between two neighbouring rings is set as 3.5 m,

338

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) t = 0.35 ms

a) t = 0.35 ms

b) t = 0.55 ms

b) t = 0.55 ms

c) t = 1.0 ms

c) t = 1.0 ms

d) t = 2.0 ms

d) t = 2.0 ms

e) t = 3.5 ms

e) t = 3.5 ms

f) t = 4.5 ms

f) t = 4.5 ms

i) Simultaneous blast

ii) Consecutive blast (1 ms delay)

Fig. 10. Rock fracture processes induced by simultaneous blast and consecutive blast with 1 ms delay time in rock masses with one free surface.

the drift below the production ring is 5 m high and the height of the ring is set as 30 m. Each borehole in the five rings has a length of about 28 m and a diameter of 100 mm. The hybrid finite-discrete element modelling focuses on the rock fracture and the resultant

fragment muck-piling caused by the blasts in the five rings, which are discretized using dense finite elements. The other 12 m boundaries close to the production rings are discretized using coarse finite elements. In order to simulate the rock fragment

339

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

200 X-velocity (m/s) 150 100 50 Time (ms) 0 -50 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -100 Point A in the crushed zone -150 Point B in the cracked zone -200 Point C in the elastic zone -250

200 Y-velocity (m/s) 150 100 50 Time (ms) 0 -50 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -100 -150 Point A in the crushed zone Point B in the cracked zone -200 Point C in the elastic zone -250

a) Velocity in the X-direction 200 X-velocity (m/s) 150 100 Time (ms) 50 0 -50 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -100 Point A in the crushed zone -150 Point B in the cracked zone -200 Point C in the elastic zone -250

a) Velocity in the X-direction 200 Y-velocity (m/s) 150 Time (ms) 100 50 0 -50 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -100 Point A in the crushed zone -150 Point B in the cracked zone -200 Point C in the elastic zone -250

b) Velocity in the Y-direction 1000

1000

X-stress (MPa)

0 -1000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Time (ms)

-6000

-1000

Point B in the cracked zone Point C in the elastic zone

-4000 -5000 -6000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Time (ms)

-1000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Time (ms)

-2000

-3000

-5000

Y-stress (MPa)

0

-2000

-4000

Point A in the crushed zone Point B in the cracked zone Point C in the elastic zone

a) Stress in the X-direction 1000

X-stress (MPa)

0 -1000

Time (ms)

-3000 Point A in the crushed zone

a) Stress in the X-direction 1000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

-2000

-3000

-5000

Y-stress (MPa)

0

-2000

-4000

b) Velocity in the Y-direction

-3000 Point A in the crushed zone Point B in the cracked zone Point C in the elastic zone

-6000

-4000 -5000 -6000

b) Stress in the Y-direction i) Simultaneous blast

Point A in the crushed zone Point B in the cracked zone Point C in the elastic zone

b) Stress in the Y-direction ii) Consecutive blast (1 ms delay)

Fig. 11. Evolution of velocity and stress at monitoring points A, B and C in the crushed, cracked and elastic zones, respectively, induced by simultaneous and consecutive blasts.

muck-piling process, a 5 m high drift below the production ring is considered in the model and the bottom boundary of the drift is fixed in both horizontal and vertical directions. The right boundary of the model is regarded as a non-reflection boundary since the rock mass may be infinite in actual mining production blasting. It is assumed that the acceleration in the vertical direction is 100 times gravity in order to reduce the calculation time. Fig. 13 illustrates the modelled rock fragmentation process and resultant fragment muck-piling process. The borehole in the left ring is firstly ignited at the bottom of the borehole at 0 ms. Then the explosives in the borehole are burning upwards with a detonation velocity of 5500 m/s. As can be seen from Fig. 13(a) and (b), as soon as the borehole is detonated, high compressive stress waves are generated in the detonated area and propagate outwards with a wave velocity depending on the rock mass properties. At the

same time, the crushed zones are formed along the length of the detonated borehole. The stress wave propagation and the formation of the crushed zone are introduced in detail in Section 4.1, which are not repeated here although the detonation here happens along the full length of the borehole. After that, the explosives in the second ring are detonated at the bottom of the borehole and so on. Thus, the boreholes in the neighbouring rings are detonated consecutively at the bottom of the boreholes with the 5 ms delay time, as shown in Fig. 13(c)–(e). As the induced stress waves propagate, interact, reflect and attenuate, and the detonation gases flow, the cracked zone and the long radial cracks discussed in detail in Section 4.1 are well developed around the boreholes. Since the burden between the neighbouring ring is small and only 3.5 m, the cracks caused by the neighbouring borehole detonation actually interact and coalesce with each other resulting in the separa-

340

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

28 m

3.5 m

a) Mining production blasting in sublevel caving

5m

Non-reflection boundary

30 m

Free surface

b) Simplified numerical model for mining production blasting in 5 rings at one level

Fig. 12. Numerical model for rock fragmentation and resultant fragment muck-piling by mining production blasting in sublevel caving.

tion of the formed fragments from the infinite rock mass at the right boundary of the model. After the separation, the resultant fragments cave in due to their gravity Fig. 13(f). At about 100 ms, the caved rock fragments collide with the fixed bottom surface of the drift causing high stress concentrations are generated in the contact area and propagate upward Fig. 13(g), which results in the further fragmentation of the rock fragments, as shown in Fig. 13(h). Due to the detonations and the collisions, the formed rock fragments present complicated irregular shape and the irregular-shaped fragments further deform and fracture resulting in secondary fragmentations. Finally, the resultant irregular and deformable fragment interact with each other, further fracture if necessary, and muck-piling along the bottom surface of the drift until a steady-state is achieved, as shown in Fig. 13(i-j). 5. Discussions 5.1. Effects of explosive parameters on rock fragmentation by blasting The gas pressure induced by the detonation of typical commercial explosives can easily exceed 10 GPa. The very rapid conversion of the chemical energy into the thermal and gaseous energy may develop a detonation pressure from 0.5 to 50 GPa [16]. Correspondingly, as shown in Fig. 14, two hybrid modellings are curried out to model the rock fracture processes induced by the blast with different initial gas pressures in order to investigate their effects on the rock fragmentation process. The detonation with the high initial gas pressure generates a much (about 1300%) bigger crushed zone than that with the low initial gas pressure, as shown in Fig. 14(i). Moreover, a bigger (about 300%) cracked zone is caused by the detonation with the high initial gas pressure, too. It should be noted from Fig. 14(ii) that most of the cracks in the cracked zone are caused by compressive failures in the case of the detonation with the high initial gas pressure while most of them are tensile failures in the case of the detonation with the low initial gas pressure. Furthermore, as shown in Fig. 14(iii), the detonation with the high initial gas pressure induces much more and finer rock fragments but less long radial cracks. In addition, a bigger and deeper crater is formed near the free surface in the case of the detonation with the high initial gas pressure. These phenomena can be easily explained by the classical theories on the rock blast reviewed in Section 2. The rock fracture and fragmentation in the rock blast are mainly caused by the strong shock wave, detonation gas flow and reflected tensile wave. The strong shock wave is responsible to form the crushed zone,

which explains why a bigger crushed zone is formed and why there are much more and finer fragments in the case of the detonation with the high initial gas pressure since the high initial gas pressure results in a stronger shock wave than the low initial gas pressure. The shock wave attenuates as it propagates outward. Then, the detonation gas penetrates into the cracks and expands them to further propagate forming long radial cracks. When the compressive shock wave propagate to a free surface, it reflects back to form a tensile stress wave which usually generates tensile failures near the free surface since the tensile strength of rock is much smaller than its compressive strength. This explains why a bigger and deeper crater is formed in the case of the detonation with the high initial gas pressure since a strong compressive stress wave induced a strong reflected tensile stress at the free surface. Moreover, the bigger crushed zone formed in the case of the detonation with the high initial gas pressure can also be explained by Eqs. (10) and (11) introduced in Section 4.1.4. On the basis of Eqs. (10) and (11), the radius of the crushed zone is a function of the initial gas pressure and the higher the initial gas pressure, the bigger the crushed zone is.

5.2. Effects of rock parameters on rock fragmentation by blasting Besides the explosive parameters, the rock fragmentation process by blast is affected by the rock mass properties, too. In order to demonstrate the hybrid finite-discrete element method is able to capture the effect of rock parameters on the rock blast, three hybrid modellings are conducted for the blasts in the rocks with the different ratios between the compressive and tensile strengths (C/T ratio) while all other parameters are the same. The modelled final fracture patterns are depicted in Fig. 15. As can be seen from Fig. 15, the higher the C/T ratio, the longer (about 140% and 200% for intermediate and highest ratios, respectively) and the more (about 150% and 260% for intermediate and highest ratios, respectively) the radial cracks are while the smaller the crushed zone and the expanded borehole, and the finer the fragments in the crater are, which are consistent with the results obtained by Szuladzinski [79] using empirical equations. The higher C/T ratio means that the rock has a higher compressive strength and a lower tensile strength. The higher compressive strength indicates a smaller crushed zone and a less expanded borehole while the smaller tensile strength indicates more tensile failures in the cracked zone and near the free surface and longer radial cracks propagating out of the cracked zone, which are clearly illustrated in Fig. 15.

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) t = 1.25ms

b) t = 5ms

c) t = 6.25ms

d) t = 11.5ms

e) t = 16.25ms

f) t = 50ms

g) t = 100ms

h) t = 125ms

i) t = 175ms

j) t = 225ms

341

Fig. 13. Rock fracture process and resultant fragment muck-piling induced by rock blast in mining sub-level caving.

5.3. Effects of controlled blast on rock fragmentation During the rock blasting operations, it is important to control the fracture initiation and propagation. For example, in underground excavation, smooth blast technique is implemented

to reduce the excavation damage/disturb zone, avoid the overbreaks and achieve the desired smoothness. In order to demonstrate the application of the hybrid finite-discrete element method in controlled blast, the rock fracture process caused by the blasts in the boreholes with pre-split notches are modelled. Fig. 16 depicts

342

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

i) t = 0.5 ms

i) t = 0.5 ms

ii) t = 1.5 ms

ii) t = 1.5 ms

iii) t = 3.5 ms

iii) t = 3.5 ms

a) Low initial gas pressure (10 GPa)

b) High initial gas pressure (50 GPa)

Fig. 14. Effects of the explosive parameters (i.e. blast-induced initial gas pressure) on sizes of the crushed zone and the cracked zone.

a) Lowest C/T ratio

b) Intermediate C/T ratio

c) Highest C/T ratio

Fig. 15. Effects of the rock mass parameters (i.e. ratio between the compressive strength and the tensile strength) on the fracture pattern (t = 4.5 ms).

the modelled stress wave propagation and corresponding fracture process. As can be seen from Fig. 16(a), pre-split notches are drilled at the borehole surfaces to control the fracture initiation and propagation. After the instantaneous detonation, stresses concentrated around the boreholes and the pre-split notches, which results in the tensile extension of the pre-split notches and the compressive failures around the boreholes due to the stress wave propagation and the detonation gas flow, as shown in Fig. 16(b). As the stress waves propagate and the detonation gases flow, long radial cracks are initiated around the borehole and the tips of the pre-split notches propagating outward Fig. 16(c). At the time of 0.45 ms, the long radial cracks propagating out of the pre-split notches coalesce with each other, as shown in Fig. 16(d). After that, the stress

waves reaches the free surface and are reflected back to result in the tensile failure near the free surface depicted in Fig. 16(e). Finally, the long radial cracks propagating out of the pre-split notches further propagate in a direction sub-parallel to the free surface, as shown in Fig. 16(f). Therefore, the hybrid finitediscrete element modelling well demonstrates the effect of the pre-split notches in controlling the fracture initiation and propagation in the controlled blast. Further studies on the effect of multiple pre-existing cracks /joints are presented separately [80]. It should be noted that the explosive charge in the controlled blast must be much lighter than that in the normal blast so that the crushed zone is not formed around the borehole. Otherwise, the formed crushed zone may get rid of the effect of the pre-split

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) t = 0.15 ms

a) t = 0.15 ms

b) t = 0.25 ms

b) t = 0.25 ms

c) t = 0.35 ms

c) t = 0.35 ms

d) t = 0.45 ms

d) t = 0.45 ms

e) t = 1 ms

e) t = 1 ms

f) t = 4.5 ms

f) t = 4.5 ms

i) Stress wave propagation

343

ii) Fracture process

Fig. 16. Stress wave propagation and fracture induced by controlled blasting with pre-existing notch.

notches. Fig. 17 depicts the effect of the explosive charges on the rock fracture patterns induced by controlled blast. As can be seen from Fig. 17(a), with the light explosive charge (i.e. the initial gas pressure is 1 GPa), the blast-induced cracks are initiated from the tips of the pre-split notches and propagate approximately parallel to the original direction of the pre-split notches. With the explosive charge increasing, the effect of the pre-split notches becomes less and less important Fig. 17(b). With the high explosive charge, the crushed zone and the cracked zone are well developed around the borehole and the effect of the pre-split notches can be ignored as illustrated in Fig. 17(c).

6. Conclusions In this study, the hybrid finite-discrete element method is implemented to model the rock fracture process and the resultant fragment muck-piling process induced by two-dimension rock blast under various conditions. A state-of-the-art review is firstly conducted to clarify the rock fracture mechanism in rock blast and the advantages and disadvantages of the various continuous and discontinuous numerical methods in modelling the blastinduced rock fracture process. Then the hybrid finite-discrete element method is further developed to consider the transition from

344

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345

a) Low pressure: 1 GPa

b) Intermediate pressure: 5 GPa

 The hybrid finite-discrete element method has well modelled the transition of rock from continuum to discontinuum through the blast-induced fracture and even the whole process of the rock fracture and resultant irregular-shaped, deformable and further breakable fragment casting, and muck-piling caused by the rock blast. These make the hybrid finite-discrete element method is superior to the traditional continuum-based finite element method and discontinuum-based discrete element method reviewed in Section 2.2 since none of them can model the whole rock fracture and resultant irregular-shaped and further breakable fragment casting, and muck-piling process.  The hybrid finite-discrete element method is able to capture the dynamic fracture behaviour of rock under the blast-induced impact loads through implementing an empirical relation between the dynamic and static strengths and fracture energy releases rates derived from the dynamic rock fracture experiments.  The hybrid finite-discrete element method is a valuable tool for studying rock blasts since it has been demonstrated that the hybrid method can well model the stress wave propagation, reflection and attenuation, the detonation gas flow through the fracturing rock, the formation of the crushed, cracked and elastic zones, the long radial crack initiation and propagation, the tensile failure caused by the stress wave reflection at the free surface, and the effect of various properties of the explosive and rock on the rock fracture process in rock blast. Further study is needed to realistically model the fragment size distribution.

Acknowledgements The research presented in this study forms part of the first author’s PhD study at the University of Tasmania under the supervision of the corresponding author. The project is partly supported by a two-year visiting PhD scholarship provided by China Scholarship Council (CSC). The CSC’s support is greatly appreciated. Moreover, the authors would like to thank the anonymous reviewers for their valuable comments and constructive suggestion.

c) High pressure: 10 GPa References Fig. 17. Effect of the explosive charges on the rock fracture pattern induced by controlled blasting with pre-existing notches.

continuum to discontinuum, the detonation-induced gas expansion and flow through fracturing rock, the dependence of the rock dynamic behaviour on the loading rate, and the non-reflection boundary involved in the rock blast modelling. After that, the hybrid finite-discrete element method is implemented to model the rock fracture process and resultant fragment casting by a single borehole blast and the modelled crushed zone, cracked zone, and long radial crack propagation are compared with those well documented in literatures to calibrate the hybrid finite-discrete element method. Subsequently, the hybrid finite-discrete element method is applied to study the rock fracture process and the resultant fragment casting and even muck-piling process caused by the simultaneous blast, the consecutive blast, the mining production blast and the controlled blast. Finally, various factors affected the rock fragmentation process by the blasts are studied using the hybrid finite-discrete element method. Although there are several limitations such as simplified pressure-time relationship, small time integration step, and mesh-dependent fragment generation in this study, the results of this study conclude that

[1] Muller B. Adapting blasting technologies to the characteristics of rock masses in order to improve blasting results and reduce blasting vibrations. Int J Blast Fragmentation 1997;1(3):361–78. [2] Hamdi E, Mouza J, Fleurisson JA. Evaluation of the part of blasting energy used for rock mass fragmentation. Int J Blast Fragmentation 2001;5(3):180–93. [3] Hall JA, Brunton ID. Critical comparison of Julius Kruttschnitt Mineral Research Centre (JKMRC) blast fragmentation models. Fragblast 2002;6(2):207–20. http://dx.doi.org/10.1076/frag.6.2.207.8670. [4] Ouchterlony F, Olsson M, Kou SQ, Evertsson M, Wang YM, Forssberg E. Kartläggning av FoU inom fragmenteringsområdet i Sverige och omvärlden. GHRR project report; 2004. [5] Zhang ZX. Increasing ore extraction by changing detonator positions in LKAB Malmberget mine. Int J Blast Fragmentation 2005;9(1):29–46. [6] Kisslinger C. The generation of the primary seismic signal by a contained explosion. DTIC Document; 1963. [7] Kuznetsov V. The mean diameter of the fragments formed by blasting rock. J Min Sci 1973;9:144–8. [8] Cunningham C. The Kuz-Ram model for prediction of fragmentation from blasting. In: Proceedings of the first international symposium on rock fragmentation by blasting, Lulea, Sweden. p. 439–53. [9] Rossin P, Rammler B. The laws governing the fineness of powdered coal. J Instit Fuel 1933;7:29–36. [10] Cunningham C. Fragmentation estimations and the Kuz-Ram model-Four years on. In: Proceedings of the 2nd international symposium on rock fragmentation by blasting. p. 475–87. [11] Gheibie S, Aghababaei H, Hoseinie S, Pourrahimian Y. Modified Kuz—Ram fragmentation model and its use at the Sungun Copper Mine. Int J Rock Mech Min Sci 2009;46:967–73. [12] Djordjevic N. Two-Component model of blast fragmentation. Fragblast 1999:213–9.

H.M. An et al. / Computers and Geotechnics 81 (2017) 322–345 [13] Kanchibotla SS, Valery W, Morrell S. Modelling fines in blast fragmentation and its impact on crushing and grinding. In: Explo ‘99–A conference on rock breaking. Kalgoorlie, Australia: The Australasian Institute of Mining and Metallurgy; 1999. p. 137–44. [14] Ouchterlony F. Hur ser styckefallet efter sprängning ut? In: Bergsprängningskommittens 50:e – diskussionsmöte BK. p. 121–38. [15] Latham JP, Munjiza A, Lu P. Rock fragmentation by blasting—a literature study of research in the 19800 s and 19900 s. Fragblast 1999;3:193–212. [16] Saharan M, Mitri H, Jethwa J. Rock fracturing by explosive energy: review of state-of-the-art. Fragblast 2006;10:61–81. [17] Hino K. Fragmentation of rock through blasting and shock wave theory of blasting. In: The 1st US symposium on rock mechanics (USRMS). American Rock Mechanics Association; 1956. p. 191–207. [18] Bhandari S. On the role of stress waves and quasi-static gas pressure in rock fragmentation by blasting. Acta Astronaut 1979;6:365–83. [19] Daehnke A, Rossmanith H, Napier J. Gas pressurisation of blast-induced conical cracks. Int J Rock Mech Min Sci 1997;34(263):e1–e17. [20] Dally J, Fourney W, Holloway D. Influence of containment of the bore hole pressures on explosive induced fracture. Int J Rock Mech Min Sci Geomech Abstr 1975:5–12. [21] Hino K. Theory of blasting with concentrated charge. J Ind Explos Soc, Jpn 1954;15. S.233 ff. [22] Duvall WI, Atchison TC. Rock breakage by explosives. Bureau of Mines: US Department of the Interior; 1957. [23] Rinehart JS. On fractures caused by explosions and impacts. Colorado School of Mines; 1960. [24] Starfield A. Strain wave theory in rock blasting. In: The 8th US symposium on rock mechanics (USRMS). American Rock Mechanics Association; 1966. [25] Porter DD. Crater formation in plaster of Paris models by enclosed charges. Colorado School of Mines; 1961. [26] Clark LD, Saluja SS. Blasting mechanics. Trans Am Instit Min Eng 1964;229:78–90. [27] Saluja SS. Mechanism of rock failure under the action of explosives. In: The 9th US symposium on rock mechanics (USRMS). New York: American Rock Mechanics Association, AIME; 1967. Paper ID: ARMA-67-0297. [28] Langefors U, Kihlström B. The modern technique of rock blasting. Wiley; 1978. [29] Kutter H, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci Geomech Abstr 1971;8(3):181–202. [30] McHugh S. Crack extension caused by internal gas pressure compared with extension caused by tensile stress. Int J Fract 1983;21:163–76. [31] Preece D, Thorne B, Baer M, Swegle J. Computer simulation of rock blasting: a summary of work from 1987 through 1993. SAND94-1027. Sandia National Laboratories; 1994. [32] Thorne BJ, Hommert PJ, Brown B. Experimental computational investigation of the fundamental mechanisms of cratering. Brisbane, Australia: Fragblast’90; 1990. p. 117–212. [33] Grady DE, Kipp ME. Continuum modelling of explosive fracture in oil shale. Int J Rock Mech Min Sci Geomech Abstr 1980;17(3):147–57. [34] Taylor LM, Chen E-P, Kuszmaul JS. Microcrack-induced damage accumulation in brittle rock under dynamic loading. Comput Methods Appl Mech Eng 1986;55:301–20. [35] Kuszmaul JS. A new constitutive model for fragmentation of rock under dynamic loading. Albuquerque, NM (USA): Sandia National Labs; 1987. [36] Fourney W, Dick R, Wang X, Wei Y. Fragmentation mechanism in crater blasting. Int J Rock Mech Min Sci Geomech Abstr 1993;30(4):413–29. [37] 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. [38] Liu H, Kou S, Lindqvist PA, Tang CA. Numerical modelling of crack propagation by blast in TASQ tunnel of Aspo Hard Rock Laboratory. SKB project report, Sweden; 2005. [39] Rouabhi A, Tijani M, Moser P, Goetz D. Continuum modelling of dynamic behaviour and fragmentation of quasi-brittle materials: application to rock fragmentation by blasting. Int J Numer Anal Meth Geomech 2005;29:729–49. [40] Liu H, Kou S, Lindqvist P, Zhang Z. Numerical simulation of production blasting in sublevel caving. APS Blast 2007;2007:23–30. [41] Wang ZL, Li YC, Shen R. 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. [42] Wang Z, Li Y, Wang J. A method for evaluating dynamic tensile damage of rock. Eng Fract Mech 2008;75:2812–25. [43] Ma G, An X. Numerical simulation of blasting-induced rock fractures. Int J Rock Mech Min Sci 2008;45:966–75. [44] Wei X, Zhao Z, Gu J. Numerical simulations of rock mass damage induced by underground explosion. Int J Rock Mech Min Sci 2009;46:1206–13. [45] Sjöberg J, Schill M, Hilding D, Yi C, Nyberg U, Johansson D. Computer simulations of blasting with precise initiation. In: Proceedings of EUROCK 2012, Stockholm, Sweden. p. 12p. [46] Saharan MR, Mitri H. Numerical procedure for dynamic simulation of discrete fractures due to blasting. Rock Mech Rock Eng 2008;41:641–70. [47] Liu HY, Williams D, Pedroso D, Liang W. Numerical procedure for modelling dynamic fracture of rock by blasting. Controlling seismic hazard and sustainable development of deep mines: 7th international symposium on rockburst and seismicity in mines (RaSiM7), vols. 1 and 2. Rinton Press; 2009. p. 1551–60. [48] Zhu Z, Mohanty B, Xie H. Numerical investigation of blasting-induced crack initiation and propagation in rocks. Int J Rock Mech Min Sci 2007;44:412–24.

345

[49] Liang WM, Liu HY, Yang XL, Williams D. Effects of decoupled charge blasting on rock fragmentation efficiency. Int Soc Rock Mech, Int Congr 2011;2011:1237–40. [50] Preece D. A numerical study of bench blast row delay timing and its influence on percent-cast.. The 8th international conference of the international association for computer methods and advance in geomechanics, Georgetown, West Virginia, USA. 8p. [51] Song J, Kim K. Micromechanical modeling of the dynamic fracture process during rock blasting. Int J Rock Mech Min Sci Geomech Abstr 1996;33 (4):387–94. [52] Potyondy DO, Cundall PA, Lee C. Modelling rock using bonded assemblies of circular particles. In: Rock mechanics tools and techniques – proceedings of the 2nd north american rock mechanics symposium, Montreal, June 1996. Rotterdam: Balkema; 1996. p. 1937–44. [53] Donze F, Bouchez J, Magnier S. Modeling fractures in rock blasting. Int J Rock Mech Min Sci 1997;34:1153–63. [54] Munjiza A, Latham J, Andrews K. Detonation gas model for combined finite – discrete element simulation of fracture and fragmentation. Int J Numer Meth Eng 2000;49:1495–520. [55] Mortazavi A, Katsabanis P. Modelling burden size and strata dip effects on the surface blasting process. Int J Rock Mech Min Sci 2001;38:481–98. [56] Zhao XB, Zhao J, Cai JG, Hefny AM. UDEC modelling on wave propagation across fractured rock masses. Comput Geotech 2008;35(1):97–104. [57] Saiang D. Behaviour of blast-induced damaged zone around underground excavations in hard rock mass PhD Thesis. Sweden: Lulea University of Technology; 2008. [58] Ning Y, Yang J, An X, Ma G. Modelling rock fracturing and blast-induced rock mass failure via advanced dicretisation within the discontinuous deformation analysis framework. Comput Geotech 2011;38(1):40–9. [59] Ning YJ, Yang J, Ma GW, Chen PW. Modelling rock blasting considering explosion gas penetration using discontinuous deformation analysis. Rock Mech Rock Eng 2011;44(4):483–90. [60] Onederra IA, Furtney JK, Sellers E, Iverson S. Modelling blast induced damage from a fully coupled explosive charge. Int J Rock Mech Min Sci 2013;58:73–84. [61] Fakhimi A, Lanari M. DEM–SPH simulation of rock blasting. Comput Geotech 2014;55:158–64. [62] Liu HY, Kang YM, Lin P. Hybrid finite-discrete element modeling of geomaterials fracture and fragment muck-piling. Int J Geotech Eng 2015;9:115–31. [63] Yilmaz O, Unlu T. An application of the modified Holmberg-Person approach for tunnel blasting design. Tunn Undergr Space Technol 2014;43:113–22. [64] Aliabadian Z, Sharafisafa M, Mortazavi A, Maarefvand P. Wave and fracture propagation in continuum and faulted rock masses: distinct element modelling. Arab J Geosci 2014;7(12):5021–35. [65] Liu HY, Kou SQ, Lindqvist PA, Tang CA. Numerical studies on the failure process and associated microseismicity in rock under triaxial compression. Tectonophysics 2004;384:149–74. [66] Tang CA, Hudson JA. Rock failure mechanisms: illustrated and explained. ISBN 9780415498517. CRC Press; 2010. [67] Liu HY. A numerical model for failure and collapse analysis of geostructures. Aust Geomech 2010;45:11–9. [68] Munjiza A. The combined finite-discrete element method. Wiley Online Library; 2004. [69] Xiang J, Munjiza A, Latham JP. Finite strain, finite rotation quadratic tetrahedral element for the combined finite–discrete element method. Int J Numer Meth Eng 2009;79:946–78. [70] Tatone BSA, Grasselli G. A calibration procedure for two-dimensional laboratory-scale hybrid finite-discrete element simulations. Int J Rock Mech Min Sci 2015;75:56–72. [71] Mahabadi O, Kaifosh P, Marschall P, Vietor T. Three-dimensional FDEM numerical simulation of failure processes observed in Opalinus Clay laboratory samples. J Rock Mech Geotech Eng 2014;6:591–606. [72] Munjiza A, Andrews KRF, White JK. Combined single and smeared crack model in combined finite-discrete element analysis. Int J Numer Meth Eng 1999;44 (1):41–57. [73] Evans R, Marathe M. Microcracing and stress-strain curves for concrete in tension. Mater Struct 1968;1:61–4. [74] Zhao J. Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock. Int J Rock Mech Min Sci 2000;37:1115–21. [75] Zhang ZX, Kou SQ, Yu J, Yu Y, Jiang LG, Lindqvist PA. Effects of loading rate on rock fracture. Int J Rock Mech Min Sci 1999;36:597–611. [76] Lisjak A, Grasselli G. Combined finite-discrete element analysis of rock slope stability under dynamic loading. In: Proceedings of 2011 Pan-Am CGS Geotechnical Conference, Paper No: GEO11Paper473; 2011. [77] Lysmer J, Kuhlemeyer R. Finite dynamic model for infinite media. J Eng Mech Div ASCE 1969;95(EM4):859–76. [78] Esen S, Onederra I, Bilgin H. Modelling the size of the crushed zone around a blasthole. Int J Rock Mech Min Sci 2003;40:485–95. [79] Szuladzinski G. Response of rock medium to explosive borehole pressure. In: Proceedings of the fourth international symposium on rock fragmentation by blasting-Fragblast-4, Vienna, Austria. p. 17–24. [80] Liu HY, Han H, An HM, Shi JJ. Hybrid finite-discrete element modelling of asperity degradation and gouge grinding during direct shearing of rough rock joints. Int J Coal Sci Technol 2016. http://dx.doi.org/10.1007/s40789-0160142-1.