J. Cent. South Univ. Technol. (2008) 15: 888-894

DOI: 10.1007/s11771-008-0162-0

Simulation of failure process of jointed rock

ZHANG Xiu-li(张秀丽)¹, JIAO Yu-yong(焦玉勇)¹, ZHAO Jian(赵 坚)²

 (1. State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China;

2. Ecole Polytechnique Federale de Lausanne, Rock Mechanics Laboratory, CH-1015 Lausanne, Switzerland)

Abstract: A modified discontinuous deformation analysis (DDA) algorithm was proposed to simulate the failure behavior of jointed rock. In the proposed algorithm, by using the Monte-Carlo technique, random joint network was generated in the domain of interest. Based on the joint network, the triangular DDA block system was automatically generated by adopting the advanced front method. In the process of generating blocks, numerous artificial joints came into being, and once the stress states at some artificial joints satisfy the failure criterion given beforehand, artificial joints will turn into real joints. In this way, the whole fragmentation process of rock mass can be replicated. The algorithm logic was described in detail, and several numerical examples were carried out to obtain some insight into the failure behavior of rock mass containing random joints. From the numerical results, it can be found that the crack initiates from the crack tip, the growth direction of the crack depends upon the loading and constraint conditions, and the proposed method can reproduce some complicated phenomena in the whole process of rock failure.

Key words: discontinuous deformation analysis; jointed rock; rock failure; Monte-Carlo technique; random joint network; advancing front method; triangular block system

1 Introduction

Natural rock mass exists in compression stress environment, and its failure is related to the initiation, propagation and coalescence of interior cracks. In the past decades, significant researches on the topic of crack growth have been done experimentally^[1-4], analytically^[5-9] and numerically^[10-18]. However, due to the mechanical and geometrical complexity of this problem, it has not been yet solved fully and perfectly.

Among all the available means, numerical approach is proven to be an effective tool for simulating rock failure because the loading and boundary conditions of the model can be controlled strictly, and some parameters, such as crack length and crack location can be easily changed so as to study their influences on the mode of crack growth. So far, by using some numerical methods embedded with the fracture theory, a number of attempts have been made to study the mechanism of crack growth^[12-18]. CHARALAMBIDES and MCMEEKING^[13] adopted a finite element method to simulate crack propagation for studying the effects of microcracking at facets as a damage process in brittle composites and the consequence for material toughness. In their study, crack propagation was achieved by successive tip node relaxation upon satisfaction of a critical energy release rate criterion. TAN et al^[14] used a displacement discontinuity method (DDM) program coupled with a modified energy criterion to simulate the development of cracks, and the numerical results showed that the crack may either propagate to the free surface or rest in the rock subsurface, which is dependent on the rock and fracture properties, loading force. XU and SAIGAL^[15] established the discrete formulation for stable crack growth in an elastic solid using the element free Galerkin (EFG) method based on the moving least-squares approximations. CHEN and LI^[16] presented a 3D distinct element method for block deformation and rupture analysis. However, few efforts on the propagation and coalescence of multi-cracks have succeeded, and no results on the failure process of rock with stochastic cracks have been reported. The main reason is that the commonly-used numerical methods, such as FEM, BEM or meshfree method, always encounter mathematical difficulties in dealing with multi-crack problem, while discrete element method requires further efforts on crack growth simulation. In this work, on the basis of the discontinuous deformation analysis (DDA) method proposed by SHI^[19], an approach was presented to simulate the propagation and coalescence of multi-cracks.

In the proposed method, some modifications were made for enabling DDA to handle semi-penetrative joints and to perform rock failure analysis. Firstly, the Monte-Carlo technique was employed to generate random joints whose geometrical parameters satisfied some stochastic distribution laws obtained from geology survey. Then, on the generated joint network, a popular FEM mesh generation method, namely the advancing front method, was applied to automatically generating triangular blocks. In the block generation process, numerous artificial joints came into being. Different from the real joints, the strength of the artificial joints was much higher, and the strength of that intact rock was employed. Through the rupture analysis of the artificial joints, the rock failure process can be simulated in a straightforward way.

2 DDA method

DDA is a numerical approach based on block theory and minimum energy principle. The movement of each block is described by Newton’s second law of motion. And, the penetration of one block into the other is not allowed. The simultaneous equilibrium equations are established by minimizing the total potential energy of the block system due to different mechanisms.

2.1 Function of displacement

In DDA, each block of arbitrary geometry has six degrees of freedom, among which three components are rigid body motion terms and the other three are constant strain terms. Therefore, the deformation variable of block i can be written as

D_i=(u₀  v₀  r₀  ε_x  ε_y  γ_xy)^T                   (1)

where  D_i is the displacement vector of block i; u₀ and v₀ are the translations of block centroid along X and Y axes, respectively; r₀ is the rigid rotation around point (x₀, y₀); ε_x, ε_y and γ_xy are the three strain components at point (x₀, y₀).

The displacement vector U=(u  v)^T of point (x, y) within block i is determined by a complete one order approximation function:

U=TD_i                                      (2)

where  T is the displacement transformation matrix and is defined as

(3)

2.2 Simultaneous equilibrium equations

According to the minimum energy principle, by taking derivatives with regard to the unknown displacement variables, the global equilibrium equations can be obtained as follows:

                                    (4)

Assuming that a block system consists of n blocks, we have

  (5)

where D_i and F_i (i=1, 2, …, n) are 6×1 sub-matrices, and D_i is the deformation variable of block i, while F_i is the load distributed to block i; K_ij (i, j=1, 2, …n) is a 6×6 sub-matrix, and K_ii is relevant to the material properties of block i, while K_ij (i≠j) is defined by the contact between blocks i and j.

2.3 Kinematics conditions

The solution of Eqn.(5) must satisfy the condition of non-penetration between blocks. If penetration happens, very stiff contact springs should be applied between the two blocks, while the number, the location and the direction of the contact springs depend on the contact type and contact status detected.

After the possible contact springs are appended between two contacting blocks, the sub-matrices of contact springs are computed and added to the matrices of the global equations at the corresponding positions, then Eqn.(5) is solved. According to the results, the requirements of non-penetration and non-tension are checked. If the requirements are satisfied, the computation of this time step is finished, or else, the contact springs should be applied where penetration occurs, or should be removed where tension occurs. This iteration goes on until the requirements mentioned above are satisfied, indicating that all the contact springs are appended correctly.

3 Proposed algorithm

3.1 Generation of joint network

According to the probability models for the geometrical parameters of the joints randomly existing in rock mass, the center locations, direction angles and trace lengths of these joints within the area of interest were obtained by using the Monte-Carlo technique^[20-21]. Commonly, the center locations of joints satisfy uniform or Poisson distribution, and the direction angles and trace lengths satisfy normal distribution. Herein, the center locations of joints are assumed to satisfy uniform distribution, and their direction angles and trace lengths satisfy normal distribution.

The steps of generating random joint network are presented as follows.

1) For each joint set, repeat the following steps.

2) Estimate the numberof the random joints in the area of interest by the formula:

                                     (6)

where S is the area of the computational model; d is the mean of the joint spacings; l is the mean of the joint trace lengths.

3) A sequence ofrandom numbers of uniform distribution is generated by

           (7)

where i=1, 2, 3, …, N; mod denotes complementation; x₀ is the initial value computed by function time( ), so it varies with time, and as a result, different number sequences can be generated at different time; r_i is the ith random number of uniform distribution.

Then, x coordinates of the center locations of N joints can be obtained by

x_i=(x_max-x_min)r_i+x_min                           (8)

where x_max is the maximum value of x coordinates of the joint center locations, and x_min is the minimum value.

In the same way, y coordinates of the joint center locations can also be obtained.

4) Based on 12N random numbers derived from Eqn.(7), the direction angles of N joints with normal distribution  can be determined by

         (9)

where σ_α is the mean square deviation; μ_α is the mean of the joint direction angles.

In the same way, the joint trace lengths can also be determined.

5) According to the obtained joint geometrical parameters, the coordinates of two end points of every joint can be computed.

6) For each joint, the two end points are checked if they are beyond the area of interest. If yes, the outside section is cut off, and correspondingly, the coordinates of the end point are changed.

3.2 Generation triangular block

In this work, the advancing front method, a widely- used FEM mesh generator^[22-24], was selected to generate triangular blocks. Beforehand, the following two preparatory tasks should be done: adding assistant lines (artificial joints) to connect the semi-penetrative joints with other joints or the model boundaries, and then searching all the simply connected domains in the area of interest.

The algorithm for generating triangular blocks is descripted as follows.

1) For each simply connected domain, repeat the following steps.

2) In terms of the assigned block size l, the domain boundaries are discretized into a series of nodes that form the initial front.

3) The shortest edge AB at the current front is selected as an active one, and at the left side of AB, a point C located at a distance of l away from points A and B, is determined.

4) An array composed of nodes N₁, N₂, …, N_p, that are at the current front in the range of 5l around point C, is constructed and ascendingly sorted in terms of their distances from point C, as shown in Fig.1(a).

5) If node N₁ cannot satisfy the following conditions: ＜1.5l, ＜1.5l, and triangle ABN₁ is counter- clockwise, point C will be chosen as N₁, the original N₁ will become N₂, and the rest will move in turn, as shown in Fig.1(b).

Fig.1 Sketch maps of generating triangular block: (a) Initial array of nodes; (b) Modified array of nodes

6) Select a node in the array from the first to the last until a valid triangle ABN_i is available. The conditions for judging the validation of the formed triangle are: any node in the array is not inside the triangle ABN_i, the midline of  does not intersect the current front, and N_i is at the left side of

7) Update the current front, and repeat steps 3)-6) until the front is null.

8) The generated mesh is smoothed by the Laplacian method, i.e. the coordinates of inside node take the centroid coordinates of the polygon which is composed of the triangular blocks adjacent to the node.

3.3 Rock failure analysis

In the process of generating triangular blocks, a lot of artificial joints come into being, and they provide the potential paths along which the cracks initiate and propagate. The artificial joints take the strength parameters of intact rock in order to represent the characteristics of the continuous region. Rock failure is induced by a series of rupture of some artificial joints, which is controlled by the maximum tensile stress criterion and Mohr-Coulomb criterion:

σ=-σ_t                                     (10)

τ=c+σtan φ                                  (11)

where σ and τ are the normal and tangential stresses on the artificial joint surface, respectively; σ_t is the tensile strength of rock mass; c and φ are the cohesive strength and friction angle of rock mass.

Once one of the above criterions is satisfied, the artificial joint will become the real joint.

4 Verification examples

The algorithm discussed above was implemented into the DDA code originally developed by SHI^[19]. In this section, the propagation and coalescence of random joints with different lengths and angles were simulated by using the modified DDA program.

4.1 Simulation of failure process of rock specimens with random joints

The rock specimens in the two examples have the same configuration, 10 cm in width and 17 cm in height. Specimen A contains a set of random joints distributing with the spacing of 5 cm, direction angles of N(45?, (5?)²) and trace lengths of N(2 cm, (0.2 cm)²). In specimen B, there exist two sets of random joints, of which one set has the same parameter distribution as that in specimen A, and the other set distributes with the spacing of 5 cm, direction angles of N(135?, (5?)²) and trace lengths of N(2 cm, (0.2 cm)²). The generated joint networks in the two specimens are shown in Fig.2, and the established computational models of triangular blocks are respectively shown in Fig.3(a) and Fig.4(a), in which the gray lines are artificial joints.

Fig.2 Random joints contained in specimens A (a) and B (b)

Fig.3 Computational model and failure process of rock specimen A: (a) Computational model; (b) Step 70; (c) Step 100; (d) Step 120; (e) Step 140; (f) Step 150

Fig.4 Computational model and failure process of rock specimen B: (a) Computational model; (b) Step 40; (c) Step 70; (d) Step 90; (e) Step 110; (f) Step 120

The mechanical properties of rock materials are as follows: density ρ=2 500 kg/m³, mean elastic modulus E=70 GPa, mean Poisson ratio γ=0.28, tensile strength σ_t=18 MPa, cohesive strength c=26 MPa. The strength parameters of joints are: friction angle φ=30?, no cohesive strength and no tensile strength.

At the bottom of the specimens, the displacement in direction Y is confined; while at the top, a line loading is applied and its value increases gradually. The left and right boundaries of the specimens are free. The numerical results of failure process of specimens A and B under uniaxial compression are shown in Figs.3 and 4, respectively.

From Figs.3 and 4, the total features of crack growth and rock failure under compression can be seen clearly. At the initial stage of loading, some wing cracks grow from the crack tips. With the gradual increment of loading, these wing cracks continuously propagate in the direction approximatively parallel to the axial loading, and some secondary cracks appear at the crack tips. Subsequently, the wing and secondary cracks interact and coalesce with each other. And finally, the coalesced cracks grow to the end surfaces of specimens. The similar phenomena are observed in many laboratorial and numerical experiments of rock specimen containing few cracks under compression^[25]. Besides, it is worth noting that although some joints are near the free boundaries, they do not grow to the boundaries, but to the direction of axial loading.

4.2 Failure process of surrounding rock mass of underground chamber

As shown in Fig.5, a chamber of 4.0 m×2.5 m is excavated in semi-continuous rock mass in which there exist two sets of random joints: one set distributes with the spacing of 3 m, direction angles of N(45?, (10?)²) and trace lengths of N(2.5 m, (0.5 m)²), and the other set distributes with the spacing of 3 m, direction angles of N(135?, (10?)²) and trace lengths of N(2.5 m, (0.5 m)²). The computational model with dimensions of 16 m×10 m is constructed, as shown in Fig.6.

Fig.5 Configuration of underground chamber

Fig.6 Computational model of chamber

The mechanical properties of rock materials are follows: density ρ=2 600 kg/m³, elastic modulus E=20 GPa, Poisson ratio γ= 0.2, tensile strength σ_t=12 MPa, cohesive strength c=22 MPa, and friction angle φ=35?. The joints are strengthless.

An increasing line loading is applied at the top boundary of the model, and the other three boundaries are simply supported, i.e. the displacement in one direction is confined. Fig.7 shows the simulated failure process of surrounding rock mass of the chamber.

From Fig.7, it can be seen that cracks initiate from the tips of joints, and grow to the nearby joints. Ultimately, with the growth of cracks, some joints around the chamber are coalesced, and several individual rock blocks come into being. Under the action of gravity, these rock blocks will crumble from the native rock and move into the chamber. Moreover, it can be observed that some cracks arise near the left and right boundaries, but they do not propagate in further steps and the surrounding rock mass in the far field is relatively intact, while coalescence occurs between the joints around the chamber, and as a result, the surrounding rock mass in the near field is fragmentized. This observation indicates that the free face may cause the joints to grow more easily. Therefore, it can be concluded that for the case of multi cracks, the constraint condition around crack is an important influencing factor to determine the direction of crack growth.

Fig.7 Failure process of surrounding rock mass of chamber: (a) Step 60; (b) Step 90; (c) Step 100; (d) Step 110

5 Conclusions

1) By modifying the original DDA method, a simple methodology for modeling the complicated failure process of semi-continuous rock mass is presented.

2) Crack initiates from the crack tip, and its growth direction depends upon the loading and constraint conditions.

3) The whole process of rock fragmentation can be easily simulated by the proposed method, and the simulated results can reproduce the complicated phenomenon of rock failure, such as new crack generation, crack branching, multi-crack interaction, which cannot be simulated by the continuum-based methods.

4) The crack path obtained from the proposed algorithm severely depends on the triangular block boundaries, for it can only grow along the artificial joints. This limitation can be eliminated in the further research.

[1] WONG R H C, CHAU K T, TANG C A, LIN P. Analysis of crack coalescence in rock-like materials containing three flaws (Part I): Experimental approach [J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(7): 909-924.

[2] NARA Y, KANEKO K. Study of subcritical crack growth in andesite using the Double Torsion test [J]. International Journal of Rock Mechanics and Mining Sciences,2005, 42(4): 521-530.

[3] XIA K W, CHALIVENDRA V B, ROSAKIS A J. Observing ideal “self-similar” crack growth in experiments [J]. Engineering Fracture Mechanics,2006, 73(18): 2748-2755.

[4] KAMAYA M. Growth evaluation of multiple interacting surface cracks (Part I): Experiments and simulation of coalesced crack [J]. Engineering Fracture Mechanics,2008, 75(6): 1336-1349.

[5] HOLZHAUSEN G R, JOHNSON A M. Analysis of longitudinal splitting of uniaxially compressed rock cylinder [J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1979, 16(3): 163-177.

[6] NEMAT-NASSER S, HORII H. Compression-induced nonlinear crack extension with application to splitting, exfoliation, rockburst [J]. Journal of Geophysical Research, 1982, 87(B8): 6805-6821.

[7] WANG E Z, SHRIVE N G. Brittle fracture in compression: Mechanisms, models and criteria [J]. Engineering Fracture Mechanics, 1995, 52(6): 1107-1126.

[8] BOBET A, EINSTEIN H H. Fracture coalescence in rock-type material under uniaxial and biaxial compression [J]. International Journal of Rock Mechanics and Mining Sciences, 1998, 35(7): 836-888.

[9] WONG R H C, CHAU K T. Crack coalescence in a rock-like material containing two cracks [J]. International Journal of Rock Mechanics and Mining Sciences, 1998, 35(2): 147-164.

[10] JIAO Y Y, FAN S C, ZHAO J. Numerical investigation of joint effect on shock wave propagation in jointed rock masses [J]. Journal of Testing and Evaluation, 2005, 33(3): 197-203.

[11] JIAO Y Y, ZHANG X L, ZHAO J, LIU Q S. Viscous boundary of DDA for modeling stress wave propagation in jointed rock [J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(7): 1070-1076.

[12] GOLSHANI A, OKUI Y, ODA M, TAKEMURA T. A micromechanical model for brittle failure of rock and its relation to crack growth observed in triaxial compression tests of granite [J]. Mechanics of Materials,2006, 38(4): 287-303.

[13] CHARALAMBIDES P G, MCMEEKING R M. Finite element method simulation of crack propagation in a brittle microcracking solid [J]. Mechanics of Materials, 1987,6(1): 71-87.

[14] TAN X C, KOU S Q, LINDQVIST P A. Application of the DDM and fracture mechanics model on the simulation of rock breakage by mechanical tools [J]. Engineering Geology, 1998, 49(3): 277-284.

[15] XU Y, SAIGAL S. An element free Galerkin formulation for stable crack growth in an elastic solid [J]. Computer Methods in Applied Mechanics and Engineering,1998, 154(3/4): 331-343.

[16] CHEN Wei, LI Shi-hai. Deformable and rupturable block model for 3D distinct element method [J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(4): 545-549. (in Chinese)

[17] JIAO Y Y, LIU Q S, LI S C. A three-dimensional numerical model for simulating deformation and failure of blocky rock structures [J]. Key Engineering Materials, 2006, 306/308: 1391-1396.

[18] JIAO Y Y, ZHANG X L, WANG S L, WU H Z. On using discrete particle approaches for simulating perforation process of concrete slab by hard projectile [J]. Key Engineering Materials, 2007, 353/358: 2973-2976.

[19] SHI G H. Discontinuous deformation analysis: A new numerical model for the statics and dynamics of block system [D]. Berkeley: University of California, 1988.

[20] S?NCHEZ S, CRIADO R, VEGA C. A generator of pseudo-random numbers sequences with a very long period [J]. Mathematical and Computer Modelling,2005, 42(7-8): 809-816.

[21] DENG L Y, GUO R, LIN D K J, BAI F S. Improving random number generators in the Monte-Carlo simulations via twisting and combining [J]. Computer Physics Communications,2008, 178(6): 401-408.

[22] UEMURA K, SAITO T. Automatic mesh generation for FEM simulation of wind flow around tall buildings [J]. Journal of Wind Engineering and Industrial Aerodynamics,1993, 46/47: 357-362.

[23] DENG Jian-hui. Adaptive finite element analysis of jointed rocks—method and implementation [D]. Wuhan: Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, 1994. (in Chinese)

[24] GU?NETTE R, FORTIN A, KANE A, H?TU J F. An adaptive remeshing strategy for viscoelastic fluid flow simulations [J]. Journal of Non-Newtonian Fluid Mechanics,2008, 153(1): 34-45.

[25] JIAO Yu-yong, ZHANG Xiu-li, LIU Quan-sheng, CHEN Wei-zhong. Simulation of rock crack propagation using discontinuous deformation analysis method [J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 26(4): 682-691. (in Chinese)

Foundation item: Projects(50479071, 40672191) supported by the National Natural Science Foundation of China; Project(SKLZ0801) supported by the Independent Research Key Project of State Key Laboratory of Geomechanics and Geotechnical Engineering; Project(SKLQ001) supported by the Independent Research Frontier Exploring Project of State Key Laboratory of Geomechanics and Geotechnical Engineering

Received date: 2008-04-12; Accepted date: 2008-06-25

Corresponding author: JIAO Yu-yong, Professor, PhD; Tel: +86-27-87198299; E-mail: yyjiao@whrsm.ac.cn

(Edited by CHEN Wei-ping)