Underground rockfall stability analysis using the numerical manifold method

Advances in Engineering Software 76 (2014) 69–85

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

Underground rockfall stability analysis using the numerical manifold method Zhijun Wu a,b, Louis Ngai Yuen Wong a,⇑ a

School of Civil and Environmental Engineering, Nanyang Technological University, Block N1, 50 Nanyang Avenue, Singapore 639798, Singapore Department of Civil Engineering, Zhejiang University, Yuhangtang Road, Hangzhou 310058, China

b

a r t i c l e

i n f o

Article history: Received 20 November 2013 Received in revised form 16 May 2014 Accepted 1 June 2014 Available online 1 July 2014 Keywords: XFEM Numerical manifold method Rockfall Rockbolt Shotcrete Layered rock Rockbolt reinforcement

a b s t r a c t The present study applies the numerical manifold method (NMM) as a tool to investigate the rockfall hazard in underground engineering. The crack evolution technique with crack initiation and propagation criterion, which has been successfully applied to handle cracking problems in rocks, is used in this study. A rockbolt element is introduced, which is first validated by a simple case. The mechanism of the rockbolt in reinforcing a layered rock mass is then investigated through a four-layered rock beam example. The developed NMM is then used to investigate the rockfall instability caused by either natural joints or mining induced fractures in an underground power station house or a tunnel. The results illustrate that the developed NMM can not only capture the entire dynamic process of the rockfall but also locate the keyblock successfully. As such, corresponding reinforcement methods can be chosen reasonably. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Rockfalls pose a great hazard to underground mine workers and infrastructure. Potentially unstable rock blocks, which are defined by natural joints or mining induced fractures, are commonly found in the roof and sidewall of excavations. The instability of those blocks, referred to as keyblocks, is driven by two successive mechanisms: the detachment of a block or rock compartment (i.e. failure) and its propagation [1]. The assessment of the crack propagation and failure within the rock mass is of main significance for a rigorous rockfall hazard assessment. In the past decades, various methods have been proposed to help practitioners analyse and assess the susceptibility to rockfall hazard [2–5]. These studies have highlighted the intensity and frequency of events as the key parameters. However, detailed accounts of the entire failure processes, especially fracture related failure, which involves the crack initiation, propagation, and coalescence, were not always available. Due to the development of the computer technology, numerical methods have become promising tools for capturing the details of the progressive fracturing. However, most of the numerical methods used in geotechnics have an implicit representation of the discontinuities, which only

⇑ Corresponding author. Tel.: +65 6790 5290; fax: +65 6791 0676. E-mail address: [email protected] (L.N.Y. Wong). http://dx.doi.org/10.1016/j.advengsoft.2014.06.001 0965-9978/Ó 2014 Elsevier Ltd. All rights reserved.

consider their influences on the physical behavior, such as deformability or strength, through equivalent constitutive laws. Since the first introduction of joint elements into FEM by Goodman [6], continuum-based methods have been progressing to deal with fracture propagation problems. Based on the Partition of Unity method (PU) [7], the Extended Finite Element Method (XFEM) [8,9] is one of the most advanced achievements in this field. However, the continuous description is still deficient in handling problems including a large sliding and a large amount of fractures. An alternative to these continuum-based methods is to use discrete-based methods (DEM) which treat the material as an assemblage of independent elements. By incorporating artificial discontinuities and corresponding constitutive failure models, both the explicit DEM UDEC [10] and implicit DEM DDA [11] were successfully used to model failures in either intact rock or jointed rock mass [12–17]. However, the cracks in these methods can only develop along the block boundaries. Therefore, crack trajectories predicted by them are not arbitrary, but depending on the block discretization configuration (mesh configuration). The Particle Flow Code (PFC) by Itasca [18,19] can model the fracturing process occurring within the rock mass without assuming where and how cracks may appear. By simply breaking the bond between two distinct elements when the interaction force overcomes its tensile or shear strength, the PFC can model the arbitrary fracturing process [20–22]. However, as mentioned by Donze et al. [23], the element size and element shape used in the model both differ from those of

70

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

real constituent grains in rocks, which give a cross effect on the simulation results. Although both continuum-based and discontinnum-based methods provide useful means to analyse rockfall hazard problems, modeling complex failures related to both detachment along pre-existing discontinuities and the creation of new ones through fracturing of intact rock are still challenging. Actually, the rock mass is neither continuous nor discontinuous but the integration of two, which can be better represented by hybrid methods. The Numerical Manifold Method (NMM) [24,25], which is also based on the PU concept, is such a hybrid method. It combines the widely used continuum-based method FEM and discontinuum-based method DDA in a uniform framework. By simply cutting the mathematical cover (MC) by the discontinuity, the physical cover (PC) will be separated and the discontinuity will, for most cases, be captured without further requirement of incorporating enrichment functions. Due to its capability in dealing with continuous–discontinuous problems, NMM has been successfully extended to simulate cracking problems [26–32] and cracking involved failure problems [33–35]. This study applies the cracking evolution technique developed in our previous work [30,31,35] to analyze the stability of rock blocks in the cavern/tunnel roof based on measured properties of fractures and joints. Using the developed NMM, the entire trajectories of the falling block intersected by natural joints and mining induced fractures under the action of gravity for both cases are reproduced. The excavation procedure is conducted to simulate the effect of in-situ geostress. The breakage of the shotcrete reinforcement is simulated. The influence of the lateral confining pressure on the failure mechanism of the rockfall around the tunnel is investigated. Most importantly, a simple truss rockbolt element is incorporated to investigate the mechanism of rockbolt in mitigating the underground rockfall hazard.

2. Background of NMM 2.1. Introduction The core and most innovative feature of the NMM is the adoption of a two cover (mesh) system, on which the nodes and elements are generated. Similar to other PUMs, the finite covering of a problem domain is also the basic construction of the NMM. The finite cover systems employed in the NMM are referred to the mathematical cover (MC) and physical cover (PC), respectively [24]. The MC, which is used for building PCs, can be either a mesh of regular pattern or a combination of some arbitrary figures. However, the whole mesh has to be large enough to cover the whole physical domain. The physical mesh, which includes the boundary of the material, joints, cracks, blocks and interfaces of material zones, is a unique portrait of the physical domain of a problem, and defines the integration fields. The intersection of the MC and the physical mesh, or the common area of the two systems, defines the region of the PCs. A common area of these overlapped PCs or an independent PC corresponds to an element in the NMM. Fig. 1 illustrates the basic construction procedures of the finite covering system adopted in the NMM. As illustrated in the figure, the MCs first form from the mathematical meshes, such as the rectangular MC M1 and the circle MC M2. From the formed MCs and the physical meshes, the PCs are then defined. For example, MC M1 intersecting with the physical boundary Cu and the discontinuity boundary CD forms the PCs P 11 ; P 21 and P31 , while the MC M2 intersecting with the physical boundary Cu forms PCs P 12 and P22 . Finally, the NMM elements are created by overlapping these PCs, e.g. elements E3 and E4 form from the overlapping of PCs P21 and P 22 ; P 31 and P 22 , respectively. The left independent areas of

PCs then form the other manifold elements, such as element E1 from PC P21 , E2 from PC P31 and E5 from PC P22 . On each PC Pi, a local approximation functions di(x) is independently defined. A convenient way for constructing a basis of local approximation spaces is by using the polynomial functions. e.g.

  c l ¼ 1; x; y; . . . ; xp ; xp1 y; . . . ; xyp1 ; yp

ð1Þ

for a two-dimensional problem, where the superscript ‘‘c’’ stands for conventional PCs. Though the polynomials can approximate smooth functions well and capture the discontinuity directly for the conventional PCs, for PCs that are not fully intersected by the discontinuities (such as PC P 22 in Fig. 1), the smooth basis polynomial local approximations can capture neither the high gradient at the crack tips nor the jumps across the discontinuity surfaces. Therefore, special singular functions are needed to be used to enrich the approximation space for capturing the singularities without refining the meshes. In this study, the linear elastic asymptotic crack-tip fields as proposed by Belytschko et al. [36] are used as suitable enrichment functions for the singular PCs. Then the enriched local approximation functions for singular PCs become

  h pffiffiffi h pffiffiffi h pffiffiffi h s c pffiffiffi l ¼ d ; r sin ; r cos ; r sin h sin ; r sin h cos 2 2 2 2

ð2Þ

where r and h denote polar coordinates in the local system at the crack tip. The superscript ‘‘s’’ stands for singular PCs. Such enrichment not only leads to a better accuracy on relatively coarse meshes pffiffiffi (with term r , the high gradient stress/strain around the crack tip can be captured), but also enables the discontinuities across the cracks that only partially cut MCs be captured (with term sin(h/2), the discontinuity across the discontinuity surfaces for h = ±p can be captured). After constructing the local approximation functions, the local approximation functions can then be connected together by the weighting function (the same as PU function) to form a global displacement function for each manifold element as follows:

uðx; yÞ ¼

n n m X X X j j wi ðx; yÞui ðx; yÞ ¼ wi ðx; yÞ li di i¼1

i

ð3Þ

j

where wi(x) is the weighting function, which satisfies

wi ðx; yÞ  0;

8ðx; yÞ 2 Pi ;

8ðx; yÞ R Pi

wi ðx; yÞ ¼ 0;

with

X wi ðx; yÞ ¼ 1

ð4Þ

x2P i

ui(x, y) is the displacements of PC i, which is defined as:

ui ¼

m X j j li di

ð5Þ

j¼1 j

where li is the basis of the approximation displacement function of j PC i, di is the DOFs related to PC i. Therefore, the displacement for each manifold element can be obtained as:

uðx; yÞ ¼

n n m X X X j j wi ðx; yÞui ðx; yÞ ¼ wi ðx; yÞ li di i¼1

i

ð6Þ

j

The characteristics of the weighting functions depend on the shape of the mathematical mesh and the order of the displacement approximation field. Detailed construction of the weighting function can be referred to the work done by Lin [37]. In this study, a regular triangular mesh and a linear displacement approximation field are adopted. The weighting functions are hence equal to the three node triangular finite element shape functions.

71

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 1. Illustration of the finite cover system in NMM.

2.2. Governing equations Consider a linear elastic dynamic domain X as shown in Fig. 2, the weak form of the governing equation and associated boundary conditions can be expressed as:

Z

rðuÞ : eðduÞdX þ

X

¼

Z Ct

Fig. 2. Notations of a typical problem.

Z

qu€  dudX þ k

X

t  dudC þ

Z

Z

Þ  dudC ðu  u

Cu

b  dudX 8du 2 V 0

ð7Þ

X

where r and e are stress and strain tensor, respectively; b is the  , t are the essential boundary conditions on Cu body force and u and traction boundary condition on Ct, respectively; q is the density of the material; k is the real penalty number which is chosen to superimpose the essential boundary conditions. In the NMM, finite-dimensional subspaces Vh  V and V h0  V 0 are used as the approximating trial and test spaces. Using the Galerkin method [37], the weak form for the discrete problem can be stated as:

72

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Find uh e Vh  V such that

Z

rðuh Þ : eðduh ÞdX þ

X

¼

Z Cþ C

Z

qu€ h  duh dX þ k

X

duh  f þ dC þ

Z

 Þ  duh dC ðuh  u

Cu

Z X

b  duh dX 8duh 2 V h0  V 0

ð8Þ

The trial function uh as well as the test functions duh are expressed as

uh ðxÞ ¼

X wi ðxÞui

duh ðxÞ ¼

i X

ð9Þ

wi ðxÞui

ð10Þ

i

€ n þ bu € nþ1 Dt u_ nþ1 ¼ u_ n þ ½ð1  bÞu

Substituting the trial and test functions from Eqs. (9) and (10) into Eq. (8), and using the arbitrariness of the test functions, the following discrete system is obtained:

€ þ C u_ þ Ku ¼ F Mu

ð11Þ

€ n þ au € nþ1 Dt2 unþ1 ¼ un þ u_ n Dt þ ½ð1=2  aÞu

Z Z M ij ¼ q W Ti W j dX; C ij ¼ cm W Ti W j dX h Xh Z X Z BTi DBj dX þ k ðW i ÞT  ðW j Þd K ij ¼ Xh Cu Z Z Z  dC Fi ¼ W i  tdC þ W i  bdX þ k Wi  u X

h

Cht

ð15Þ

where a and b are two Newmark’s integration parameters, which should satisfy the following relation:

b  0:5; a  0:25ðb þ 0:5Þ2

where

Cht

integration of weak form in most of the cracked elements (such as element E3 and E4 in Fig. 3). For the enriched elements (whose associated MCs are all partially cut by the crack, such as elements E1 and E2 in Fig. 3), the elements of arbitrary shapes are first partitioned into several sub-triangles, as illustrated in Fig. 3. 7-point Gaussian quadrature is then carried out on each sub-triangle and summation is performed to obtain the integration for the whole element. For time integration schemes, Newmark’s integration scheme is adopted. It relates displacement and velocity at time tn+1 to known values at time tn, and acceleration at time tn+1 can then be expressed as:

ð16Þ

Substituting Eq. (15) into Eq. (11) gives

ð12Þ

eu ¼ e K F

ð13Þ

where

ð14Þ

e¼ K

where Wi is the weighting (interpolation) matrix, Bi is the strain displacement matrix, D is the constitutive matrix for an isotropic linear elastic material.

1

a Dt 2

ð17Þ

Mþ

b

aDt

CþK

   1 1 €n þ F~ ¼ F þ u_ n M 1 u 2a aDt      b Dt b €n C þ  1 u_ n þ 2 u 2 a a

2.3. Integration scheme For space integration, in the conventional NMM, the simplex integration method (see Li et al. [38] for more details) is used to analytically evaluate the element matrices. Mapping to convert an arbitrary shape into a regular pattern is not required and the physical boundary can be of an arbitrary shape. However, simplex integration method cannot integrate the non-polynomial terms induced by the enriched functions for singular PCs (Eq. (2)). In order to overcome this limitation, in this study, the simplex-Gaussian integration method is adopted. For the elements without enriched terms, the simplex integration method is used, which can avoid the element subdivision for the numerical

2.4. Rockbolt element Anchors or rockbolts are often used in geotechnical engineering design for stabilizing potential unstable soil or rock masses. To model the reinforcing effect of the anchors or rockbolts, this study introduces the rockbolt element, which has been adopted by Kim et al. [39] for DDA, to the NMM. Consider a rockbolt with an initial length of l, a stiffness of s, and a pretension of Ft, which is pinned at two points (x1, y1) and (x1, y1) within elements i and j, respectively (Fig. 4).where

l¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx1  x2 Þ2 þ ðy1  y2 Þ2

Fig. 3. Partition of enriched elements into sub-triangles by connecting each vertex to the centroid.

ð18Þ

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

73

A further development of Eq. (26) leads to the following:

Y b

¼

T

T



s

T

T



s deðiÞ EeðiÞ EeðiÞ deðiÞ  deðiÞ EeðiÞ GeðjÞ deðjÞ 2l l   

T

T





T

T lx s deðiÞ GeðiÞ GeðiÞ deðiÞ þ F t deðiÞ weðiÞ þ 2l ly  

T

T lx  deðjÞ weðjÞ ly

ð27Þ

where



T EeðiÞ ¼ weðiÞ GeðjÞ ¼ weðjÞ

 lx

ly 

T lx ly

ð28Þ ð29Þ

Minimizing the total potential energy of the rockbolt element, the following matrices can be found.

s ½EeðiÞ ½EeðiÞ T is added to sub-matrix K ii : l

Fig. 4. Illustration of a rockbolt element.

The displacement for each rockbolt end point can be expressed

s ½EeðiÞ ½GeðjÞ T is added to sub-matrix K ij : l

as:

dx1 ¼ u1 ¼ uðx1 ; y1 Þ; dy1 ¼ m1 ¼ mðx1 ; y1 Þ

ð19aÞ

dx2 ¼ u2 ¼ uðx2 ; y2 Þ; dy2 ¼ m2 ¼ mðx2 ; y2 Þ

ð19bÞ

where u1, v1 and u2, v2 are x, y displacements of points i and j respectively. The displacement of the rockbolt element can be obtained as:

1 dl ¼ ½ðx1  x2 Þðu1  u2 Þ þ ðy1  y2 Þðm1  m2 Þ l

ð20Þ

Assuming lx = (x1  x2)/l and ly = (y1  y2)/l, Eq. (20) can then be rewritten as:



dl ¼ ½u1

m1 

  lx ly

 ½u2

m2 

  lx

ð21Þ

ly

Eq. (21) can be expressed in terms of the PC (element node) displacement vector, dl becomes:



T

T dl ¼ deðiÞ weðiÞ

  lx ly



T

T  deðjÞ weðjÞ

  lx ly

ð22Þ

where de(i), we(i) and de(j), we(j) are PCs (element nodes) displacements and weighting functions vectors of element i and j respectively. Assuming the maximum tensile resistant force of the rockbolt to be Ff, the rockbolt force (f) due to the rockbolt displacement can be expressed as:

dl f ¼ s  F t ; f  F f l f ¼ 0; f > F f

ð23Þ ð24Þ

where the compressive force is positive. The total potential energy of the rockbolt element is:

Y b

1 s 2 ¼  fdl ¼ dl þ F t dl 2 2l

ð25Þ

Substituting Eq. (23) into Eq. (25) gives

Y b

¼

     2

T

T lx

T

T lx s  deðjÞ weðjÞ deðiÞ weðiÞ 2l ly ly     

T

T lx

T

T lx  weðjÞ weðjÞ þ F t deðiÞ weðiÞ ly ly

ð26Þ

s ½GeðjÞ ½EeðiÞ T is added to sub-matrixK ji : l s ½GeðjÞ ½GeðjÞ T is added to sub-matrix K jj : l F t ½EeðiÞ  is added to sub-matrix ½DF i : F t ½GeðjÞ  is added to sub-matrix ½DF j : 2.5. Verification of the rockbolt element To verify the proposed rockbolt element in the NMM, an example as shown in Fig. 5, with 1238 MCs, 1238PCs and 2274 manifold elements, is presented. In this example, a block of known weight is hung from a fixed block using three rockbolts. The hanging block is 5 m long and 2 m wide. It has a density of 2.5  103 kg/m3, and hence a weight of 245 kN. The top block is of the same dimensions but fixed at several points. All the rockbolts have the same stiffness and length (3 m). The rockbolt force for each rockbolt is calculated to be 81.667 kN. It agrees well with the weight of the hanging block (81.677 kN  3 = 245 kN). 2.6. Contact mechanics with cohesion The success of the discontinuous numerical analysis depends heavily on the accurate treatment of contacts between blocks. Wrong detection or description may lead to wrong analysis, or breakdown of the computation. The NMM adopts a rigorous contact searching method as used in DDA to account for the interactions between loops. A loop is an independent closed domain, which can be either a block (an independent closed physical area) or a joint (crack). As shown in Fig. 6, loop one is defined as the closed domain encompassed by the outside boundary vertices, while loop two is defined as the closed domain encompassed by vertices from vertex 1 to vertex 8 and vertex 9 to vertex 16. When the distance of either two vertices or the vertex with the edge of the loops is smaller than a pre-defined value d, they are treated as a possible contact pair and the contact status will be judged during the iteration of computation.

74

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 5. Verification example for the rockbolt elements.

Fig. 6. Representation of loops and loop vertices in the NMM.

The two-dimensional NNM adopts the following three types of contacts: angle-to-angle, angle-to-edge, edge-to-edge (Fig. 7). An edge-to-edge contact can be transformed into two angle-to-edge contacts, i.e. angle P2 to edge P1P3 and angle P3 to edge P2P4. Based on the contact states, the interactive normal and shear forces Fn and Fs can be obtained from the following criteria: For the ‘‘open’’ status (w1 > 0), the normal and shear forces Fn and Fs are:

Fn ¼ Fs ¼ 0

ð30Þ

and for the ‘‘closed’’ status (w1 6 0), the normal force Fn is:

F n ¼ kn w1

ð31Þ

The tangential force Fs under the ‘‘closed’’ status is dependent on the Coulomb friction law: For the ‘‘sticking’’ state, i.e. Fs 6 Fn  tan (/j) + cj  lc,

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

75

Fig. 7. Three types of contacts (a) angle-to-angle, (b) angle-to-edge and (c) edge-to-edge.

F s ¼ ks w2

ð32Þ

For the ‘‘sliding’’ state, i.e. Fs > Fn  tan (/j) + cj lc,

F s ¼ F n  tanð/j Þ

analysis which have been successfully developed in our previous work [30,31] are described briefly.

ð33Þ

where kn and ks are the normal and tangential contact spring stiffness; w1 and w2 are the relative normal and tangential displacements of the contact pair; cj and /j are the cohesion and the friction angle of the discontinuity, respectively; and lc is the calculated discontinuity length of the contact. 3. Extended NMM for rock failure analysis It is well known that failure processes of rock structures involve the nucleation or activation of pre-existing cracks within the rock matrix and their coalescence which leads to the creation of critical fractures. It is thus essential for a model to reproduce the crack initiation, propagation and coalescence processes which eventually transect the material by a through-going discontinuity. In this section, the basic procedures of extending the NMM to rock failure

3.1. Crack initiation criterion Our prior work has been able to model the crack initiation or propagation in a rock as well as to determine the crack growth direction. Based on the fracture mechanics, several popular crack propagation criteria such as the maximum circumferential tensile stress (MTS) [40], the maximum energy release rate [41], the minimum strain energy release rate [42] and the principle of local symmetry [43] have been developed. The crack propagation of these criteria is determined by the degree of stress concentration around the crack tip, namely the stress intensity factor or strain energy, which can only handle crack initiation and propagation emanating from pre-existing cracks. In the physical tests performed by Wong and Einstein [44] on Carrara marble and gypsum specimens containing pre-existing flaws, however, some of the first cracks did not always emanate from the flaw crack tips. To

Fig. 8. Illustration of crack initiation for elements without crack tip (a) tensile crack and (b) shear crack (O represents the center of the element) [30,31].

Fig. 9. Illustration of crack propagation in elements containing one or more crack tips (a) crack propagation in elements containing single crack tips, (b) newly formed crack in element containing two crack tips and (c) linking the crack tips which have a minimum included angle h1 with the originally predicted crack trajectory in elements containing more than two crack tips [31].

76

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 10. Updating of loops after crack initiation and propagation: (a) Type I: a new loop forms, (b) Type II: a crack is added to the pre-existing loop, (c) Type III: the pre-existing loop is cut into two new loops and (d) Type IV: the pre-existing two loops combine into one [30,36].

Fig. 11. Failure process in a four-layered rock beam without rockbolts predicted by the NNM (a) Initial configuration, (b) cracking of the first layer, (c) cracking of the second layer, (d) cracking of the third layer and (e) final status.

77

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 12. Failure process in a four-layered rock beam with rockbolts predicted by the NMM (a) initial configuration, (b) cracking of the first layer, (c) cracking of the second layer, (d) cracking of the third layer and (e) final status.

Fig. 13. Failure test results of a four-layered rock beam (a) beam without rockbolt and (b) beam with rockbolts [47].

model the development of a new crack in the intact flawless rock material, an additional criterion is needed. In this study, the Mohr–Coulomb with a tensile cut-off [45] criterion is adopted for the intact elements [30]. New crack initiation occurs in the intact elements when the local stresses reach the tensile or shear strength of the material, while the MTS criterion is used to consider the effect of pre-existing flaw on the crack propagation. The following section provides a detailed description of the Mohr–Coulomb crack initiation criterion and the MTS crack growth criterion.

To model a new crack initiation in the intact material, a tensile or shear crack initiates when one of the following two conditions is satisfied. For shear crack,

r1  r3 2

¼ ci cos /i þ

r1 þ r3 2

sin /i ;

r3 > rt

ð34Þ

For tensile crack,

r3 ¼ rt

ð35Þ

78

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 14. A very shallow tunnel failure due to excavation (rh = 0 MPa) (a) initial geometry, (b) computational step 4 and (c) computational step 5000.

where ci and /i are the cohesion and friction angle of the intact material; r1 and r3 are the maximum principal stress and the minimum principal stress, respectively. For the crack propagation cases, the pre-existing crack will propagate when the following condition is satisfied.

cos

  h0 h0 3 K I cos2  K II sin h0 ¼ K IC 2 2 2

ð36Þ

where KI and KII are the mode I and mode II stress intensity factors respectively. KIC is the mode I fracture roughness and the propagation angle h0 can be obtained by the following equation:

 8 pffiffiffiffiffiffiffiffiffiffiffiffi < 2 arctan jK I j K 2I þ8K 2II K II ¼ 0 4K II h0 ¼ : 0 K II – 0

ð37Þ

3.2. Treatment of manifold elements and relevant covers during cracking process For the crack initiation case, when the element without a crack tip satisfies the failure criterion Eq. (34) or Eq. (35), a crack will initiate. We further assume that the newly initiated crack passes

through the center of the element. Then if the failure criterion Eq. (35) is met, a tensile crack perpendicular to r3 as shown in Fig. 8(a) will develop. If the principal stresses of the elements satisfy the condition expressed in Eq. (34), a shear crack at an angle of (p/4 + /i/2)with the direction of r3 as shown in Fig. 8(b) will develop [30,31]. For crack propagation cases, complications arise when an element contains more than one crack tip satisfying the failure criterion in Eq. (36). Three different cases are described below. If an element contains only one crack tip that reaches the failure criterion, the crack will propagate in a direction predicted by Eq. (37), as shown in Fig. 9(a). If an element, that contains two crack tips, reaches the failure criterion, the newly formed crack is assumed to link the original two cracks together, as illustrated in Fig. 9(b). If an element, that contains more than two crack tips, reaches the failure criterion, the newly formed crack will link up the crack tip that has a minimum included angle h1 with the crack propagation direction predicted according to Eq. (29) (Fig. 9c). After crack initiation and propagation, the newly formed cracks then cut the MCs into certain sub-domains. The PCs and the manifold elements will be subsequently updated automatically according to the procedure discussed in Section 2.1.

Table 1 Material properties used in the present study. Intact rock

Joint

Density q (kg/m3)

Young’s modulus E (GPa)

Poisson’s ratio v

Tensile strength rt (MPa)

Friction angle /i (°)

Cohesion ci (MPa)

Friction angle /j (°)

Cohesion cj (MPa)

Fracture toughness KIC (MPa m1/2)

3100

30.0

0.15

13.1

23.5

18.1

17

0.1

1.0

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

79

Fig. 15. A very shallow tunnel with rockbolts reinforcement (rh = 0 MPa) (a) initial geometry, (b) computational step 5000 for case g = 1.0g0 and (c) computational step 5000 for case g = 1.82g0. (g0 is the acceleration of gravity in the initial state).

3.3. Treatment of loops after crack initiation and propagation In the NMM, cracks are treated as loops which are the bases of the contact techniques. Therefore, after the crack initiation or propagation, the loops need to be updated. Based on the relative position between the newly developed crack and the original loops, four types of loop updating as shown in Fig. 10 can be defined. Interested readers can refer to literatures [30,31,34] for further details. 3.4. Calibrations Layered geological structures are commonly encountered in rock engineering projects. The failure of a layered rock structure is commonly attributed to the weak interfaces between the rock layers. Rockbolts are usually used to maintain the integrity of the layered structure as a whole. In this section, an example of a four-layered rock beam is presented to demonstrate the effectiveness of the rockbolt elements and the cracking evolution techniques described in Section 2. In addition, the mechanisms of the rockbolt in reinforcing the layered structure are modeled and discussed. To investigate the effect of the rockbolt, a four-layered rock beam with or without rockbolts is simulated. The length of the beam is 448 mm and the thickness of each layer is 32 mm, as shown in Figs. 11a and 12a. The beam is simply supported and under two initial applied forces F = 4 kN. The applied force F increases by 0.1 kN after each step. The material used is man-made gypsum and its properties are 10 GPa for Young’s modulus, 0.25 for

Poisson’s ratio, 0.5 MPa for tensile strength and 0.03 kN cm1/2 for the fracture toughness. Each rockbolt is assumed to be an elastic bar. The diameter and the length of the elastic rockbolt used are 3.0 mm and 127.0 mm, respectively, while its Young’s modulus is 20 GPa. The contact surfaces between each beam layer are assumed to be smooth, without friction and cohesion. As shown in Figs. 11 and 12, for both beams, when the applied forces F increases, the bottom layer cracks first. When the load keeps on increasing, the top layers also crack. However, for the beam without any rockbolt reinforcement, the bottom two layers collapse when the applied forces F increases to 14.4 kN; while for the beam with rockbolt reinforcement, no collapse occurs after a similar deformation has attained. This is because with the rockbolt reinforcement, the slip between the layered surfaces can be better controlled. As a result, the layered beam with rockbolts behaves as an integrated unit, hence achieving a higher overall strength. Compared with the physical results obtained by Yang and Xia [46], for both beams, not only the failure pattern but also the strength can be satisfactorily predicted by the NMM (Fig. 13).

4. Tunnel stability analysis Rockfall from the tunnel/cavern roof is highly undesirable during the underground excavation. The potential unstable block can be defined by natural joints or mining induced fractures. The fallen blocks may be fatal and possibly lead to a more massive collapse subsequently. Therefore, to identify these critical blocks during excavation is very important. Since rockbolts are commonly used

80

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 16. A tunnel failure due to excavation (rh = 0 MPa) (a) initial geometry, (b) computational step 20 and (c) computational step 2000.

to reinforce these unstable blocks, it is highly beneficial to investigate the effect of the rockbolt reinforcement. In this section, the rockfall from the tunnel/cavern roof due to excavation under different conditions is investigated. The effect of the release of in situ geostress due to the excavation is simulated by updating the physical covers to delete the elements in the tunnel after the initial equilibrium is reached. 4.1. Rockfall due to mining induced fracture initiation A very shallowly-buried horse-shoe-shaped tunnel located in an intact rock mass, as shown in Fig. 14a, is simulated. The initial NMM model contains 997 MCs, 1033 PCs and 1870 Manifold elements. The material properties are listed in Table 1. No confining stress is considered in this simulation. A time step size of 0.0004 s is selected. As illustrated in Fig. 14b, immediately after the excavation, fractures initiate at the roof of the tunnel due to the stress release at the free surfaces. Without reinforcement, the initiated fractures are unstable and keep propagating until the unstable rock blocks form finally (Fig. 14c). To reinforce this shallowly-buried tunnel, a reinforcing scheme consisting of 4 m long rockbolts at a spacing of 4 m as shown in Fig. 15a is adopted. As illustrated in Fig. 15b, fractures still initiate when the rock mass inside the tunnel is excavated. However, due to the rockbolt reinforcement, the

displacement of the rock mass around the free surfaces is limited, which prohibits the initiated fractures from propagating. Therefore, with the rockbolt reinforcement, the very shallow tunnel can be stable in spite of the fractures induced by the excavation associated with the stress release. To investigate the failure mode of the shallow tunnel after rockbolt reinforcement, a loading scheme involving a gradual increase of the gravity acceleration g until the formation of rockfall blocks is implemented. As illustrated in Fig. 15c, with the rockbolt reinforcement, the final failure of the shallow tunnel is different from the one without reinforcement. With the rockbolt reinforcement, the roof of the shallow tunnel maintains a high integrity and the tensile failure (fracture) is mostly prohibited. With the increase of the gravity acceleration g (g = 1.82g0), two inclined shear bands trending away from the rockbolt reinforcing region form. The shallow tunnel thus fails in a shear sliding mode. 4.2. Rockfall due to mining induced fracture propagation In this section, the same horse-shoe-shaped tunnel as discussed in Section 4.1 but with different boundary conditions and the overburden thickness located in a rock mass containing two major preexisting joints (Fig. 16a) is simulated. The initial NMM model contains 763 MCs, 795 PCs and 1431 Manifold elements. The material properties are listed in Table 1. A dynamic control parameter of 1.0

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

81

numerical simulation findings of a rectangular rock specimen containing a single pre-existing flaw under biaxial compression [30]. As such, no independent block can form to lead to rockfall events. To avoid rockfall during tunnel excavation under a low lateral confining stress (rh 6 2 MPa), rockbolts are usually used for reinforcing the unstable rock mass. Since only one independent block is induced for the scenario shown in Fig. 16a, it is easy to determine where the key block is. To prevent the key block from falling, two elastic rockbolts as shown in Fig. 17b are applied. Each rockbolt is 5 m in length and 2.5 cm in diameter. The Young’s modulus of each rockbolt is chosen to be 200 GPa and a pre-stress 20 kN is applied. Simulation results shown in Fig. 17b illustrate their positive reinforcing effect. After the excavation and the crack evolution at both pre-existing joints, most of the stresses around the roof of the tunnel are released. Due to the applied rockbolts, though an independent block still forms during the excavation, it is tightly fastened to the surrounding rock mass and no rackfall can happen. The bolt force in both rockbolts is around 120 kN. 4.3. Rockfall due to cutting of natural joints

Fig. 17. Final state at different scenarios (a) final principal stress vectors for case rh = 5 MPa and (b) final principal stress vector for case rh = 0 MPa with rockbolt reinforcement.

is chosen, which means the initial velocity of each element in each time step inherits 100% of the element final velocity of the previous time step. A time step size of 0.0004 s is selected. As for the case without the lateral confining stress (rh = 0 MPa), due to the immediate re-distribution of stresses in surrounding rock mass during tunnel excavation, cracks initiate at both preexisting joints. At t = 0.008s after the excavation, the crack propagated from joint 1 coalesces with the pre-existing joint 2, which eventually transects the rock mass as a through-going joint to produce an independent block (Fig. 16b). The induced block is then driven by the gravity force to detach from the roof of the tunnel to fall onto the floor. Since the weight of the block is around 535 kN and the height of the tunnel is around 29 m, its failure is therefore fatal to any workers beneath. To investigate the effect of the lateral confining stress on the rock stability especially the prohibition of the pre-existing fracture from propagation, the following lateral confining stresses are applied and studied: rh = 0.5 MPa, 1 MPa, 2 MPa, 5 MPa, 10 MPa. Simulation results illustrate that when the confining stress rh 6 2 MPa, a failure mode similar to that without confining stress is predicted. For the higher lateral confining stress rh P 5 MPa, which leads to the arch effect (Fig. 17a), the crack emanating from the pre-existing joints is prohibited. This agrees with our previous

The scenario discussed in Section 4.2 is relatively simple, in which the keyblock is easy to be identified. However, the actual underground engineering condition is usually more complicated due to the presence of multiple potential unstable blocks. Therefore, it is always challenging to precisely determine the location of keyblock(s). In this section, through a more complicated case, the capability of the developed NMM in locating the keyblock and providing insights on suitable reinforcement is illustrated. Fig. 18 shows a typical cross section of the underground cavern housing a hydropower station [47]. The geometry of the power house (cavern) is the same as the tunnel illustrated in Fig. 16a. The only difference is the presence of multiple natural joints. Due to the cutting of the natural joints, a number of independent blocks are formed. Some of them are potentially unstable blocks (Blocks 1 to 5 as shown in Fig. 18). However, before excavation, it is really hard to assess the stability of these blocks upon excavation. It is even more challenging to determine which block is the keyblock. To investigate the stability and locate the keyblocks of the underground power house, an NMM model with 824 MCs, 1022 PCs, 1648 manifold elements and 15 loops is developed. The material properties used are the same as those listed in Table 1 except for the fracture toughness of joints. For simplicity, the crack propagation which is not the key factor in this scenario is not considered. Therefore, the joint fracture toughness KIC is set to a value of 100 MPa/m1/2. A full dynamic analysis with a time step size of 0.0004 s is performed. The failure processes predicted by the NMM for the case without lateral confining stress (rh = 0 MPa) and no support are illustrated in Fig. 19. As shown in the figure, immediately after the rock mass in the power station house is removed to release the constraint in radial direction, block 5 begins to detach from its upper part (Fig. 19a). When block 5 is falling down, since the low confining stress cannot resist the gravity load, the whole roof (blocks 1 to 4) begins to detach from the upper rock mass and falls down (Fig. 19b). As a result, the free surface of the second layer blocks is exposed, which leads to their detachment from the top of the rock mass (Fig. 19c). In a similar way, the third layer blocks subsequently detach and fall down (Fig. 19d). According to Bray [48] and Sofianos et al. [49], the factor of safety for rock blocks is defined as:

FS ¼

P0 jWj

ð38Þ

82

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 18. Geometry of one cross section of the power station.

Fig. 19. Failure processes of the power station predicted by the NMM under no lateral confining stress (rh = 0 MPa) and without support (a) computational step 750, (b) computational step 1800, (c) computational step 2500 and (d) computational step 5000.

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

83

Fig. 20. Failure state of the power station predicted by the NMM under a very high lateral confining stress (rh = 10 MPa) without support. (a) Principal stress vectors at computational step 3000 and (b) failure state at computational step 3000.

Fig. 21. Failure processes of the power station predicted by the NMM under no lateral confining stress (rh = 0 MPa) with 20 cm thick shotcrete support (a) computational step 300 and (b) computational step 5000.

Fig. 22. Failure state of the power station predicted by the NMM under a very high lateral confining stress (rh = 10 MPa) with a 20 cm thick shotcrete support. (a) Principal stress vectors at computational step 3000 and (b) failure state at computational step 3000.

where P0 is the pullout resistance of the rock block, which can be evaluated as follows:

P0 ¼ 2MH0

ð39Þ

in which M is a function of the mechanical properties of the fractures and the rock block inclination angle, and H0 is the horizontal confining force applied to the rock block by the surrounding rock mass. As illustrated in Eqs. (38) and (39), the horizontal confining stress can significantly increase the stability of rock blocks against falling. For deeply buried-underground structures, very high horizontal stresses are usually present in the rock mass. To investigate the positive effect of such horizontal confining stress on the stability of the station house, the following lateral confining stresses are

applied: rh = 1 MPa, 2 MPa, 5 MPa, 10 MPa. Simulation results illustrate that when the confining stress is low (rh 6 2 MPa), a failure mode similar to that of no confining stress is observed. Under a high confining stress (rh P 5 MPa), all the blocks at the roof are stable due to the arch effect (Fig. 20a) except for block 5. Even under a very high confining stress (rh = 10 MPa), block 5 is still unstable (Fig. 20b). The above results illustrate that the power house is unstable without support regardless of what level the confining stress is. Therefore, additional support is needed. To investigate the feasibility of different reinforcement schemes, in the following, the effect of a 20 cm thick shotcrete support under different confining stresses is first investigated. Secondly, the combined effect of the 20 cm thick shotcrete support and the rockbolt reinforcement is studied.

84

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85

Fig. 23. Final state of the power station predicted by the NMM under no lateral confining stress (rh = 0 MPa) in scheme 2 with combined support. (a) Final state and (b) final principal stress vectors.

The material properties of the shotcrete support are 24 kN m3 for unit weight, 26 GPa for Young’s modulus, 0.167 for Poisson’s ratio, 1.8 MPa for cohesion, 60o for friction angle and 1.7 MPa for tensile strength. In the simulation of the stability of the power house with a 20 cm thick shotcrete support, two conditions without and with a high confining stress are considered. For the first case without confining stress, due to the relatively large gravitational pulling force, the shotcrete fails immediately after the excavation has finished (Fig. 21a). After the breakage of the shotcrete, the remaining structure behaves similar to the case without support. As a result, a failure pattern similar to that without shotcrete support is observed (Fig. 21b). For the second case with a very high confining stress (rh = 10 MPa), although the associated arch effect can keep most parts of the roof stable (Fig. 22a), the weak support still cannot resist the downward falling force caused by block 5. After the shotcret breaks, block 5 detaches from the roof and falls onto the floor (Fig. 22b). These two results clearly illustrate the requirement of a stronger reinforcement in order to keep the power house stable. Previous results illustrate that block 5, which fails in all of the analyzed scenarios, is always the first one to fall down. Therefore, block 5 is determined as the keyblock, which must be secured by suitable support. Since rockbolts have been proven to be a useful type of reinforcement, a combined reinforcement scheme consisting of rockbolts and shotcrete is further considered below. Three rockbolts and a 20 cm thick shotcrete support are installed immediately after the excavation (Fig. 23a). The lengths of the rockbolts are 3 m, 6 m and 9 m for bolts 1, 2 and 3 respectively. The Young’s modulus of each rockbolt is chosen to be 200 GPa and a pre-stress of 30 kN is applied. All the rockbolts are inclined at 45° to the horizontal plane. To investigate the stiffness effect of the rockbolts on the reinforcement performance, two schemes adopting different rockbolt diameters are simulated. In scheme 1, the diameter of each rockbolt is 2.5 cm. In scheme 2, the diameter of rockbolt 1 is 2.5 cm and the diameter of rockbolts 2 and 3 are both 3 cm. To be conservative, the stability of the power house under the worst condition, i.e. without confining stress (rh = 0 MPa) is investigated. Simulation results illustrate the excellent reinforcing performance for both schemes. The power house can be kept stable with a maximum deformation of around 14 mm for scheme 1 and around 10 mm for scheme 2. Despite the similarly good reinforcing results achieved in both schemes, the final stresses in the rockbolts of these two schemes are different. In scheme 1, the final stresses in rockbolts 1, 2 and 3 are 345.48 MPa, 457.99 MPa and 302.68 MPa respectively. In scheme 2, the final stresses in rockbolts 1, 2, and 3 are 332.02 MPa, 341.18 MPa and 272.68 MPa respectively. The rockbolt stresses of

scheme 2 are more uniformly distributed and have a much lower overall maximum stress than those in scheme 1. Therefore, scheme 2 is a more preferable choice for reinforcing the power house. 5. Conclusions As demonstrated in this study, the developed NMM appears to be a useful tool for modeling the rockfall hazard in underground engineering. By combining the widely used FEM and DDA in a uniform framework, not only the deformation, but also the block detachment and movement can be accurately modeled. The entire processes of rockfall have been captured dynamically and valuable information including falling time, pattern, weight and velocity have been obtained. Simulation results illustrate that the keyblock which is very critical in rockfall assessment can be located by the developed NNM. The lateral confining stress in the underground structure not only helps stabilize the dissected blocks but also prohibits the pre-existing fractures from propagating to form new unstable blocks. The simulation also illustrates that rockbolts are usually very effective in reinforcing underground rock structures, especially the layered structures. To summarize, the NMM developed shows its promising capability as a tool to help locate possible rockfalls and assist the design of the corresponding reinforcement support. However, the NMM is still at its infant stage of development, a further development such as incorporating more complicated rockbolt element to consider the breakage and sliding effect between the rockblot and rock mass is needed to target at practical rockfall assessments. Acknowledgements The research was supported by the Academic Research Fund Tier 1 (RG19/10) from the Ministry of Education, Singapore and the Nanyang Technological University Start Up Grant (M4080115.030). The authors would also like to thank Dr. Gen-Hua Shi, Dr. Youjun Ning and Professor Guowei Ma for valuable discussion, which inspired the authors to conduct the present research. References [1] Jaboyedoff M, Dudt JP, Labiouse V. An attempt to refine rockfall hazard zoning based on the kinetic energy, frequency and fragmentation degree. Nat Hazards Earth Syst Sci 2005;5:621–32. [2] Alejano LR, Stockhausen HW, Alonso E, Bastante FG, Oyanguren PR. ROFRAQ: A statistics-based empirical method for assessing accident risk from rockfalls in quarries. Int J Rock Mech Min Sci 2008;45:1252–72. [3] Pierson L, Davis S, van Vickle R. The rockfall hazard rating system implementation manual. Technical report, Oregon State Highway Division, 1993.

Z. Wu, L.N.Y. Wong / Advances in Engineering Software 76 (2014) 69–85 [4] Crosta GB, Agliardi F. A methodology for physically based rockfall hazard assessment. Nat Hazards Earth Syst Sci 2003;3:407–22. [5] Peila D, Patrucco M, Falanesca M. Quantification and management of rockfall risk in opencast quarrying activities. Environ Eng Geosci 2011;17:39–51. [6] Goodman RE. Methods of geological engineering in discontinuous rocks San Francisco. CA: West Publishing Company; 1976. [7] Babuska I, Melenk JM. The partition of unity method. Int J Numer Meth Eng 1997;40:727–58. [8] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [9] Moes N, Dolbow J, BELytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. [10] Cundall PA, Strack ODL. Discrete numerical-model for granular assemblies. Geotechnique 1979;29:47–65. [11] Shi GH, Goodman RE. 2 dimensional discontinuous deformatio analysis. Int J Numer Anal Meth Geomech 1985;9:541–56. [12] Choi SO, Chung SK. Stability analysis of jointed rock slopes using the Barton– Bandis constitutive model in UDEC. Int J Rock Mech Min Sci 2004;41:469. [13] Alejano LR, Ferrero AM, Ramirez-Oyanguren P, Fernandez MIA. Comparison of limit-equilibrium, numerical and physical models of wall slope stability. Int J Rock Mech Min Sci 2011;48:16–26. [14] 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:671–94. [15] 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:483–90. [16] Xia CC, Xu CB, Zhao X. Study of the strength reduction DDA method and its application to mountain tunnel. Int J Comput Methods 2012. http://dx.doi.org/ 10.1142/S021987621250041. [17] Salciarini D, Tamagnini C, Conversini P. Numerical approaches for rockfall analysis: a comparison. In: Anderssen RS, Braddock RD, Newham LTH, editors. 18th World IMACS/MODSIM Congress Cairns, Australia, 2009. p. 2706–12. [18] Itasca. PFC 2D (Partical Flow Code in 2 Dimension). Version 1.1. Itasca Consulting Group, Minneapolis, 1999. [19] Itasca. PFC 3D (Partical Flow Code in 3 Dimension). Version 1.1. Itasca Consulting Group, Minneapolis, 1999. [20] Donze FV, Bouchez J, Magnier SA. Modeling fracturing in rock blasting. Int J Rock Mech Min Sci 1997;34(8):1153–63. [21] Camones LAM, Vargas Jr EdA, de Figueiredo RP, Velloso RQ. Application of the discrete element method for modeling of rock crack propagation and coalescence in the step-path failure mechanism. Eng Geol 2013;153:80–94. [22] Zhang XP, Wong LNY. Cracking processes in rock-like material containing a single flaw under uniaxial compression: a numerical study based on parallel bonded-particle model approach. Rock Mech Rock Eng 2012;45:711–37. [23] Donze FV, Richefeu V, Magnier SA. Advances in disctrete element method applied to soil rock and concrete mechanics. Electron J Geotech Eng 2009;44:1–44. [24] Shi GH. Manifold method of material analysis. Trans 9th Army Conf on Applied Mathematics and Computing Minneaplolis: Minesota, 1991. p. 57–6. [25] Shi GH. Modeling rock joints and blocks by manifold method. In: Proceedings of the 33th US Rock Mechanics Symposium New Mexico: SanTa Fe, 1992. p. 639–48.

85

[26] Tsay RJ, Chiou YJ, Chuang WL. Crack growth prediction by manifold method. J Eng Mech–ASCE 1999;125:884–90. [27] Chiou YJ, Lee YM, Tsay RJ. Mixed mode fracture propagation by manifold method. Int J Fract 2002;114:327–47. [28] Li SC, Cheng YM. Enriched meshless manifold method for two-dimensional crack modeling. Theoret Appl Fract Mech 2005;44:234–48. [29] Ma GW, An XM, Zhang HH, Li LX. Modeling complex crack problems using the numerical manifold method. Int J Fract 2009;156:21–35. [30] Wu ZJ, Wong LNY. Frictonal crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39:38–53. [31] Wu ZJ, Wong LNY. Elastic–plastic cracking analysis for brittle-ductile rocks using manifold method. Int J Fract 2013;180:71–91. [32] Zhang HH, Li LX. Modeling inclusion problems in viscoelastic materials with the extended finite element method. Finite Elem Anal Des 2009;45:721–9. [33] Zhang GX, Zhao Y, Peng XC. Simulation of toppling failure of rock slope by numerical manifold method. Int J Comput Methods 2010;7:167–89. [34] Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. Int J Rock Mech Min Sci 2011;48:964–75. [35] Wu ZJ, Wong LNY. Dynamic study on fracture problems in viscoelastic sedimentary rocks using the numerical manifold method. Rock Mech Rock Eng 2013;46:1415–27. [36] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Meth Eng 1994;37:229–56. [37] Lin JS. A mesh-based partition of unity method for discontinuity modeling. Comput Methods Appl Mech Eng 2003;192:1515–32. [38] Li S, Cheng Y, Wu YF. Numerical manifold method based on the method of weighted residues. Comput Mech 2005;35:470–80. [39] 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:949–70. [40] Erdogan F, Sih SC. On the crack extension in plates under plane loading and transverse shear. J Basic Eng ASME 1963;85:519–25. [41] Hussian MA, Pu SL, Underwood J. Strain energy release rate for a crack under combined mode I and mode II. Fract Anal ASTM STP 1974;560:2–28. [42] Sih GC. Strain-energy-density factor applied to mixed crack problems. Int J Fract 1974;10:305–21. [43] Cotterell B, Rice JR. Slightly curved or kinked cracks. Int J Fract 1980;16:155–69. [44] Wong LNY, Einstein HH. Crack coalescence in molded gypsum and carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech Rock Eng 2009;42:475–511. [45] Brady BHG, Brown ET. Rock mechanics for underground mining. 2nd ed. London, UK: Chapama&Hall; 1993. p. 571. [46] Yang JH, Xia JZ. Test study on the complete deformation property of bolted stratified surrounding rockmasses. J China Coal Soc 2005;30:414–7 [In Chinese]. [47] ReS Luis, Eurípedes VJ, Fernandes MM, Azevedo R. Innovative numerical modelling in geomechanics. CRC Press; 2012. [48] Bray JW. Unpublished note. London: Imperial College, 1977. [49] Sofianos AI, Nomikos P, Tsoutrelis CE. Stability of symmetric wedge formed in the roof of a circular tunnel: nonhydrostatic natural stress field. Int J Rock Mech Min Sci 1999;36:687–91.