Simulation of failure process of jointed rock

J. Cent. South Univ. Technol. (2008) 15: 888- -894包SpringerDOI: 10.1007/11771-008-0162-0Simulation of failure process of jointed rockZHANG 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 modifed discontinuous deformation analysis (DDA) algorithm was proposed to simulate the failure behavior of jointedrock. 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. Inthe process of generating blocks, numerous artificial joints came into being, and once the stress states at some artificial joints satisfythe failure criterion given beforehand, artificial joints will turn into real joints. In this way, the whole fragmentation process of rockmass can be replicated. The algorithm logic was described in detail, and several numerical examples were carried out to obtain someinsight into the failure behavior of rock mass containing random joints. From the numerical results, it can be found that the crackinitiates from the crack tip, the growth direction of the crack depends upon the loading and constraint conditions, and the proposedmethod 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; advancingfront method; triangular block systemrelaxation upon satisfaction of a critical energy release1 Introductionrate criterion. TAN et al14] used a displacementdiscontinuity method (DDM) program coupled with aNatural rock mass exists in compression stressmodified energy criterion to simulate the development ofenvironment, and its failure is related to the initiation,cracks, and the numerical results showed that the crackpropagation and coalescence of interior cracks. In themay cither propagate to the free surface or rest in thepast decades, significant researches on the topic of crackrock subsurface, which is dependent on the rock andgrowthhavebeendone experimentallyh- ,fracture properties, loading force. XU and SAIGAL!S1analytically'- ”and numericallyt-i8. However, due toestablished the discrete formulation for stable crackthe mechanical and geometrical complexity of thisgrowth in an elastic solid using the element free Galerkinproblem, it has not been yet solved fully and perfectly.(EFG) method based on the moving least-squaresAmong all the available means, numerical approachapproximations. CHEN and LIt6 presented a 3D distinctis proven to be an effective tool for simulating rockelement method for block deformation and rupturefailure because the loading and boundary conditions ofanalysis. However, few efforts on the propagationion andthe model can be controlled strictly, and some parameters, coalescence of multi-cracks have succeeded, and nosuch as crack length and crack location can be easilyresults on the failure process of rock with stochasticchanged so as to study their influences on the mode 0cracks have been reported. The main reason is that thecrack growth. So far, by using some numerical methodscommonly-used numerical methods, such as FEM, BEMembedded with the fracture theory, a number of attemptsor meshfree method, always encounter mathematicalhave been made to study the mechanism of crackdifficulties in dealing with multi-crack problem, whilegrowth!I2- 18]. CHARALAMBIDES and MCMEEKING(I3]discrete element method requires further efforts on crackadopted a finite element method to simulate crackgrowth simulation. In this work, on the basis of thepropagation for studying the effects of microcracking atdiscontinuous deformation analysis (DDA) methodfacets as a damage process in brittle composites and theproposed by SHI，an approach was presented toconsequence for material toughness. In their study, cracksimulate the propagation and coalescenceofpropagation was achieved by successive tip nodemulti-cracks.Foundation item: Projects(50479071, 40672191) supported by the National Natural Science Found中国煤化工uportedbotheIndependent Research Key Project of State Key Laboratory of Geomechanicroject(SKLQ001)supported by the Independent Research Frontier Exploring Project of StateMYHC N M H Gand GeotechnicalEngineeringReceived date: 2008- 04- 12; Accepted date: 2008 - 06 -25Corresponding author: JIAO Yu-yong, Professor, PhD; Tel: +86- 27- -87 198299; E-mail: yyjiao@whrsm.ac.cn.J. Cent. South Univ. Technol. (2008) 15: 888- 894889In the proposed method, some modifications werecan be obtained as follows:made for enabling DDA to handle semi-penetrative jointsKD=F(4)andperform rock failure analysis. Firstly, thMonte-Carlo technique was employed to generateAssuming that a block system consists of n blocks,random joints whose geometrical parameters satisfiedwe havesome stochastic distribution laws obtained from geology[K1 K12.... KInDj 1F]survey. Then, on the generated joint network, a popularD2|F2|FEM mesh generation method, namely the advancingK=(5)front method, was applied to automatically generatingtriangular blocks. In the block generation process,_Kn Kn2 ... Kmm_D,」Fn」numerous artificial joints came into being. Differentfrom the real joints, the strength of the artificial jointswhere D; and F:(i=l, 2,”，n) are 6X 1 sub-matrices,was much higher, and the strength of that intact rock wasand Di is the deformation variable of block i, while F isemployed. Through the rupture analysis of the artificialthe load distributed to block i; Kj (i,j=l, 2, ..n) is ajoints, the rock failure process can be simulated in a6X6 sub-matrix, and Ki; is relevant to the materialproperties of block i, while K; (i≠j) is defined by the .straightforward way.contact between blocks i andj.2 DDA method2.3 Kinematics conditionsThe solution of Eqn.(5) must satisfy the condition ofDDA is a numerical approach based on block theoryand minimum energy principle. The movement of eachnon-penetration between blocks. If penetration happens,block is described by Newton's second law of motion.very stiff contact springs should be applied between theAnd, the penetration of one block into the other is nottwo blocks, while the number, the location and theallowed. The simultaneous equilibrium equations aredirection of the contact springs depend on the contactestablished by minimizing the total potential energy oftype and contact status detected.the block system due to different mechanisms.After the possible contact springs are appendedbetween two contacting blocks， the sub-matrices of2.1 Function of displacementcontact springs are computed and added to the matricesIn DDA, each block of arbitrary geometry has sixof the global equations at the corresponding positions,degrees of freedom, among which three components arehen Eqn.(5) is solved. According to the results, therigid body motion terms and the other three are constantrequirements of non-penetration and non-tension arestrain terms. Therefore, the deformation variable of blockchecked. If the requirements are satisfied, thei can be written ascomputation of this time step is finished, or else, thecontact springs should be applied where penetrationD=(uovoroεxεy2(1occurs, or should be removed where tension occurs. Thiswhere D; is the displacement vector of block i; uo anditeration goes on until the requirements mentioned aboveVo are the translations of block centroid along X and Yare satisfied, indicating that all the contact springs areaxes, respectively; ro is the rigid rotation around point (x0,appended correctly.yo); Eex, Ey and Yxy are the three strain components at point(x0, yo).3 Proposed algorithmThe displacement vector U=(u ) of point (x, y)within block i is determined by a complete one order3.1 Generation of joint networkapproximation function:According to the probability models for theU=TD, .(2geometrical parameters of the joints randomly existing inwhere T is the displacement transformation matrix androck mass, the center locations, direction angles andis defined asrace lengths of these joints within the area of interestwere obtainedII by using the Monte-Carlo technique'r=|[1 0 -(v-yo) x-x0 0 (y-yo)/27Commonly, the center locations of joints satisfy uniform0 1 x-xo0 y-yo (x-xo)/2」or Poisson distribution, and the direction angles and trace(3)lengths satisfy normal distribution. Herein, the centerlocations of joints are assumed to satisfy uniform2.2 Simultaneous equilibrium equationsdistribution, anc中国煤化工trace lengthsAccording to the minimum energy principle， bysatisfy normal CCNMHGtaking derivatives with regard to the unknownThe steps oI generating ranaom JoInI network aredisplacement variables, the global equilibrium equationspresented as follows..890J. Cent. South Univ. Technol. (2008) 15: 888- 8941) For each joint set, repeat the following steps.The algorithm for generating triangular blocks is2) Estimate the number N of the random joints indescripted as follows.the area of interest by the formula:1) For each simply connected domain, repeat thefollowing steps.N=S(62) In terms of the assigned block size l, the domainboundaries are discretized into a series of nodes thatwhere S is the area of the computational model; d is theform the initial front.mean of the joint spacings; l is the mean of the joint trace3) The shortest edge AB at the current front islengths.selected as an active one, and at the left side of AB, a3) A sequence of N random numbers of uniformpoint C located at a distance of l away from points A anddistribution is generated byB, is determined.[x; =(2 053x;-1 + 13 849) (mod 65 536)4) An array composed of nodes N, N2, ", Np, that{r; =x;/65 536are at the current front in the range of 5l around point C,is constructed and ascendingly sorted in terms of their|xo=time()distances from point C, as shown in Fig.1(a).where i=1, 2, 3, ., N; mod denotes complementation;5) If node N cannot satisfy the following conditions:x0 is the initial value computed by function time( ), so itAN1<1.5l, BN1<1.5Il, and triangle ABN is counter-varies with time, and as a result, different numberclockwise, point C will be chosen as N, the original N;sequences can be generated at different time; r; is the ithwill become N2, and the rest will move in turmn, as shownrandom number of uniform distribution.in Fig.1(b).Then, x coordinates of the center locations of Njoints can be obtained by__Nsx=(Xxmax~ -Xmi)r/+xmin(8Vswhere Xmax is the maximum value of x coordinates of theN2joint center locations, and Xmin is the minimum value.51In the same way, y coordinates of the joint center~C fBlocations can also be obtained.4) Based on 12N random numbers derived fromNNEqn.(7), the direction angles of N joints with normal(adistribution N(μa,σa) can be determined by'NNoQ; =σa2(2i-)+2j- 2-1-+2)-)+ Ha(9)N4j=where σa is the mean square deviation; μa is the mean ofV3the joint direction angles.CIN tIn the same way, the joint trace lengths can also bedetermined.NsAccording to the obtained jointgeometrical(b)parameters, the coordinates of two end points of everyFig.1 Sketch maps of generating triangular block: (a) Initialjoint can be computed.array of nodes; (b) Modified array of nodes6) For each joint, the two end points are checked ifthey are beyond the area of interest. If yes, the outside6) Select a node in the array from the first to the lastsection is cut off, and correspondingly, the coordinates ofuntil a valid triangle ABN; is available. The conditionsthe end point are changed.for judging the validation of the formed triangle are: anynode in the array is not inside the triangle ABN, the3.2 Generation triangular blockIn this work, the advancing front method, a widely-midline of AB does not intersect the current front, andused FEM mesh generatorl22- , was selected to generateN; is at the left side of AB.7) Update the current front, and repeat steps 3)- -6)triangular blocks. Beforehand, the following twpreparatory tasks should be done: adding assistant linesuntil the front is null.(artificial joints) to connect the semi- penetrative joints8) The中国煤化工hed by thewith other joints or the model boundaries, and thenLaplacian me:rYHC NMH G。二inside nodesearching all the simply connected domains in the area oftake the centrol-.gon which isinterest.composed of the triangular blocks adjacent to the node..J. Cent. South Univ. Technol. (2008) 15: 888- -894891 _3.3 Rock failure analysisThe rock specimens in the two examples have theIn the process of generating triangular blocks, a lotsame configuration, 10 cm in width and 17 cm in height.of artificial joints come into being, and they provide theSpecimen A contains a set of random joints distributingpotential paths along which the cracks initiate anwith the spacing of 5 cm, direction angles of N(45*, (5°)2)propagate. The artificial joints take the strengthand trace lengths of N(2 cm, (0.2 cm). In specimen B,parameters of intact rock in order to represent ththere exist two sets of random joints, of which one setcharacteristics of the continuous region. Rock failure ishas the same parameter distribution as that in specimen A,induced by a series of rupture of some artificial joints,and the other set distributes with the spacing of 5 cm,which is controlled by the maximum tensile stressdirection angles of N(135", (5")}) and trace lengths ofcriterion and Mohr-Coulomb criterion:N(2 cm, (0.2 cm). The generated joint networks in theσ=-σt(10)two specimens are shown in Fig.2, and the establishedtFc+otan φ(11)computational models of triangular blocks arewhere σ and T are the normal and tangential stresses onrespectively shown in Fig.3(a) and Fig.4(a), in which thethe artificial joint surface, respectively; σ is the tensilegray lines are artificial joints.strength of rock mass; c and φ are the cohesive strengthand friction angle of rock mass.(aOnce one of the above criterions is satisfied, theartificial joint will become the real joint.///4 Verification examplesThe algorithm discussed above was implementedinto the DDA code originally developed by SHr19. Inthis section, the propagation and coalescence of randomjoints with different lengths and angles were simulatedby using the modified DDA program.4.1 Simulation of failure process of rock specimenswith random jointsFig.2 Random joints contained in specimens A (a) and B (b)a)_2 (b(c//d)F (e(f中国煤化工CHCNM HGFig.3 Computational model and failure process of rock specimen A: (a) Computational Tnroucl, (0) siep 1v, (心) Dicp 1U0; (d) Step 120;(e) Step 140; () Step 150.892J. Cent. South Univ. Technol. (2008) 15: 888- -894a).(b)>)(d)(e)(0.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 120The mechanical properties of rock materials are asSubsequently, the wing and secondary cracks interactfollows: density p=2 500 kg/m', mean elastic modulus .and coalesce with each other. And finally, the coalescedE=70 GPa, mean Poisson ratio y=0.28, tensile strengthcracks grow to the end surfaces of specimens. The0=18 MPa, cohesive strength c=26 MPa. The strengthsimilar phenomena are observed in many laboratorial andparameters of joints are: friction angle φ =30"，nonumerical experiments of rock specimen containing fewcohesive strength and no tensile strength.cracks under compression's. Besides, it is worth notingAt the bottom of the specimens, the displacement inthat although some joints are near the free boundaries,direction Y is confined; while at the top, a line loading isthey do not grow to the boundaries, but to the directionapplied and its value increases gradually. The left andof axial loading.right boundaries of the specimens are free. Thnumerical results of failure process of specimens A and B4.2 Failure process of surrounding rock mass ounder uniaxial compression are shown in Figs.3 and 4,underground chamberrespectively.As shown in Fig.5, a chamber of 4.0 mX2.5 m isFrom Figs.3 and 4, the total features of crackexcavated in semi-continuous rock mass in which theregrowth and rock failure under compression can be seenexist two sets of random joints: one set distributes withclearly. At the initial stage of loading, some wing cracksthe spacing of 3 m, direction angles of N(45", (10°)3) and .grow from the crack tips. With the gradual increment oflengths中国煤化工the other setloading, these wing cracks continuously propagate in thedistributes withion angles ofdirection approximatively parallel to the axial loading,N(135°, (10*)2)YHCNMHGm,(0.5m).and some secondary cracks appear at the crack tips.The computational model with dimensions of 16 mX 10.J. Cent. South Univ. Technol. (2008) 15: 888- -894893 .m is constructed, as shown in Fig.6. .the tips of joints, and grow to the nearby joints.The mechanical properties of rock materials areUltimately，with the growth of cracks, some jointsfollows: density ρ=2 600 kg/m', elastic modulus E=20around the chamber are coalesced, and several individualGPa, Poisson ratio F 0.2, tensile strength σ=12 MPa,rock blocks come into being. Under the action of gravity,cohesive strength c=22 MPa, and friction angle φ=35°.these rock blocks will crumble from the native rock andThe joints are strengthless.move into the chamber. Moreover, it can be observedAn increasing line loading is applied at the topthat some cracks arise near the left and right boundaries,boundary of the model, and the other three boundariesut they do not propagate in further steps and thare simply supported, i.e. the displacement in onesurrounding rock mass in the far field is relatively intact,direction is confined. Fig.7 shows the simulated failurewhile coalescence occurs between the joints around theprocess of surrounding rock mass of the chamber.chamber, and as a result, the surrounding rock mass inFrom Fig.7, it can be seen that cracks initiate fromthe near field is fragmentized. This observation indicatesthat the free face may cause the joints to grow moreeasily. Therefore, it can be concluded that for the case ofmulti cracks, the constraint condition around crack is animportant influencing factor to determine the direction ofcrack growth.5 Conclusionsmethodology for modeling the complicated failureFig.5 Configuration of underground chamberprocess of semi-continuous rock mass is presented.2) Crack initiates from the crack tip, and its growthdirection depends upon the loading and constraintconditions.3) The whole process of rock fragmentation can beeasily simulated by the proposed method, and thesimulated results can reproduce the complicatedphenomenon of rock failure， such as new crackgeneration，crack branching， multi-crack interaction,which cannot be simulated by the continuum-basedmethods.Fig.6 Computational model of chamber4) The crack path obtained from the proposed(a)(b(e)中国煤化工MYHCNMH GFig.7 Failure process of surrounding rock mass of chamber: (a) Step 60; (b) Step 90; (c) Step 100; (d) Step 110.894J. Cent. South Univ. Technol. (2008) 15: 888- 894algorithm severely depends on the triangular blockcrack growth observed in triaxial compression tests of granite []Mechanics of Materials, 2006, 38(4): 287 -303.boundaries, for it can only grow along the artificial joints.13] CHARALAMBIDES P G, MCMEEKING R M. Finite elementThis limitation can be eliminated in the further research.method simulation of crack propagation in a brittle microcrackingsolid [I]. Mechanics of Materials, 1987, 6(1): 71-87.References[14] TAN x C, KOU s Q, LINDQVIST P A. Application of the DDM andfracture mechanics model on the simulation of rock breakage by[1] WONGR HC, CHAU K T, TANGC A, LIN P. Analysis of crackmechanical tools [J]. Engineering Geology, 1998, 49(3): 277-284.coalescence in rock-like materials containing three flaws (Part D):15] XU Y, SAIGAL S. An element free Galerkin formulation for stableExperimental approach []. International Journal of Rock Mechanicscrack growth in an elastic solid [D]. Computer Methods in Appliedand Mining Sciences, 2001, 38(7): 909 -924.Mechanics and Engineering, 1998, 154(3/4): 331-343.[2NARA Y, KANEKO K. Study of subcritical crack growth in andesite16] CHEN Wei, LI Shi-hai. Deformable and rupturable block model forusing the Double Torsion test [J]. Intemnational Journal of Rock3D distinct clement method [U]. Chinese Journal of Rock MechanicsMechanics and Mining Sciences, 2005, 42(4): 521-530.and Engineering, 2004, 23(4): 545 -549. (in Chinese)[3] XIA K W, CHALIVENDRA V B, ROSAKIS A J. Observing ideal17] JIAO Y Y, LIU QS, LI S C. A three-dimensional numerical model“self- similar”crack growth in experiments [0]. Engineering Fracturefor simulating deformation and failure of blocky rock stuctures [].Mechanics, 2006, 73(18): 2748- -2755.Key Engineering Materials, 2006, 306/308: 1391-1396.[4] KAMAYA M. Growth evaluation of multiple interacting surface18] JIAO Y y, ZHANG x L, WANG s L, WU H z. On using discretecracks (Part D): Experiments and simulation of coalesced crack [].particle approaches for simulating perforation process of concreteEngineering Fracture Mechanics, 2008, 75(6): 1336- :1349.slab by hard proijctile []. Key Engineering Materials, 2007, 353/358:[5] HOLZHAUSEN G R, JOHNSON A M. Analysis of longitudinal2973- 2976.spilting of uniaxially compressed rock cylinder [] InternationalI9] SHI G H. Discontinuous deformation analysis: A new numericalJournal of Rock Mechanics and Mining Sciences and Geomechanicsmodel for the statics and dynamics of block system [D] Berkeley:Abstracts, 1979, 16(3): 163-177.University of California, 1988.[6] NEMAT-NASSER s, HORII H. Compression induced nonlinear20] SANCHEZ S, CRIADO R, VEGA C. A generator of pseudo-randomcrack extension with application to spltting, exfoliation, rockburstnumbers sequences with a very long period []. Mathematical and[]. Journal of Geophysical Research, 1982, 87(B8): 6805- -6821.Computer Mdelling, 2005, 42(7- 8): 809- -816.[7] WANG E z, SHRIVE N G Brittle fracture in compression:[21] DENG L y, GUO R, LIN D K J, BAI F s. Improving random numberMechanisms, models and criteria []. Engineering Fracturegenerators in the Monte-Carlo simulations via twisting andMechanics, 1995, 52(6): 1107-1126.combining [] Computer Physics Communications, 2008, 178(6):[8] BOBET A, EINSTEIN H H. Fracture calasecence in rock-type401-408.matrial under uniaxial and biaxial compression [J]. Intermational2] UEMURA K, SAITO T. Automatic mesh generation for FEMJournal of Rock Mechanics and Mining Sciences, 1998, 35(7):836- 888. .simulation of wind flow around tall buildings [U]. Journal of WindEngineering and Industrial Aerodynamies, 1993, 46147: 357-362.[9] WONG R H C, CHAU K T. Crack coalescence in a rocklike23] DENG Jian-hui. Adaptive finite element analysis of jointed rocks-material containing two cracks []. International Journal of RockMechanics and Mining Sciences, 1998, 35(2): 147-164.method and implementation [D]. Wuhan: Institute of Rock and Soil[10] JIAO Y y, FAN s C, ZHAO J. Numerical ivstiation of joint efeMechanics, Chinese Academy of Sciences, 1994. (in Chinese)on shock wave propagation in jointed rock masses []. Journal of24] GUENETTE R, FORTIN A, KANE A, HETU J F. An adaptiveTesting and Evaluation, 2005, 33(3): 197- 203.remeshing strategy for viscoelastic fluid flow simulations [0]. Journal[11] JIAO Y Y, ZHANG X L, ZHAO J, LIU Q s. Viscous boundary ofof Non-Newtonian Fluid Mechanics, 2008, 153(1): 34-45.DDA for modeling stress wave propagation in jointed rock [].[25] JIAO Yu-yong, ZHANG Xiu-i, LIU Quan sheng. CHEN Wei- zhong.Intermational Journal of Rock Mechanics and Mining Sciences, 2007,Simulation of rock crack propagation using discontinuous44(7): 1070-1076.deformation analysis method []. Chinese Journal of Rock Mechanics[12] GOLSHANI A， OKUI Y, ODA M, TAKEMURA T. Aand Engineering, 2006, 26(4): 682- 691. (in Chinese)micromechanical model for brittle failure of rock and its relation to(Edited by CHEN Wei-ping)中国煤化工MHCNMH G.