Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework

Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework

Computers and Geotechnics 38 (2011) 40–49 Contents lists available at ScienceDirect Computers and Geotechnics journal homepage: www.elsevier.com/loc...

3MB Sizes 0 Downloads 86 Views

Computers and Geotechnics 38 (2011) 40–49

Contents lists available at ScienceDirect

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

Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework Youjun Ning a,b, Jun Yang b, Xinmei An a, Guowei Ma a,c,⇑ a

School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China c School of Civil and Resource Engineering, University of Western Australia, 35 Stirling Highway, Crawley WA 6009, Australia b

a r t i c l e

i n f o

Article history: Received 20 May 2010 Received in revised form 21 September 2010 Accepted 21 September 2010 Available online 16 October 2010 Keywords: Advanced discretisation Rock mass failure Explosion gas pressure Rock blasting Discontinuous deformation analysis

a b s t r a c t Rock mass failure is a particularly complex process that involves the opening and sliding of existing discontinuities and the fracturing of the intact rock. This paper adopts an advanced discretisation approach to simulate rock failure problems within the discontinuous deformation analysis (DDA) framework. The accuracy of this approach in continuum analysis is verified first. Then, the advanced discretisation approach for fracturing modelling is presented, and the discretisation strategy is discussed. Sample rock static failures are simulated and the results are compared with experimental results. Thereafter, with a generalised definition of the artificial joints, this approach is further extended and applied in the simulation of blast-induced rock mass failures in which the instant explosion gas pressure obtained by the detonation pressure equation of state is loaded on the main blast chamber walls and the induced surrounding connected fracture surfaces. In the simulation instance of rock mass cast blasting, the whole process, including the blast chamber expansion, explosion gas penetration, rock mass failure and cast, and the formation of the final blasting pile, is wholly reproduced. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Unlike most engineering materials, rock mass is an inhomogeneous and anisotropic geological material consisting of both continuous rock medium and discontinuous components, such as joints, faults, and bedding planes. In rock mass failure problems, breaks occurring in relatively intact rock medium are known as cracks and fractures, and those occurring along discontinuous planes are known as discontinuity openings and slidings. In a rock mass containing a small number of discontinuities, cracking and fracturing may play a major role in its failure [1,2], while in a rock mass with a large proportion of discontinuities, opening and sliding along existing discontinuities will be more significant. Meanwhile, in rock failures associated with rock slopes [3] and underground excavations [4], breaks always appear in limited places and isolated rock blocks are produced and become unstable. However, in some other complex rock failure problems, such as rock blasting [5], a great number of breaks will occur, and the rock mass will be fractured into pieces and cast away. ⇑ Corresponding author at: School of Civil and Resource Engineering, University of Western Australia, 35 Stirling Highway, Crawley WA 6009, Australia. Tel.: +61 8 64883102; fax: +61 8 64881018. E-mail addresses: [email protected] (Y. Ning), [email protected] (J. Yang), [email protected] (X. An), [email protected] (G. Ma). 0266-352X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2010.09.003

To numerically simulate opening and sliding along discontinuities in a rock mass, discontinuous numerical methods, such as the distinct element method (DEM) [6] and discontinuous deformation analysis (DDA) [7], are always used. These methods were originally derived for discontinuous medium analysis and are applicable in modelling the behaviour of discontinuous blocky systems. By assembling discrete blocks together with a special algorithm to identify and update the contacts between the blocks, the detachment between blocks, block rotations, and the largescale displacement of individual blocks can be accurately predicted. To simulate fracturing in intact rock, continuum-based numerical methods are theoretically more applicable. In these methods, the stress distribution and the deformation of continuum can be calculated accurately. Then, by using fracture mechanics and/or other failure criteria, the process from continuum to discontinuum can be simulated; however, some difficulties may likely be encountered. In the popularly used finite element method (FEM), because the element edges need to coincide with the external boundaries and internal discontinuous planes of the physical model, great effort is required for the meshing and re-meshing process, which makes the realisation of the modelling procedure for complex failure problems complicated and sometimes even infeasible. Alternatively, to overcome the unexpected tedious meshing and re-meshing process, the conventional FEM has been modified using the partition of unity. The extended finite element method (XFEM)

41

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

[8] and the generalised finite element method (GFEM) [9] are two representatives. In these methods, the mesh does not need to conform to the crack surfaces anymore; however, the number of discontinuities that can be accounted for is limited, and large deformations and displacements of the block cannot be treated properly. Meanwhile, when fracture mechanics are adopted in the FEM, mathematical difficulties will be encountered in the crack branching and coalescent processes. Another typical continuous numerical method, namely, the boundary element method (BEM), has difficulty in treating material heterogeneity and nonlinearity. It is thus not suitable for solving problems in discontinuous and heterogeneous rocks. The meshless methods, such as the element-free Galerkin method (EFG) [10], are also promising for rock failure simulation, but they present difficulties in the numerical integration of weak forms and the imposition of essential boundary conditions. The application of discontinuous numerical methods to model fracturing is another option for simulating rock failure. One approach to implementing such a method is the continuous-discontinuous hybrid approach; basically this approach combines the FEM and the DEM together [11]. This approach is complicated to apply, and the difficulties that are encountered in the continuous methods will also be encountered in fracturing modelling. Another simplified approach is to model the fracturing of blocks based on approximate individual block stress states [12]. Because an intact block is broken into two blocks across the geometric centre once the failure criterion is fulfilled, the simulated fracturing route may be different from that observed in reality, especially when the block size is relatively large. An effective and practical approach to simulating rock failure problems with discontinuous numerical methods is to bond discrete elements or particles together and restrict failures to take place along the initially defined discontinuities based on certain failure criteria [13,14]. In the present study, by introducing artificial joints to divide a rock mass domain into small rock blocks beforehand and accompanying that with a joint strength reduction consideration during the fracturing along the joint planes, an advanced discretisation approach is adopted to simulate rock failure problems within the framework of the DDA method. After verifying this approach in continuum deformation analysis, the advanced discretisation approach is successfully extended for fracturing modelling. The advance discretisation strategy is discussed and a few examples are carried out to validate this fracturing modelling method in static problems. With a generalised definition of the artificial joints, the advanced discretisation approach is then further developed and applied in a rock blasting simulation in which the penetration effect of the explosion gas pressure is considered.

2. Continuum analysis via advanced discretisation 2.1. The DDA method The DDA method computes the static and dynamic behaviours of discrete blocky systems based on block kinematics. It gives a unique solution for large displacement and failure computations of block structures. The difference between static and dynamic computations is that in static problems, the block velocities are assumed to be zero at the beginning of each time step, while in dynamic problems the block velocities are inherited from the previous time step, either fully or partially, and the inherited block velocity proportion is defined as the dynamic parameter. Blocks of any shapes can be involved in a DDA block system, and a block system kinematics is developed to make blocks fulfill the constraints of no tension and no penetration between one another. Meanwhile, the Mohr–Coulomb criterion can also be applied be-

tween blocks. The basic theories of the DDA method can be found in Shi’s thesis [7]. Since it was first proposed, the DDA method has been widely used in rock mechanics and rock engineering problems, such as slope and landslide stability analysis [15–17], underground structure failures and reinforcement [18,19], and particle and masonry structure behaviours [20,21]. 2.2. Continuum analysis In the DDA method with first-order displacement functions, an individual block has constant stresses and constant strains throughout in each time step. A sub-block approach [12,22] can be utilised to improve the deformation ability and the stress distribution accuracy of a DDA block. In such an approach, an intact block is cut into some smaller blocks, called sub-blocks, and internal contact springs are applied between adjacent blocks in addition to the contact springs used in the original DDA method to guarantee displacement compatibility between blocks. In the present study, joints assigned with the same strength (three parameters including the friction angle, cohesion, and tensile strength are considered) of the intact block, called artificial joints here, are introduced to divide a continuous rock domain into small blocks beforehand to use the original DDA method to analyse the deformation and stress distribution of continuous rock medium. In this advanced discretisation approach, the displacement compatibility across the artificial joint planes is supposed to be satisfied by the bonding effect of the joint strength rather than the introduction of additional treatment techniques like that in the sub-block approach, i.e., the application of additional contact springs. The advantage of this approach is that the whole realisation procedure is within the frame work of the original DDA method and is easily implemented. Here, a cantilever model is used to verify the continuum analysis accuracy of this advanced discretisation approach. Accurate deformation and stress distribution analysis is still the precondition for the subsequent modelling from continuum to discontinuum. As shown in Fig. 1, the cantilever has dimensions of 8.0 m  1.0 m and is divided into 1340 triangular blocks. The origin is located at the bottom of the fixed end. A static downward point force of 5.0 kN is loaded at the point (7.9 m, 1.0 m). The density, Young’s modulus, and the Poisson’s ratio of the intact rock block are selected to be q = 2500 kg/m3, E = 35 GPa, and l = 0.3, respectively. The friction angle, cohesion, and tensile strength of the artificial joints between the discrete blocks are set as u = 70°, c = 30 MPa, and rt = 15 MPa, respectively. The calculated stresses in x direction (rX) at the upper boundary of the cantilever, the shear stresses (s) along the axis of the cantilever, and the cantilever deflection (f) are shown in Table 1. The corresponding theoretical solutions are also given in the table. At these given points, the maximum error in rX is 3.86% (at x = 6.0 m), and the maximum error in s is 2.67% (at x = 1.0 m). At x = 7.9 m, the error in the deflection is 1.42%. The numerically and the theoretically derived rX and f curves along the cantilever are shown in Fig. 2. It can be found that the simulated curves closely match the corresponding theoretical curves. A key factor that affects the continuum analysis accuracy is the contact spring stiffness applied between the artificial joint planes.

y

P

fixed

1m x O(0, 0)

8m Fig. 1. Cantilever model.

42

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

Table 1 Cantilever simulation result and theoretical solutions. x/m 5

1.0

2.0

3.0

4.0

5.0

6.0

7.0

7.9

rX/(10 Pa)

Theory DDA

2.070 2.020

1.770 1.739

1.470 1.424

1.170 1.144

0.870 0.863

0.570 0.548

0.270 0.272

– –

s/(103 Pa)

Theory DDA

7.50 7.30

7.50 7.35

7.50 7.35

7.50 7.36

7.50 7.37

7.50 7.36

7.50 7.44

– –

f/mm

Theory DDA

0.006 0.007

0.025 0.025

0.053 0.054

0.090 0.091

0.133 0.135

0.181 0.184

0.233 0.236

0.281 0.285

0.225

0.00

theory DDA

0.200 0.175

-0.10

0.150 0.125

f /mm

σX /MPa

theory DDA

-0.05

0.100 0.075

-0.15 -0.20 -0.25

0.050

-0.30

0.025 1

2

3

4

5

6

7

-0.35

0

1

2

3

4

5

6

7

8

x /m

x /m

(a)

(b) Fig. 2. rX and f curve comparisons.

Using the penalty contact method [7], the shear contact force at each angle-edge contact treatment between DDA blocks varies with different contact spring stiffness selections when the spring stiffness is relatively small; however, with a relatively large spring stiffness, stable and quite accurate contact forces can be obtained [23,24], and, meanwhile, the contact displacements can be neglected compared with the block deformations [7]. The fracturing along the artificial joints that is simulated later in this paper is also based on accurate contact forces between blocks. Because a toolarge spring stiffness will lead to an ill condition of the system equilibrium equations and increase the computer run-time significantly, a contact stiffness of 40E (E is the block Young’s modulus) is used in the simulation examples in this paper. Because the DDA method with the first-order displacement functions takes the small deformation assumption strictly in each time step, the selected time step displacement ratio (the allowed step maximum displacement throughout the analysed domain divided by the half of the length of the analysed domain in the y direction) should be small enough. A value between 0.001 and 0.01 is recommended by Shi (Shi GH. User’s manuals of discontinuous deformation analysis codes). In this paper, a time step displacement ratio smaller than 0.002 is used, and the step time is selected automatically by the DDA program itself based on the indicated time step displacement ratio. The above simulation example shows the accuracy of the advanced discretisation approach in continuum analysis, which indicates the possibility of extending this approach for the modelling of continuum failure. Of course, a large-enough block number is required to obtain a sufficiently accurate continuum analysis result. For an actual problem, the analysis results will converge by dividing the physical domain with an increasing number of small blocks, and the smallest block number required for the problem can be derived. In addition, second- or even higher-order displacement functions can help to improve the accuracy of the DDA method in the

bending problem analysis [25]. It can be deduced that less blocks will be required for a problem if higher-order displacement functions are applied; however, in the present paper, we conduct the study within the framework of the first-order DDA.

3. Fracturing modelling via advanced discretisation 3.1. Fracturing modelling method In the DDA method, the rigid body movement and deformation occur in each individual block while the opening and sliding occur between blocks. The failure of blocky systems between blocks can be calculated effectively. In the advanced discretisation approach, after dividing a continuous rock domain into small blocks with artificial joints, each discrete block is treated the same way as a regular block in the original DDA method, except that the artificial joints between the blocks are assigned with the same strength as that of the intact block to represent the continuous state across the joints, including both the continuity of the stress distributions and the displacements. With certain failure criteria, fracturing may take place along the pre-set artificial joints between blocks, and then the corresponding joint strength should be reduced. There are three types of contacts in the two-dimensional DDA method, including angle-to-angle contacts, angle-to-edge contacts, and edge-to-edge contacts [7], as shown in Fig. 3. In the original DDA method, the three contact types are all treated as contacts between angles and edges, and both tensile and shear contact forces can be obtained in each angle-edge contact treatment. Based on the angle-edge contact forces derived, in the fracturing judgment procedure of the advanced discretisation approach, the contact is directly processed as one of the above three contact types. Angleto-angle contacts are assumed to resist neither tensile nor shear contact forces; angle-to-edge contacts are assumed to resist shear

43

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

block i

block i

block i

block j

block j

block j

(a) angle-to-angle contact

(b) angle-to-edge contact

(c) edge-to-edge contact

Fig. 3. Three contact types in 2D-DDA.

but not tensile contact forces; edge-to-edge contacts resist both tensile and shear contact forces. When the analysed domain is divided into blocks in the style depicted in Fig. 1, adjacent blocks all take edge-to-edge and angle-to-angle contacts, or just angle-to-angle contacts before fracturing, as illustrated in Fig. 4. Angle-to-edge contacts can only appear after fracturing due to the deformation and movement of blocks. When two blocks (e.g., blocks i and m) take both edgeto-edge contacts (i.e., along P0P4) and angle-to-angle contacts (i.e., at P0), the fracturing between the two blocks can be identified based on the edge-to-edge contact (i.e., along P0P4) forces; on the other hand, when two blocks just take angle-to-angle contacts (e.g., blocks i and j), the two blocks will detach from one another after the fracturing takes to a group of blocks (e.g., block i to m and block m to j), which is also identified based on edge-to-edge contact (i.e., along P0P4 and P0P3) forces. Thus, all fracturing can be identified by edge-to-edge contact forces. For blocks i and m, after fracturing takes place along edge P0P4, the angle-to-angle contact at P0 may still exist; however, because angle-to-angle contact resists no contact forces, blocks i and m can detach each other freely when possible (e.g., after the fracturing between block i and n along P0P1). Both tensile and shear failures are considered in this advanced discretisation fracturing modelling approach by using the maximum tensile strength criterion and the Mohr–Coulomb criterion, respectively. The maximum tensile strength criterion is given by.

r ¼ rt

ð1Þ

and the Mohr–Coulomb criterion by

s ¼ c þ r tan u

ð2Þ

where r and s are the normal and shear contact stresses between the two sides of each pair of contacts, respectively, and rt, c, and u are the tensile strength, cohesion, and friction angle of the intact rock, respectively. Once any of the two criteria is fulfilled for a contact pair, a crack will initiate or propagate between the two sides of the contact. The corresponding artificial joints are then transformed into real joints, and the joint strength is reduced based on a proper joint strength model. In this fracturing modelling approach, the

P4 P1 block i block m

P0 block n block j

P3

P2

Fig. 4. Illustration of fracturing between blocks.

fracturing in both crack initiation and propagation is identified and processed with the same procedure based on the above failure criteria. Due to the deformation and movement of the pre-divided blocks, fracturing may take place between any two adjacent blocks. For a new fracture between two blocks, if there is no existing crack connected to this new crack, it is defined as an initiated crack; if there are one or more existing cracks connected to this new crack, it is a propagated or coalesced crack. In the present study, a simple joint strength reduction model is applied. After the fracturing of an artificial joint, the friction angle is supposed to be unchanged, the cohesion is reduced to a certain small value, and no tension will be endured anymore. In a real fracturing process, a closed crack may be able to endure shear forces due to the irregularity of the fracturing surfaces without normal compression. Using the advanced discretisation fracturing modelling approach with triangular small blocks, the crack between each pair of blocks remains straight and the irregularity within this crack length cannot be considered directly. Thus, a small cohesion, called a residual cohesion, is introduced to be applied to the generated crack surfaces in the strength reduction model to reflect the effect of the crack surface irregularity that exists in reality. The proper value of this residual cohesion needs special study to be derived, and, in the present study, an indicated value of 0.05 MPa is temporarily used in the simulation examples. 3.2. Advance discretisation strategy Because cracks can only initiate and propagate along pre-set artificial joints, the crack routes may be affected by the advance discretisation strategy. Regarding the angularity characteristics of two-dimensional polygons, triangles have the smallest average interior angles, which may have the least influence on the crack evolution paths, thus triangular blocks are recommended for the discretisation. Moreover, by dividing a continuous domain into enough small blocks, the unwanted influence of the discretisation strategy on the modelling result can become insignificant and negligible. Like that in the continuum analysis in Section 2, a numerical testing can be conducted for an actual problem. By dividing an analytical domain into more blocks with more reasonable local refinement, a convergent fracturing process and failure pattern can be obtained as the final result. Fig. 5 gives an edge fractured plate model under static uniaxial tension and its failure simulation result with several different discretisation strategies using the advanced discretisation fracturing modelling method. When pre-set artificial joints exist along the theoretical fracturing route, as shown in Fig. 5b, an ideal failure result can be obtained; otherwise, the fracturing route may be more or less different from the theoretical solution. As shown in Fig. 5c and d, the obtained fracturing route error is quite obvious. However, by dividing the region near the theoretical fracturing route into enough smaller blocks through global or local discretisation refinement, as shown in Fig. 5e and f, respectively, the simulated fracturing route converges to the theoretical solution.

44

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

(a) model

(d) strategy 3

(b) theoretical result (strategy 1)

(e) global refinement

(c) strategy 2

(f) local refinement

Fig. 5. Modelling of edge fractured plates under uniaxial tension with different discretisation strategies.

In the advanced discretisation fracturing modelling approach, the propagation of cracks along pre-set artificial joints between blocks might appear different from the cracks propagating in continuous rock. However, similar to that in the bonded-particle model (BPM) for rock fractures [13], based on the accurate computation of the stress distribution and the contact forces between blocks, the equivalence between the numerical modelling and linear elastic fracture mechanics (LEFM) can be established in terms of observable behaviours and mathematical descriptions. With the proper size of discrete blocks, the match of the fracture toughness is also possible to obtain. An extensive comparison between the advanced discretisation fracturing modelling approach and fracture mechanics will be carried out in future. Because the crack tip singularity that has been extensively discussed in fracture mechanics is naturally avoided in this advanced discretisation approach, even in particularly complex fracturing problems, including crack initiation, propagation, and coalescence, no technical or mathematical difficulties will be encountered.

Table 2 Joint mechanical parameters. Joint type

Friction angle (u/°)

Cohesion (c/MPa)

Tensile strength (rt/MPa)

Artificial joint Real joint

35 35

4 0.05

0.5 0

are applied through plates at both ends of the loaded diameter for each disc. The density, Young’s modulus, and the Poisson’s ratio of the intact rock are q = 2500 kg/m3, E = 10 GPa, and l = 0.2, respectively. The joint mechanical parameters are shown in Table 2. The pre-notched surfaces in Model (2) and the contact surfaces between the loading plates and the discs are smooth, without friction. The diametrical loading is chosen to be P = 40 kN in all the models. (1) The flattened intact disc

3.3. Simulation of rock sample failures To investigate the applicability of the advanced discretisation approach in the modelling of rock static fracturing problems, the failure process of several types of rock discs, including (1) a flattened intact disc, (2) a flattened pre-notched disc, and (3) two flattened discs each with an eccentric circular hole under diametrical pressure loadings are simulated, and the simulation results are discussed by comparing them with experimental results. Here, the comparison and validation are focused on the fracturing process and the final failure patterns, so all the discs are assumed to take the same diameter and physical parameters. The diameter of the discs is 5.0 cm and the flattened length is 0.16 cm (forming a 7.2° central angle at the disc centre). Loadings

Fig. 6 shows the simulated failure process of the intact disc model (divided into 1432 triangular blocks). Due to the tensile stress concentration at the disc centre perpendicular to the loaded diameter, a tensile crack is initiated there and propagates along the loaded diameter to the two ends. Scattered small cracks also appear near the loaded ends due to the stress concentrations there. The cracks finally break the disc into two semicircular halves, similar to the experimental results [26]. (2) The flattened pre-notched disc Fig. 7 shows the simulated failure process of the pre-notched disc model. The disc is divided into 1280 triangular blocks, and

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

45

(3) The flattened discs with eccentric circular holes

(a) model

(c) final failure

(b) crack initiation

(d) experimental result [26]

Fig. 6. Failure of a flattened intact disc.

the pre-notched crack length is 1.5 cm with an incline angle of 45°. Wing cracks first appear at the two ends of the pre-notched crack and then propagate towards the loaded ends under the effect of the tensile stress perpendicular to the crack route. The disc finally fractures into two halves, similar to the experimental results [27], but, due to the restrictions of the advance discretisation strategy, the actual crack route is a little different.

(a) model

(c) final failure

Fig. 8 is the configuration of a flattened disc with an eccentric circular hole [28]. The centre of the hole is located on the diameter perpendicular to the loaded diameter. As the hole radius a and its eccentricity e vary, the tensile stress gradients near the top and bottom of the hole will change. Although the maximum tensile stresses always appear near the top and bottom of the hole, due to the variety of their gradients, the dominant cracks may appear near the top and bottom of the hole or along the loaded diameter. Two actual models are simulated, including (3-I) in which a = 3.5 mm, e = 12 mm, where the disc splits along the loaded diameter and (3-II) in which a = 7 mm, e = 8 mm, where the disc fractures by intersecting the hole. Fig. 9 shows the simulated failure process of model 3-I (divided into 1502 triangular blocks). Because the hole is relatively small, the tensile stresses near the top and bottom of the hole decrease rapidly, so although small cracks appear there, the cracks fail to propagate far away. Along the loaded diameter, the existence of the hole makes the tensile stress concentration at the disc centre perpendicular to the loaded diameter relaxed, and the maximum tensile stress concentrations appear near the middle of the upper and lower radius, where cracks initiate. The propagation and coalescence of the cracks finally break the disc into two pieces along the loaded diameter. Due to the stress concentration near the loaded ends, scattered small cracks also appear there. All these results are consistent with the corresponding experimental results [28], except that the small offset of the failure route to the hole centre in the experiment is not revealed by the simulation due to the influence of the advance discretisation strategy. Fig. 10 shows the simulated failure process of model 3-II (divided into 1464 triangular blocks). Because the hole is relatively large, the maximum tensile stress concentrations near the top and bottom sides of the hole decrease slowly, so cracks are initiated there and propagate away. The propagation of the cracks towards the loaded ends finally breaks the disc into two parts. These results are also consistent with the corresponding experimental results [28]. In the above simulation examples, by comparing the simulation results of rock sample failures with related experimental results, it can be concluded that the static fracturing processes including crack initiation, propagation, and coalescence can be simulated relatively well with the advanced discretisation approach within the DDA framework.

(b) wing cracking

(d) experimental result [27]

Fig. 7. Failure of a flattened pre-notched disc.

Fig. 8. Flattened disc with a circular hole [28].

46

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

extended and applied in the simulation of a rock mass dynamic failure process in a rock blasting application. 4. Blast-induced rock mass failure simulation 4.1. Simulation method

(a) model

(b) crack initiation

(c) final failure

(d) experimental result [28]

Fig. 9. Failure of a flattened disc with a circular hole (3-I).

(a) model

(c) final failure

(b) crack initiation

(d) experimental result [28]

Fig. 10. Failure of a flattened disc with a circular hole (3-II).

The stress wave propagation characteristics study of DDA block systems shows that even with particularly high strength applied between blocks, stress waves attenuate rapidly while propagating across block boundaries [23], hence the advanced discretisation approach put forward above is more suitable for the static and quasi-static rather than dynamic modelling of continuous rock behaviours. Based on the proposed fracturing modelling method, in the following section, the advanced discretisation approach is further

4.1.1. Rock blasting simulation methods In the dynamic failure process of rock blasting, rock masses may be fractured and cast under the effect of explosion shock waves and explosion gas pressure. In a continuous rock mass, the explosion shock waves play a major role in the rock mass failure. In a jointed rock mass, because the stress waves attenuate rapidly from being reflected by discontinuities, the explosion gas pressure plays the dominant role in the rock mass fragmentation and rock block cast. Continuous numerical methods, such as the FEM, which employ damage constitutive models and erosion algorithms, are usually used to simulate the damage and fracture of the rock mass near the blast chamber to study the expansion of the blast chambers [29] and the fragmentation mechanisms of the rock mass under blasting loading [30]. For the simulation of rock mass failure, rock block cast and the formation of blasting piles, discontinuous numerical methods, however, have obvious advantages and are more applicable. Preece [31] studied the bench blasting row delay timing and its influence on percent-cast using a computer program known as the Distinct Motion Code (DMC). Minchinton [14] simulated the process from continuous to discontinuous rock medium in rock blasting successfully using the distinct element-based program MBM-2D. Mortazavi [32] simulated the rock block movement of the fragmented rock mass in rock-bench blasting using the DDA method considering the block’s rigid-body damping and the explosion gas penetration. These studies indicated the advantage of discontinuous numerical methods in rock blasting simulations, although the models were too idealised and the results deviated in a certain degree from the real process in rock blasting engineering. 4.1.2. Advanced discretisation modelling of blast-induced rock mass failure Due to the rapid attenuation of stress waves propagating in DDA block systems, the traditional DDA method essentially can only be used in the simulation of jointed rock mass blasting but not in continuous rock mass blasting [23]. In the present study, a generalised artificial joint concept is introduced to enhance the applicability of the DDA method in rock blasting simulation, in which the explosion gas pressure loading is applied directly and the effect of the explosion shock waves is considered indirectly by the discretisation of the rock mass beforehand. In such a method, the effect of the explosion shock waves and the explosion gas pressure are to be considered separately, and the explosion gas pressure is assumed to be acted on the jointed rock mass that was created previously by the explosion shock waves. This assumption is based on the fact that the explosion shock waves travel much faster than the explosion gas. Because the analysed domain is already a jointed rock mass, its behaviour under the explosion gas pressure loadings and the propagation of the induced stress waves in it can be simulated by the DDA method. In such a way, not only can jointed rock mass blasting be simulated, but also the blasting process in continuous rock mass can be simulated approximately. In an initial rock blasting model, discontinuities that do not violate the failure criteria are all defined as artificial joints, namely the generalised artificial joints, which may include both the discontinuities that are introduced for the advance discretisation of the rock mass domain, and the actual existing natural discontinuities with strength yet meeting with the failure criteria. The blastinduced rock mass failure along the generalised artificial joints is

47

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

simulated with the similar fracturing modelling algorithm and criteria that are used in Section 3.1. However, here the mechanical parameters of the artificial joints are selected according to actual conditions, but not necessarily to be the same as that of the continuous rock material like that in the modelling of continuum fracturing. After failure, an artificial joint will also be transformed into a real joint (or crack/fracture), and the joint strength is reduced with the same method that is used in Section 3.1. Discontinuities violating the failure criteria in the initial model are treated as cracks/ fractures directly with actual low joint strength.

RP roc rock ock mass maass

co connected onn onn neccted nec ctedd fractures fracctu frac ctuure uress

R0

main chamber

unccon unc unconnecte connne nneecte ect d fractures fract ffr aacttur ac turres res

4.1.3. Explosion gas pressure loading Fig. 11 shows a rock blasting model with one column borehole, as viewed along the borehole axis direction. The left, right, and lower boundaries are viscous non-reflection boundaries [33,34] and the upper boundary is a free surface. Taking rock mass gravity into consideration, this model can be comprehended as horizontalcolumn-borehole cast blasting. The advance discretisation of the rock mass is carried out by dividing the rock mass domain into triangular rock blocks, while the discretisation refinement is done around the borehole to account for the intensive fragmentation of the rock mass under the effect of the explosion shock waves. Assuming that the joints among the discrete blocks have sufficient strength that they do not yet meet the failure criteria, these joints are the generalised artificial joints. Under the effect of the explosion gas pressure, the initial bore hole will expand to form the main blast chamber and fractures will appear in the surrounding rock mass. Based on the fracturing modelling approach described in Section 4.1.2, the explosion gas is supposed to expand in the main blast chamber and penetrate into the surrounding connected fractures with a constant velocity, which we denoted as VP here. Fig. 12 is a sketch of the bore hole expanding as the blast chamber, where R0 is the initial borehole radius and RP is the dispersing radius of the explosion gas at time t. Then

RP ¼ V P t

ð3Þ

Thus, the blast chamber consists of the main chamber and the surrounding fractures connected to the main chamber within the range of RP. In each time step, the blast chamber profile can be searched by starting from the block boundaries along the blast chamber wall at the end of the previous time step to determine the newly generated connected fractures in the current time step according to crack surface connectivity, then the blast chamber profile is recorded for the application of the explosion gas pressure loadings and the search of the blast chamber in the next time step. As revealed in Fig. 12, the blast chamber may have an especially complicated geometrical shape, so the cross-sectional blast chamber area is derived indirectly by subtracting the area that does not belong to the chamber from the area within the range of RP, i.e., the area of the rock blocks, the artificial joints that haven’t been fractured, and the fractures that are not connected to the main chamber within RP is to be subtracted. As polygons, the area of each block, artificial joint, and fracture is precisely calculated by the simplex integration [7].

rocck mass rock mas ass

Fig. 12. Blast chamber expanding sketch.

Taking a plane stress assumption, the blast chamber has a unit length depth in the axis direction, so the blast chamber volume is equal to the cross-sectional blast chamber area in quantity. The instant explosion gas pressure in the blast chamber is calculated based on the detonation pressure equation of state

P ¼ P0

 c V0 V

ð4Þ

where P0 and V0 are the initial blast chamber pressure and volume, respectively, P and V are, respectively, the blast chamber pressure and volume at time t, and c is a constant related to the properties of the charged explosive. Once the explosion gas disperses to the free surface and leakage takes place, the loading of the explosion gas pressure vanishes. In each time step during the application of the explosion gas pressure, the instant blast chamber pressure is loaded on all the block boundaries along the blast chamber wall. In the DDA method with first-order displacement functions, the block boundaries always remain straight. The explosion gas pressure is applied as line loading [7,23] on each corresponding block boundary perpendicularly towards the block interior uniformly. As illustrated in Fig. 13, the pressure is loaded on main blast chamber wall, the connected fracture surfaces within RP, and the boundaries of isolated blocks inside the blast chamber. 4.2. Cast blasting simulation The cast blasting model shown in Fig. 11 has dimensions of 15.0 m  5.0 m. The initial borehole radius is R0 = 0.1 m. The

RP

aarti ar artificial rrtiific ific ficial fi cial joint join jo nt

P

P

expl eex explosion xp x plosi pl p lo los lo osi os ssiio on n gas gaaass g

P

P rocck mass rock mas ass

P

borehole

Fig. 11. Cast blasting model.

isola is olaate ateed ed block blo blo ockk ock P isolated

Fig. 13. Explosion gas pressure loading illustration.

48

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

Table 3 Joint mechanical parameters. Joint type

Friction angle (u/°)

Cohesion (c/MPa)

Tensile strength (rt/MPa)

Artificial joint Real joint

45 45

2 0.05

0.5 0

burden (distance from the borehole centre to the upper free surface) is w = 3.0 m. The initial blast chamber pressure is assumed to be P0 = 1.0 GPa, and we select VP = 100 m/s, and c = 1.4. The density, Young’s modulus, and Poisson’s ratio of the rock material are chosen to be q = 2600 kg/m3, E = 55 GPa, and l = 0.3, respectively. The joint mechanical parameters are shown in Table 3. A dynamic parameter adjustment [23,34] for DDA dynamic computations is also applied in the simulation, in which the dynamic parameter was improved as a function of the step time and block velocity. Fig. 14 gives the blast chamber volume expanding and pressure attenuating histories. Under the effect of the explosion gas pressure, the blast chamber begins to expand, and the failure of the rock mass around the blast chamber makes it easier for the chamber to expand further, so the chamber expands faster and faster as time progresses. The blast chamber pressure attenuates especially quickly at the beginning. After about 5.0 ms, the attenuating speed becomes much smaller. At 26.57 ms, the explosion gas expands to the upper free surface and leaks, so the pressure drops to zero, and its loading on the rock mass stops. At that time, the blast chamber volume is 4.073 m3, which is 131.4 times the initial borehole volume (0.0031 m3), and the blast chamber pressure correspondingly decreases from 1.0 GPa to 1.04 MPa. Fig. 15 gives several instances within the simulated blasting process. The two big blocks at the two sides of the model are used to receive the blasting pile. Fractures appear around the borehole under the effect of the explosion gas pressure first, and the dispersion of the explosion gas along fractures accelerates the further failure of the rock mass. Fig. 15a shows the explosion gas dispersing state at 0.02 s (joints in red represent where the explosion gas has arrived). At that moment, many fractures have been generated, and an obvious bulge can be seen from the free surface. After the leakage of the explosion gas (or the disappearing of the explosion gas pressure loading), the rock mass fractures further, and rock blocks are thrown up under the inertia effect. At 0.74 s, the cast blocks get the maximum height and then begin to fall back to form the final blasting pile.

0.8

chamber volume /m

3

chamber pressure /GPa

4

1

(b)

= 0.74 s

(c)

= 1.65 s

(d)

= 3.62 s

In the above simulation example, the blast chamber volume expanding and the explosion gas pressure attenuating histories are obtained based on the proposed blast chamber calculating procedure and the application method of the explosion gas pressure loadings. The rock mass fracturing process is represented. The movement of the generated rock blocks and the formation of the final blasting pile are wholly reproduced. Compared with the crater modelling result by the FEM employing the erosion algorithm [29], in which the elements of the fractured rock mass were removed from the model, the simulation result in the present study shows an obvious advantage of the discontinuous numerical methods in

1.0

2

= 0.02 s

Fig. 15. Simulated cast blasting process.

5

3

(a)

0.6 0.4 0.2 0.0

0 0

5

10

15

20

25

30

0

5

10

15

20

25

time /ms

time /ms

(a) chamber volume

(b) chamber pressure

Fig. 14. Blast chamber volume and pressure histories.

30

Y. Ning et al. / Computers and Geotechnics 38 (2011) 40–49

the modelling of blast-induced rock mass cast and the formation and characteristics of the blasting piles, which are of great importance in rock blasting engineering, e.g., in rock quarrying by traditional drill and blast methods. This simulation instance indicates the capability of the newly developed DDA method in rock blasting simulation. More parametric and comparison studies need to be carried out subsequently to validate this new method further. 5. Conclusions and remarks In this paper, by dividing a rock mass domain into small rock blocks with artificial joints beforehand and considering the joint strength reduction during the fracturing along artificial joint planes, rock fracturing and blast-induced rock mass failures are simulated via an advanced discretisation approach within the framework of the DDA method. The continuity across the artificial joint planes is guaranteed by the bonding effect of the joint mechanical parameters in the continuum analysis and accurate deformation and stress distribution results are obtained. By the global or local refinement of the advance discretisation, the influence of the discretisation strategy on the continuum fracturing modelling result can become negligible. Comparisons between the rock sample static failure simulation results and the experimental results indicate that the crack initiation, propagation, and coalescence processes can be predicted with the advanced discretisation fracturing modelling approach without encountering any technical or mathematical difficulties. Using a similar fracturing modelling algorithm with a generalised definition of the artificial joints, the advanced discretisation approach is then extended for the simulation of blast-induced rock failures. By tracking the blast chamber profile dynamically, the instant blast chamber volume is calculated and the explosion gas pressure that is obtained from the detonation pressure equation of state is loaded on the main blast chamber walls and the surrounding connected fracture surfaces. In the cast blasting simulation instance, the fracturing process in the rock mass is represented; thereafter, the movement of the generated blocks and the formation of the final blasting pile continue to be reproduced successfully. In addition to providing a new fracturing modelling technique, the fracturing modelling examples in this paper also verify the ability and correctness of the DDA method in the failure treatment of complex block systems. A further refinement of the advance discretisation method can help to gain more realistic fracturing modelling results. This work can be accomplished due to the optimisation of the program and the rapid development of computer technology. After the validation of the blasting simulation method proposed in the present study, more practical geological factors, such as the actual distribution of discontinuities, can be considered so that more realistic rock blasting models can be simulated. References [1] Wong LNY, Einstein HH. Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int J Rock Mech Min Sci 2009;46(2):239–49. [2] Wong LNY, Einstein HH. Crack coalescence in molded gypsum and carrara marble: part 1 – macroscopic observations and interpretation. Rock Mech Rock Eng 2009;42(3):475–511. [3] Hoek E, Bray JW. Rock slope engineering. 3rd ed. London: Institution of Mining and Metallurgy; 1981.

49

[4] Hoek E, Brown ET. Underground excavation in rock. London: Institution of Mining and Metallurgy; 1980. [5] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8(3):181–8. [6] Cundall PA. A computer model for simulating progressive, large-scale movements in blocky rock systems. In: Proceedings of the international symposium on rock mechanics. Nancy, France; 1971. [7] Shi GH. Discontinuous deformation analysis – a new numerical model for the static and dynamics of block systems. Berkeley: UC Berkeley; 1988. [8] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46(1):131–50. [9] Duarte CA, Reno LG, Simone A. High-order generalized FEM for through-thethickness branched cracks. Int J Numer Meth Eng 2007;72(3):325–51. [10] Belytschko T, Lu YY, Gu L. Element-free Galerkin method. Int J Numer Meth Eng 1994;37(2):229–56. [11] Munjiza A, Owen DRJ, Bic´anic´ N. A combined finite-discrete element method in transient dynamics of fracturing solids. Eng Comput 1995;12(2):145–74. [12] Lin CT, Amadei B, Jung J, Dwyer J. Extensions of discontinuous deformation analysis for jointed rock masses. Int J Rock Mech Min Sci Geomech Abstr 1996;33(7):671–94. [13] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [14] Minchinton A, Lynch PM. Fragmentation and heave modeling using a coupled discrete element gas flow code. FFAGBLAST 1997;1(1):41–57. [15] Hatzor YH, Arzi AA, Zaslavsky Y, Shapira A. Dynamic stability analysis of jointed rock slopes using the DDA method: King Herod’s Palace, Masada, Israel. Int J Rock Mech Min Sci 2004;41(5):813–32. [16] Wu JH, Wang WN, Chang CS, Wang CL. Effects of strength properties of discontinuities on the unstable lower slope in the Chiu-fen-erh-shan landslide. Eng Geol 2005;78(3–4):173–86. [17] Kveldsvik V, Einstein HH, Nilsen B, Blikra LH. Numerical analysis of the 650,000 m2 åknes rock slope based on measured displacements and geotechnical data. Rock Mech Rock Eng 2009;42(5):689–728. [18] Kim Y, Amadei B, Pan E. Modeling the effect of water, excavation sequence and rock reinforcement with discontinuous deformation analysis. Int J Rock Mech Min Sci 1999;36(7):949–70. [19] Tsesarsky M, Hatzor YH. Tunnel roof deflection in blocky rock masses as a function of joint spacing and friction-A parametric study using discontinuous deformation analysis (DDA). Tunn Undergr Space Technol 2006;21(1):29–45. [20] Pearce CJ, Thavalingam A, Liao Z, Bic´anic´ N. Computational aspects of the discontinuous deformation analysis framework for modelling concrete fracture. Eng Fract Mech 2000;65(2–3):283–98. [21] Bic´anic´ N, Stirling C, Pearce CJ. Discontinuous modelling of masonry bridges. Comput Mech 2003;31(1–2):60–8. [22] Cheng YM, Zhang YH. Rigid body rotation and block internal discretization in DDA analysis. Int J Numer Anal Methods Geomech 2000;24(6):567–78. [23] Ning YJ. Study on dynamic and failure problems in DDA method and its application. Beijing: Beijing Institute of Technology; 2008. [24] Ning YJ, Yang J, Ma GW, Chen PW. Contact algorithm modification of DDA and its verification. In: Ma GW, Zhou YX, editors. Analysis of discontinuous deformation – new developments and applications. Singapore: Research Publishing; 2009. p. 73–81. [25] Grayeli R, Mortazavi A. Discontinuous deformation analysis with second-order finite element meshed block. Int J Numer Anal Meth Geomech 2006;30(15): 1545–61. [26] Liu HY. Numerical modelling of the rock fragmentation process by mechanical tools. Luleå, Sweden: Luleå University of Technology; 2004. [27] Al-Shayea Naser A. Crack propagation trajectories for rocks under mixed mode I–II fracture. Eng Geol 2005;81(1):84–97. [28] Van De Steen B, Vervoort A, Napier JAL. Observed and simulated fracture pattern in diametrically loaded discs of rock material. Int J Fract 2005;131(1):35–52. [29] Wang ZL, Li YC, Shen RF. Numerical simulation of tensile damage and blast crater in brittle rock due to underground explosion. Int J Rock Mech Min Sci 2007;44(5):730–8. [30] Zhu ZM, Xie HP, Mohanty B. Numerical investigation of blasting-induced damage in cylindrical rocks. Int J Rock Mech Min Sci 2008;45(2):111–21. [31] Preece DS. A numerical study of bench blast row delay timing and its influence on percent-cast. In: Proceedings of the eighth international conference of the international association for computer methods and advances in geomechanics (IACMAG). West Virgina University, USA; 1994. [32] Mortazavi A, Katsabanis PD. Modelling burden size and strata dip effects on the surface blasting process. Int J Rock Mech Min Sci 2001;38(4):481–98. [33] Jiao YY, Zhang XL, Zhao J, Liu QS. Viscous boundary of DDA for modeling stress wave propagation in jointed rock. Int J Rock Mech Min Sci 2007;44(7): 1070–6. [34] Ning YJ, Yang J, An XM, Ma GW. Simulation of blast induced crater in jointed rock mass by DDA method. Front Archit Civ Eng China 2010;4(2):, 223–232.