A bonded-particle model for rock

ARTICLE IN PRESS International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 www.elsevier.com/locate/ijrmms A bonded-particle model for rock D.O. Potyondya,,1, P.A. Cundallb a University of Toronto, Department of Civil Engineering and Lassonde Institute, Toronto, Ont., Canada M5S 1A4 b Itasca Consulting Group Inc., 111 Third Avenue South, Suite 450, Minneapolis, MN 55401, USA Accepted 9 September 2004 Available online 27 October 2004 Abstract A numerical model for rock is proposed in which the rock is represented by a dense packing of non-uniform-sized circular or spherical particles that are bonded together at their contact points and whose mechanical behavior is simulated by the distinctelement method using the two- and three-dimensional discontinuum programs PFC2D and PFC3D. The microproperties consist of stiffness and strength parameters for the particles and the bonds. Damage is represented explicitly as broken bonds, which form and coalesce into macroscopic fractures when load is applied. The model reproduces many features of rock behavior, including elasticity, fracturing, acoustic emission, damage accumulation producing material anisotropy, hysteresis, dilation, post-peak softening and strength increase with confinement. These behaviors are emergent properties of the model that arise from a relatively simple set of microproperties. A material-genesis procedure and microproperties to represent Lac du Bonnet granite are presented. The behavior of this model is described for two- and three-dimensional biaxial, triaxial and Brazilian tests and for two-dimensional tunnel simulations in which breakout notches form in the region of maximum compressive stress. The sensitivity of the results to microproperties, including particle size, is investigated. Particle size is not a free parameter that only controls resolution; instead, it affects the fracture toughness and thereby influences damage processes (such as notch formation) in which damage localizes at macrofracture tips experiencing extensile loading. r 2004 Elsevier Ltd. All rights reserved. Keywords: Rock fracture; Distinct-element method; Numerical model; Micromechanics 1. Introduction In this paper, we argue that rock behaves like a cemented granular material of complex-shaped grains in which both the grains and the cement are deformable and may break, and that such a conceptual model can, in principle, explain all aspects of the mechanical behavior. Various numerical models have been proposed that mimic such a system with differing levels of fidelity. The bonded-particle model for rock (referred to hereafter as the BPM) directly mimics this system and Corresponding author. Tel.: +1 416 978 3115; fax: +1 416 978 6813. E-mail addresses: [email protected] (D.O. Potyondy), [email protected] (P.A. Cundall). 1 Formerly with Itasca. 1365-1609/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2004.09.011 thus exhibits a rich set of emergent behaviors that correspond very well with those of real rock. The BPM provides both a scientific tool to investigate the micromechanisms that combine to produce complex macroscopic behaviors and an engineering tool to predict these macroscopic behaviors. The mechanical behavior of rock is governed by the formation, growth and eventual interaction of microcracks. Microscopic observations of rock [1–5] reveal detailed information about initial defects and loadinduced cracks, such as length, density, aspect ratio and orientation. Acoustic-emission (AE) based observations of rock [6] record the acoustic signals that are generated spontaneously from this microcracking, thereby providing information about the size, location and deformation mechanisms of the acoustic events as well as properties of the medium through which the acoustic ARTICLE IN PRESS 1330 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 waves travel (e.g., velocity, attenuation and scattering). Experimental observations reveal that most of the compression-induced cracks nucleate at the initial defects, such as grain boundaries or crack-like, lowaspect-ratio cavities, and that all compression-induced cracks are almost parallel to the direction of the maximum compression. The micromechanism responsible for the formation of these cracks is not understood fully; however, many possible models exist (e.g., sliding deformation between the faces of initial defects inclined to the direction of maximum compression leading to formation of a wing crack at each of the two defect tips) that can reproduce many of the essential features of the brittle fracture phenomenon [1,7–15]. Kemeny [10] asserts that: Although the actual growth of cracks under compression may be due to many complicated mechanisms, as revealed by laboratory tests, it appears that they can all be approximated by the crack with a central point load. The origins of these ‘‘point’’ loads in rock under compression are small regions of tension that develop in the direction of the least principal stress. Kemeny and Cook [11] emphasize the importance of producing compression-induced tensile cracks: Laboratory testing of rocks subjected to differential compression have revealed many different mechanisms for extensile crack growth, including pore crushing, sliding along pre-existing cracks, elastic mismatch between grains, dislocation movement, and hertzian contact. Because of the similarity in rock behavior under compression in a wide range of rock types, it is not surprising that micromechanical models have many similarities. This may explain the success of models based on certain micromechanisms (such as the sliding crack and pore models) in spite of the lack of evidence for these mechanisms in microscopic studies. One mechanism [16] for the formation of compression-induced tensile cracks is shown in Fig. 1(c), in which a group of four circular particles is forced apart by axial load, causing the restraining bond to experience tension. These axially aligned ‘‘microcracks’’ occur during the early loading stages of compression tests on bonded assemblies of circular or spherical particles. Similar crack-inducing mechanisms occur even when different conceptual models for rock microstructure are used. For example, ‘‘wedges’’ and ‘‘staircases’’ also induce local tension if angular grains replace circular grains—see Fig. 1(a) and (b). In addition to the formation and growth of microcracks, the eventual interaction of these cracks is necessary to produce localization phenomena such as axial splitting or rupture-zone formation during unconfined or confined compression tests. Thus, any model intending to Fig. 1. Physical mechanisms for compression-induced tensile cracking (a and b) and idealization as bonded assembly of circular particles (c). reproduce these phenomena must allow the microcracks to interact with one another. Rock can be represented as a heterogeneous material comprised of cemented grains. In sedimentary rock, such as sandstone, a true cement is present, whereas in crystalline rock, such as granite, the granular interlock can be approximated as a notional cement. There is much disorder in this system, including locked-in stresses produced during material genesis, deformability and strength of the grains and the cement, grain size, grain shape, grain packing and degree of cementation (i.e., how much of the intergrain space is filled with cement). All of these items influence the mechanical behavior, and many of them evolve under load application. The mechanical behavior of both rock and a BPM is driven by the evolution of the force-chain fabric, as will be explained here with the aid of Fig. 2. An applied macroscopic load is carried by the grain and cement skeleton in the form of force chains that propagate from one grain to the next across grain contacts, some of which may be filled with cement. The force chains are similar to those that form in a granular material [17]. The cement-filled contacts experience compressive, tensile and shear loading and also may transmit a bending moment between the grains, whereas the empty contacts experience only compressive and shear loading. Thus, applied loading produces heterogeneous force transmission and induces many sites of tension/ compression oriented perpendicular to the compression/ tension direction—for the reasons illustrated in Fig. 1(c). In addition, the force chains are highly nonuniform, with a few high-load chains and many lowload chains. The chain loads may be much higher than the applied loads, such that a few grains will be highly loaded while others will be less loaded or carrying no load (see Fig. 5(b))—because the forces will arch around these grains, thereby forming chains with a fabric that is larger than the grain size, as seen in the compression ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Fig 2. Force and moment distributions and broken bonds (in magenta) in a cemented granular material with six initial holes idealized as a bonded-disk assembly in the post-peak portion of an unconfined compression test. (Blue is grain–grain compression, while black and red are compression and tension, respectively, in the cement drawn as two lines at the bond periphery. Line thicknesses and orientations correspond with force magnitude and direction, respectively, and the moment in the cement contributes to the forces at the bond periphery.) rings in Fig. 6.2 (The existence of such large-scale features in the force chains provides evidence that, in general, the mechanism shown in Fig. 1(c) will be operative at a length scale larger than the grain size.) It is these microforces and micromoments that provide the local loading to produce grain and/or cement breakage, which, in turn, induces global force redistribution (because damaged material is softer and sheds load to stiffer, undamaged material) and the eventual formation of macroscopic fractures and/or rupture zones. As more and more grains and cement are broken, the material behaves more and more like a granular material with highly unstable force chains. The key to explaining material behavior is the evolving force-chain fabric, which is related in a complex way to the deformability and strength of the grains and the cement, the grain size, grain shape, grain packing and degree of cementation. Computational models of rock can be classified into two categories, depending on whether damage is represented indirectly, by its effect on constitutive relations or, directly, by the formation and tracking of many microcracks. Most indirect approaches idealize the material as a continuum and utilize average measures of material degradation in constitutive rela2 Figs. 5(b) and 6 are discussed in detail in Section 2.5, Materialgenesis procedure. 1331 tions to represent irreversible microstructural damage [18], while most direct approaches idealize the material as a collection of structural units (springs, beams, etc.) or separate particles bonded together at their contact points and utilize the breakage of individual structural units or bonds to represent damage [19–22]. Most computational models used to describe the mechanical behavior of rock for engineering purposes are based on the indirect approach, while those used to understand the behavior in terms of the progress of damage development and rupture are based on the direct approach. The BPM is an example of a direct modeling approach in which particles and bonds are related to similar objects observed microscopically in rock [23–30]. Alternative rock models in which the material is represented as a continuum include those in Refs. [31,32], where a network of weak planes is superimposed on an otherwise homogeneous elastic continuum, and those in Refs. [33,34], where the stiffness and strength of elements in a continuum with initial heterogeneous strength are allowed to degrade based on a strength failure criterion in the form of an elastic–brittle–plastic constitutive relation. These rock models exhibit more realistic responses in terms of shear and post-failure behavior than most lattice models, because they can carry compressive and shear forces arising from loading subsequent to bond breakage, whereas most lattice models retain no knowledge of the particles after each bond has broken. Comprehensive review articles [35,36] summarize rock models and approaches for simulating crack growth, and [37] provides additional discussion of the BPM and its relation to other rock models. The micromechanisms occurring in rock are complex and difficult to characterize within the framework of existing continuum theories. Microstructure controls many of these micromechanisms. The BPM approximates rock as a cemented granular material with an inherent length scale that is related to grain size and provides a synthetic material that can be used to test hypotheses about how the microstructure affects the macroscopic behavior. This is a comprehensive model that exhibits emergent properties (such as fracture toughness, which controls the formation of macroscopic fractures) that arise from a small set of microproperties for the grains and cement. The BPM does not impose theoretical assumptions and limitations on material behavior as do most indirect models (such as models that idealize the rock as an elastic continuum with many elliptical cracks—in the BPM, such cracks form, interact and coalesce into macroscopic fractures automatically). Behaviors not encompassed by current continuum theories can be investigated with the BPM. In fact, continuum behavior itself is simply another of the emergent properties of a BPM when averaged over the appropriate length scale. ARTICLE IN PRESS 1332 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 The remainder of the paper is organized as follows. The formulation of the BPM, which includes a microproperty characterization expressed in terms of grain and cement properties and a material-genesis procedure, is presented in the next section. In Sections 3 and 4, the measured macroscopic properties of a BPM for Lac du Bonnet granite are presented and discussed. In Section 3, the focus is on laboratory-scale behavior (during biaxial, triaxial and Brazilian tests), and Section 4 focuses on field-scale behavior (involving damage formation adjacent to an excavation). In both sections, the effect of particle size on model behavior is investigated. In Section 5, the emergent properties of the BPM are listed, and the general source of some of these behaviors is identified. As a particular example, a formal equivalence between the mechanisms and parameters of the BPM and the concepts and equations of linear elastic fracture mechanics (LEFM) is established and then used to explain the effect of particle size on model behavior. The capabilities and limitations of the BPM, as well as suggestions for overcoming these limitations, are presented in the conclusion. For completeness, the algorithms used to measure average stress and strain within the BPM, to install an initial stress field within the BPM, to create a densely packed BPM with low locked-in forces to serve as the initial material and to perform typical laboratory tests on the BPM are described in Appendix A. Throughout the paper, all vector and tensor quantities are expressed using indicial notation. 2. Formulation of the BPM The BPM simulates the mechanical behavior of a collection of non-uniform-sized circular or spherical rigid particles that may be bonded together at their contact points. The term ‘‘particle’’, as used here, differs from its more common definition in the field of mechanics, where it is taken as a body of negligible size that occupies only a single point in space; in the present context, the term ‘‘particle’’ denotes a body that occupies a finite amount of space. The rigid particles interact only at the soft contacts, which possess finite normal and shear stiffnesses. The mechanical behavior of this system is described by the movement of each particle and the force and moment acting at each contact. Newton’s laws of motion provide the fundamental relation between particle motion and the resultant forces and moments causing that motion. The following assumptions are inherent in the BPM: 1. The particles are circular or spherical rigid bodies with a finite mass. 2. The particles move independently of one another and can both translate and rotate. 3. The particles interact only at contacts; because the particles are circular or spherical, a contact is comprised of exactly two particles. 4. The particles are allowed to overlap one another, and all overlaps are small in relation to particle size such that contacts occur over a small region (i.e., at a point). 5. Bonds of finite stiffness can exist at contacts, and these bonds carry load and can break. The particles at a bonded contact need not overlap. 6. Generalized force–displacement laws at each contact relate relative particle motion to force and moment at the contact. The assumption of particle rigidity is reasonable when movements along interfaces account for most of the deformation in a material. The deformation of a packed-particle assembly, or a granular assembly such as sand, is described well by this assumption, because the deformation results primarily from the sliding and rotation of the particles as rigid bodies and the opening and interlocking at interfaces—not from individual particle deformation. The addition of bonds between the particles in the assembly corresponds with the addition of real cement between the grains of a sedimentary rock, such as sandstone, or notional cement between the grains of a crystalline rock, such as granite. The deformation of a bonded-particle assembly, or a rock, should be similar, and both systems should exhibit similar damage-formation processes under increasing load as the bonds are broken progressively and both systems gradually evolve toward a granular state. If individual grains or other microstructural features are represented as clusters of bonded particles, then grain crushing and material inhomogeneity at a scale larger than the grain size can also be accommodated by this model [38–40]. 2.1. Distinct-element method (DEM) The BPM is implemented in the two- and threedimensional discontinuum programs PFC2D and PFC3D [41] using the DEM. The DEM was introduced by Cundall [42] for the analysis of rock-mechanics problems and then applied to soils by Cundall and Strack [43]. Thorough descriptions of the method are given in Refs. [44,45]. The DEM is a particular implementation of a broader class of methods known as discrete-element methods, which are defined in Ref. [46] as methods that allow finite displacements and rotations of discrete bodies, including complete detachment, and recognize new contacts automatically as the calculation progresses. Current development and applications of discrete-element methods are described in Ref. [47]. In the DEM, the interaction of the particles is treated as a dynamic process with states of equilibrium ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 developing whenever the internal forces balance. The contact forces and displacements of a stressed assembly of particles are found by tracing the movements of the individual particles. Movements result from the propagation through the particle system of disturbances caused by wall and particle motion, externally applied forces and body forces. This is a dynamic process in which the speed of propagation depends on the physical properties of the discrete system. The calculations performed in the DEM alternate between the application of Newton’s second law to the particles and a force–displacement law at the contacts. Newton’s second law is used to determine the translational and rotational motion of each particle arising from the contact forces, applied forces and body forces acting upon it, while the force–displacement law is used to update the contact forces arising from the relative motion at each contact. The dynamic behavior is represented numerically by a time-stepping algorithm in which the velocities and accelerations are assumed to be constant within each time step. The solution scheme is identical to that used by the explicit finite-difference method for continuum analysis. The DEM is based on the idea that the time step chosen may be so small that, during a single time step, disturbances cannot propagate from any particle farther than its immediate neighbors. Then, at all times, the forces acting on any particle are determined exclusively by its interaction with the particles with which it is in contact. Because the speed at which a disturbance can propagate is a function of the physical properties of the discrete system (namely, the distribution of mass and stiffness), the time step can be chosen to satisfy the above constraint. The use of an explicit, as opposed to an implicit, numerical scheme provides the following advantages. Large populations of particles require only modest amounts of computer memory, because matrices are not stored. Also, physical instability may be modeled without numerical difficulty, because failure processes occur in a realistic manner—one need not invoke a non-physical algorithm, as is done in some implicit methods. 2.2. Damping of particle motion Because the DEM is a fully dynamic formulation, some form of damping is necessary to dissipate kinetic energy. In real materials, various microscopic processes such as internal friction and wave scattering dissipate kinetic energy. In the BPM, local non-viscous damping is used by specifying a damping coefficient, a: The damping formulation is similar to hysteretic damping [48], in which the energy loss per cycle is independent of the rate at which the cycle is executed. The damping force applied to each particle is given by F d ¼ ajF j signðV Þ; (1) 1333 where jF j is the magnitude of the unbalanced force on the particle and signðV Þ is the sign (positive or negative) of the particle velocity. This expression is applied separately to each degree-of-freedom. A common measure of attenuation or energy loss in rock is the seismic quality factor, Q. The quality factor is defined as 2p times the ratio of stored energy to dissipated energy in one wavelength Q ¼ 2pðW =DW Þ; (2) where W is energy [49]. For a single degree-of-freedom system, or for oscillation in a single mode, the quality factor can be written as [41] Q ¼ p=2a: (3) All BPMs described in this paper were run with a damping coefficient of 0.7, which corresponds with a quality factor of 2.2. The quality factor of Lac du Bonnet granite in situ is about 220 [50]; therefore, the models were heavily damped to approximate quasistatic conditions. Hazzard et al. [49] ran similar models using a damping coefficient corresponding with a quality factor of 100 and found that the energy released by crack events (i.e., bond breakages) may have a significant influence on the rock behavior, because the waves emanating from cracks are capable of inducing more cracks if nearby bonds are close to failure. Hazzard and coworkers showed that the effect of this dynamically induced cracking was to decrease the unconfined compressive strength by up to 15%. In addition, if low numerical damping is used, then seismic source information (such as magnitude and mechanism) can be determined for every modeled crack [51]. 2.3. Grain-cement behavior and parameters The BPM mimics the mechanical behavior of a collection of grains joined by cement. The following discussion considers each grain as a PFC2D or PFC3D particle and each cement entity as a parallel bond. The total force and moment acting at each cemented contact is comprised of a force, F i ; arising from particle–particle overlap, denoted in Fig. 3 as the grain behavior, and a force and moment, F̄ i and M̄ i ; carried by the parallel bond and denoted as the cement behavior. These quantities contribute to the resultant force and moment acting on the two contacting particles. The force–displacement law for this system will now be described, first for the grain behavior and then for the cement behavior. Note that if a parallel bond is not present at a contact, then only the grain-based portion of the force–displacement behavior occurs. The grain-based portion of the force–displacement behavior at each contact is described by the following six parameters (see Fig. 3): the normal and shear stiffnesses, kn and ks ; and the friction coefficient, m; of ARTICLE IN PRESS 1334 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Fig. 3. Force–displacement behavior of grain-cement system. the two contacting particles, which are assumed to be disks of unit thickness in PFC2D and spheres in PFC3D. Whenever two particles overlap, a contact is formed at the center of the overlap region along the line joining the particle centers (xðcÞ i in Fig. 3), and two linear springs are inserted (with stiffnesses derived from the particle stiffnesses assuming that the two particle springs act in series) along with a slider in the shear direction. The contact force vector, F i (which represents the action of particle A on particle B), can be resolved into normal and shear components with respect to the contact plane as F i ¼ F n ni þ F s ti ; n (4) s where F and F denote the normal and shear force components, respectively, and ni and ti are the unit vectors that define the contact plane. The normal force is calculated by F n ¼ K nU n; (5) where U n is the overlap and K n is the contact normal stiffness given by Kn ¼ ðBÞ kðAÞ n kn ; ðBÞ kðAÞ n þ kn ðBÞ kðAÞ n and kn (6) with being the particle normal stiffnesses. The shear force is computed in an incremental fashion. When the contact is formed, F s is initialized to zero. Each subsequent relative shear–displacement increment, DU s ; produces an increment of elastic shear force, DF s ; that is added to F s (after F s has been rotated to account for motion of the contact plane). The increment of elastic shear force is given by DF s ¼ ks DU s ; (7) where ks is the contact shear stiffness given by ks ¼ ðBÞ kðAÞ s ks ðBÞ kðAÞ s þ ks ; (8) with kðAÞ and kðBÞ being the particle shear stiffnesses. s s The relative displacement increment during the time step Dt is given by DU i ¼ V i Dt; where V i is the contact velocity

 V i ¼ x_ ðcÞ  x_ ðcÞ i i B

A  ðBÞ ðBÞ ðBÞ ¼ x_ i þ eijk oj xðcÞ  x k k

 ðAÞ ðAÞ ðcÞ  x_ i þ eijk oj xk  xðAÞ ; ð9Þ k x_ i and oj are the particle translational and rotational velocities, respectively, and eijk is the permutation symbol. The relative shear–displacement increment vector is DU si ¼ V si Dt ¼ ðV i  V ni ÞDt ¼ ðV i  V j nj ni ÞDt: (10) If U n p0 (a gap exists), then both normal and shear forces are set to zero; otherwise, slip is accommodated by computing the contact friction coefficient   m ¼ min mðAÞ ; mðBÞ ; (11) with mðAÞ and mðBÞ being the particle friction coefficients, and setting F s ¼ mF n if F s 4mF n : ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Note that K n is a secant stiffness in that it relates total displacement and force, whereas ks is a tangent stiffness in that it relates incremental displacement and force. In this paper, an upper-case K denotes a secant stiffness and a lower-case k denotes a tangent stiffness. Use of a secant stiffness to compute normal force makes the computations less prone to numerical drift and able to handle arbitrary placement of particles and changes in particle radii after a simulation has begun. The cement-based portion of the force–displacement behavior at each cemented contact is described by the following five parameters that define a parallel bond (see Fig. 3): normal and shear stiffnesses per unit area, n s k̄ and k̄ ; tensile and shear strengths, s̄c and t̄c ; and bond-radius multiplier, l̄; such that parallel-bond radius   R̄ ¼ l̄ min RðAÞ ; RðBÞ ; (12) with RðAÞ and RðBÞ being the particle radii. A parallel bond approximates the mechanical behavior of a brittle elastic cement joining the two bonded particles. Parallel bonds establish an elastic interaction between these particles that acts in parallel with the grain-based portion of the force–displacement behavior; thus, the existence of a parallel bond does not prevent slip. Parallel bonds can transmit both force and moment between particles, while grains can only transmit force. A parallel bond can be envisioned as a set of elastic springs uniformly distributed over a rectangular crosssection in PFC2D and a circular cross-section in PFC3D lying on the contact plane and centered at the contact point. These springs behave as a beam whose length, L̄ in Fig. 3, approaches zero to approximate the mechanical behavior of a joint. The total force and moment carried by the parallel bond are denoted by F̄ i and M̄ i ; respectively, which represent the action of the bond on particle B. The force and moment vectors can be resolved into normal and shear components with respect to the contact plane as n s F̄ i ¼ F̄ ni þ F̄ ti ; n s M̄ i ¼ M̄ ni þ M̄ ti ; n s ð13Þ n s where F̄ ; F̄ and M̄ ; M̄ denote the axial- and sheardirected forces and moments, respectively, and ni and ti are the unit vectors that define the contact plane. n (For the PFC2D model, the twisting moment, M̄ ¼ 0; s and the bending moment, M̄ ; acts in the out-of-plane direction.) When the parallel bond is formed, F̄ i and M̄ i are initialized to zero. Each subsequent relative displacement- and ðDU n ; DU s ; Dyn ;
rotation increment  ðBÞ ðAÞ s Dy with Dyi ¼ oi  oi Dt produces an increment of elastic force and moment that is added to the current values (after the shear-component vectors have been rotated to account for motion of the contact plane). The increments of elastic force and moment 1335 are given by n n DF̄ ¼ k̄ A DU n ; s s DF̄ ¼ k̄ A DU s ; n s s n DM̄ ¼ k̄ J Dyn ; DM̄ ¼ k̄ I Dys ; ð14Þ where A; I and J are the area, moment of inertia and polar moment of inertia of the parallel bond crosssection, respectively. These quantities are given by ( 2R̄t; t ¼ 1; PFC2D; A¼ 2 pR̄ ; PFC3D; 8 3 2 < R̄ t; t ¼ 1; PFC2D; 3 I¼ :1 4 pR̄ ; PFC3D; (4 NA; PFC2D; J¼ ð15Þ 4 1 PFC3D: 2 pR̄ ; The maximum tensile and shear stresses acting on the parallel-bond periphery are calculated from beam theory to be n s F̄ jM̄ jR̄ ; þ ¼ s̄ I A s n jF̄ j jM̄ jR̄ þ : ð16Þ t̄max ¼ A J If the maximum tensile stress exceeds the tensile strength ðs̄max Xs̄c Þ or the maximum shear stress exceeds the shear strength ðt̄max Xt̄c Þ; then the parallel bond breaks, and it is removed from the model along with its accompanying force, moment and stiffnesses. A simplified form of the BPM represents the cement using contact bonds instead of parallel bonds. A contact bond approximates the physical behavior of a vanishingly small cement-like substance lying between and joining the two bonded particles. The contact bond behaves, essentially, as a parallel bond of radius zero. Thus, a contact bond does not have a radius or shear and normal stiffnesses, as does a parallel bond, and cannot resist a bending moment or oppose rolling; rather, it can only resist a force acting at the contact point. Also, slip is not allowed to occur when a contact bond is present. A contact bond is defined by the two parameters of tensile and shear strengths, fn and fs ; expressed in force units. When the corresponding component of the contact force exceeds either of these values, the contact bond breaks, and only the grainbased portion of the force–displacement behavior occurs. max 2.4. Microproperty characterization In general, the BPM is characterized by the grain density, grain shape, grain size distribution, grain ARTICLE IN PRESS 1336 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 packing and grain-cement microproperties. Each of these items influences the model behavior. The grain density, r; does not affect the quasi-static behavior but is included for completeness. In this paper, we focus on circular or spherical grains comprised of individual PFC2D or PFC3D particles, although the effect on the strength envelope of introducing particle clusters of complex interlocking shapes into the PFC2D particle assembly is discussed in Section 3.4. The particle diameters satisfy a uniform particle-size distribution bounded by Dmin and Dmax ; and a dense packing is obtained using the material-genesis procedure of Section 2.5. The packing fabric, or connectivity of the bonded assembly, is controlled by the ratio ðDmax =Dmin Þ; for a fixed ratio, varying Dmin changes the absolute particle size but does not affect the packing fabric. Such a microproperty characterization separates the effects of packing fabric and particle size on material behavior and clearly identifies Dmin as the controlling length scale of the material. The final item that characterizes the BPM is the grain-cement microproperties:   E c ; kn =ks ; m ; grain microproperties;  n s l̄; Ē c ; k̄ =k̄ ; s̄c ; t̄c ; cement microproperties; ð17Þ where E c and Ē c are the Young’s moduli of the grains n s and cement, respectively; ðkn =ks Þ and ðk̄ =k̄ Þ are the ratios of normal to shear stiffness of the grains and cement, respectively; l̄ is the radius multiplier used to set the parallel-bond radii via Eq. (12); m is the grain friction coefficient; and s̄c and t̄c are the tensile and shear strengths, respectively, of the cement. In the analysis below, the grain and cement moduli are related to the corresponding normal stiffnesses such that the particle and parallel-bond stiffnesses are assigned as ( 2tE c ; t ¼ 1; PFC2D; kn :¼ 4RE c ; PFC3D; kn ks :¼ ; ðkn =ks Þ Ē c n ; k̄ :¼ ðAÞ R þ RðBÞ n k̄ s k̄ :¼ n s ; ð18Þ ðk̄ =k̄ Þ where R is particle radius. The usefulness of these modulus–stiffness scaling relations is confirmed in Section 3.5, where it is shown that the macroscopic elastic constants are independent of particle size for the PFC2D material and exhibit only a minor size effect for the PFC3D material. The deformability of an isotropic linear elastic material is described by two elastic constants. These quantities are emergent properties of the BPM and cannot be specified directly. It is possible, however, to relate the grain and cement moduli, E c and Ē c ; Fig. 4. Equivalent continuum material of grain-cement system. respectively, to the normal stiffnesses by envisioning the material at each contact as an elastic beam with its ends at the particle centers, as shown in Fig. 4. The axial stiffness of such a beam is K ¼ AE=L; where A, E and L are the cross-sectional area, modulus and length, respectively, of the beam. For the grain-based behavior, kn ðLtÞE c kn ¼ ¼ E c t ) E c ¼ ; t ¼ 1; PFC2D; 2 L 2t kn ðL2 ÞE c kn kn ¼ ¼ EcL ) Ec ¼ ¼ ; PFC3D ð19Þ 2 L 2L 4R by assuming that kn ¼ knðAÞ ¼ kðBÞ in Eq. (6). For the n cement-based behavior, n AĒ c AĒ c ¼ ðAÞ L R þ RðBÞ  n ) Ē c ¼ k̄ RðAÞ þ RðBÞ : k̄ A ¼ ð20Þ The cement modulus is dependent on particle size; to achieve a constant cement modulus, parallel-bond stiffnesses must be scaled with the particle radii. For the PFC3D material, the grain modulus is dependent on particle size; to achieve a constant grain modulus, particle stiffnesses must be scaled with particle radius. This analysis does not yield a relation between Poisson’s ratio and particle stiffness at the microlevel; however, a macroscopic Poisson’s ratio will be observed, and its value will be affected by grain shape, grain packing and n s the ratios ðkn =ks Þ and ðk̄ =k̄ Þ: For a fixed grain shape and packing, increasing these ratios increases the Poisson’s ratio. 2.5. Material-genesis procedure The BPM represents rock as a dense packing of nonuniform-sized circular or spherical particles that are joined at their contact points with parallel bonds. The material-genesis procedure produces this system such that the particles are well connected and the locked-in forces are low. Both the packing and the locked-in ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1337 forces are arbitrary and isotropic at a macroscale—i.e., when averaged over approximately one dozen adjacent particles. (This would not occur if the material were created by gravity compaction—in which case, the force chains would be aligned vertically and their magnitudes would increase with depth.) The material-genesis procedure employs the following five-step process (illustrated in Figs. 5 and 6 for the PFC2D Lac du Bonnet granite model of Table 1). 1. Compact initial assembly: A material vessel consisting of planar frictionless walls is created, and an assembly of arbitrarily placed particles is generated to fill the vessel. The vessel is a rectangle bounded by four walls for PFC2D and a rectangular parallelepiped bounded by six walls for PFC3D. The wall normal stiffnesses are made just larger than the average particle normal stiffness to ensure that the particle–wall overlap remains small. The particle diameters satisfy a uniform particle-size distribution bounded by Dmin and Dmax : To ensure a reasonably tight initial packing, the number of particles is determined such that the overall porosity in the Fig. 6. Force and moment distributions in PFC2D material after removal from material-genesis vessel followed by relaxation (color convention described in Fig. 2). Fig. 5. Material-genesis procedure for PFC2D model: (a) particle assembly after initial generation but before rearrangement; (b) contact-force distribution (All forces are compressive, thickness proportional to force magnitude.) after step 2; (c) floating particles (with less than three contacts) and contacts after step 2; (d) parallel-bond network after step 4. ARTICLE IN PRESS 1338 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Table 1 Model microproperties for Lac du Bonnet granite Grains r ¼ 2630 kg=m3 ðDmax =Dmin Þ ¼ 1:66; Dmin varies 62 GPa; PFC2D Ec ¼ 72 GPa; PFC3D ðkn =ks Þ ¼ 2:5 m = 0.5 vessel is 16% for PFC2D and 35% for PFC3D. The particles, at half their final size, are placed randomly such that no two particles overlap. Then, the particle radii are increased to their final values, as shown in Fig. 5(a), and the system is allowed to rearrange under zero friction. The ballgeneration procedure is described in more detail in Appendix A. 2. Install specified isotropic stress: The radii of all particles are reduced uniformly to achieve a specified isotropic stress, s0 ; defined as the average of the direct stresses. These stresses are measured by dividing the average of the total force acting on opposing walls by the area of the corresponding specimen cross-section. Stresses in the PFC2D models are computed assuming that each particle is a disk of unit thickness. When constructing a granite material, s0 is set equal to approximately 1% of the uniaxial compressive strength. This is done to reduce the magnitude of the locked-in forces that will develop after the parallel bonds are added and the specimen is removed from the material vessel and allowed to relax, as shown in Fig. 6. The magnitude of the locked-in forces (both tensile and compressive) will be comparable to the magnitude of the compressive forces at the time of bond installation. The contact-force distribution at the end of this step is shown in Fig. 5(b). The isotropic stress installation procedure is described in more detail in the Appendix A. 3. Reduce the number of ‘‘floating’’ particles: An assembly of non-uniform-sized circular or spherical particles, placed randomly and compacted mechanically, can contain a large number (perhaps as high as 15%) of ‘‘floating’’ particles that have less than N f contacts, as shown in Fig. 5(c), where N f ¼ 3: It is desirable to reduce the number of floating particles such that a denser bond network is obtained in step 4. In our granite models, we wish to obtain a densely packed and well-connected assembly to mimic a highly interlocked collection of grains. By setting N f ¼ 3 and the allowed number of floating particles Cement l̄ ¼ 1 62 GPa; 72 GPa; n s ðk̄ =k̄ Þ ¼ 2:5 Ē c ¼ PFC2D PFC3D s̄c ¼ t̄c ¼ ðmean  std: dev:Þ¼ 157  36 MPa; 175  40 MPa; PFC2D PFC3D to zero, we obtain a bonded assembly for which nearly all particles away from the specimen boundaries have at least three contacts. The floaterelimination procedure is described in more detail in Appendix A. 4. Install parallel bonds: Parallel bonds are installed throughout the assembly between all particles that are in near proximity (with a separation less than 106 times the mean radius of the two particles), as shown in Fig. 5(d). The parallel-bond properties are assigned to satisfy Eqs. (17) and (18), with s̄c and t̄c picked from a Gaussian (normal) distribution. The grain property of m is also assigned. 5. Remove from material vessel: The material-genesis procedure is completed by removing the specimen from the material vessel and allowing the assembly to relax. This is done by deleting the vessel walls and stepping until static equilibrium is achieved. During the relaxation process, the material expands and generates a set of self-equilibrating locked-in forces, as shown in Fig. 6. These are similar to locked-in stresses that may exist in a free specimen of rock. The locked-in forces can be divided into two types, depending upon their existence at a macro- or a microscale. The macroscale forces exist within selfcontained clusters of approximately one dozen particles, in which the inner particles are in tension and the outer particles are in compression, or vice versa. Such compression rings can be seen in Fig. 6. The microscale forces exist within individual parallel bonds and are mostly tension to equilibrate the tendency of the contact springs to repel the compressed particles. There is a net energy stored in the assembly in the form of strain energy in both the particle–particle contacts and the parallel bonds. These locked-in forces can have a significant effect upon behavior (e.g., feeding energy into a rockburst or causing strain to occur when the assembly is cut). Holt [52] used a BPM to investigate stress-relief effects induced during coring of sandstone by setting s0 equal to the compaction stresses existing at the time of grain cementation. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 3. Measured macroscopic-properties (laboratory scale) The BPM consists of a procedure for material genesis and two independent sets of microproperties for shortand long-term behavior. A set of short-term microproperties are described in Section 2.4, and a set of longterm microproperties are described in Ref. [53]. No model is complete or fully verifiable [54], but the validity of the BPM is demonstrated by comparing model behavior with measured and observed responses of Lac du Bonnet granite at both laboratory and field scales in this and the next section. 3.1. Choosing microproperties For continuum models, the input properties (such as modulus and strength) can be derived directly from measurements performed on laboratory specimens. For the BPM, which synthesizes macroscale material behavior from the interactions of microscale components, the input properties of the components usually are not known. The relation between model parameters and commonly measured material properties is only known a priori for simple packing arrangements. For the general case of arbitrary packing of arbitrarily sized particles, the relation is found by means of a calibration process in which a particular instance of a BPM—with a particular packing arrangement and set of model parameters—is used to simulate a set of material tests, and the BPM parameters then are chosen to reproduce the relevant material properties measured in such tests. The BPM for Lac du Bonnet granite is described by the microproperties in Table 1. These microproperties were chosen to match the macroproperties of Lac du Bonnet granite discussed in Section 3.3. The ratio ðDmax =Dmin Þ is set greater than 1 to produce an arbitrary isotropic packing. As ðDmax =Dmin Þ ! 1; the packing tends toward a crystalline arrangement. Such a uniform packing exhibits anisotropic macroproperties and is not representative of the more complex grain connectivity of a granite. l̄ is set equal to 1 to produce a material with cement that completely fills the throat between cemented particles. As l̄ ! 0; the material behavior approaches that of a granular material. The grain and cement moduli and ratios of normal to shear stiffness are set equal to one another to reduce the number of free parameters. Then, the moduli are chosen to match the Young’s modulus, and the ratios of normal to shear stiffness are chosen to match the Poisson’s ratio. Next, the cement strengths are set equal to one another so as not to exclude mechanisms that may only be activated by microshear failure (see below). The ratio of standard deviation to mean of the cement strengths is chosen to match the crack-initiation stress (defined in Section 3.3), and the mean value of the cement strengths is chosen to match the unconfined compressive strength. The parti- 1339 cle-friction coefficient appears to affect only post-peak response, and it is not clear to what it should be calibrated; thus, m ¼ 0:5 is used as a reasonable non-zero value. By setting s̄c ¼ t̄c ; both tensile and shear microfailures are possible. If microtensile failure is excluded (by setting s̄c to infinity), then cracking does not localize onto a single macrofracture plane under macroscopic extensile loading. Because this mechanism does occur in granite, microtensile failure is allowed to occur in the granite model. The alternative cases for which microshear failure is excluded fully (by setting t̄c to infinity) or allowed (by setting s̄c ¼ t̄c ) are investigated for the PFC2D material by Potyondy and Autio [39], who study damage formation adjacent to a circular hole subjected to far-field compressive loading. Both materials produce breakout notches in compressive regions and tensile macrofractures in tensile regions that are generally similar. The material with s̄c ¼ t̄c allows more well-defined notches to form than does the material that can only fail in the microtensile mode (see Fig. 7). The former material accommodates the shearing deformation along the notch outline by forming a string of shear microcracks, whereas the latter material must form several sets of en echelon tensile microcracks (which requires that additional deformation and accompanying damage occur within the notch region). For the granite model, s̄c ¼ t̄c so as not to exclude mechanisms that may only be activated via microshear failure. 3.2. PFC2D model behavior during biaxial and Brazilian tests Typical stress–strain response and damage patterns are shown in Fig. 8 for the PFC2D granite model with an average particle diameter of 0.72 mm. The testing procedures are described in Appendix A. The crack distributions in this and all such figures are depicted as bi-colored lines (in which red represents tension-induced parallel-bond failure for which bond tensile strength has been exceeded and blue represents shear-induced parallel-bond failure for which bond shear strength has been exceeded) lying between the two previously bonded particles with a length equal to the average diameter of the two previously bonded particles. The cracks are oriented perpendicular to the line joining the centers of the two previously bonded particles. The damage mechanisms during biaxial tests are summarized in observations 1–3, and those during Brazilian tests in observation 4. Similar observations apply to the PFC3D granite model subjected to triaxial and Brazilian tests. 1. During the early stages of biaxial loading, relatively few cracks form, and those that do form tend to be aligned with the direction of maximum compression. These cracks are distributed throughout the specimen ARTICLE IN PRESS 1340 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Fig. 7. Damage in the notch region for the same compressive load (acting vertically) for a PFC2D material with (a) microshear failure allowed (3267 cracks) and (b) microshear failure excluded (3800 cracks). Fig. 8. Stress–strain response and damage patterns formed during biaxial tests at confinements of 0.1, 10 and 70 MPa for a 63.4 126.8 mm2 specimen of the PFC2D granite model with Davg ¼ 0:72 mm: (a) Stress-strain response; (b) Post-peak crack distribution, s3 ¼ 0:1 MPa (1716 cracks, tensile/shear=1452/264); (c) Post-peak crack distribution, s3 ¼ 10 MPa (2443 cracks, tensile/shear=2035/408); (d) Post-peak crack distribution, s3 ¼ 70 MPa (3940 cracks, tensile/shear=2800/1140). ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 and do not seem to be interacting. The onset of this early cracking is quantified by the crack-initiation stress. 2. Near peak load in a biaxial test, one or more macroscopic failure planes develop that cut across the specimen. When confinement is low, a set of secondary macrofractures oriented parallel with the loading direction form on either side of these failure planes, as seen in Fig. 8(b). The secondary macrofractures are comprised of tensile microcracks and suggest an axial-splitting phenomenon. As the confinement increases, the number of failure planes and their widths increase, as seen in Figs. 8(c) and (d). 3. Damage at peak load is quite similar for all confinements, and most damage occurs after peak load as the specimen expands laterally—volumetric strain is dilational [53]. This observation, combined with observation 2, suggests that the post-peak damage-formation process is affected more by confinement than is the pre-peak damage-formation process. This change in behavior is related in part to the fact that, as confinement increases, the ratio of shear microcracks to tensile microcracks increases. The confinement reduces the tensile forces that develop in a direction perpendicular to the specimen axis and thereby causes more shear microcracks to form. 4. During Brazilian tests, a wedge of cracks forms inside of the specimen below each platen. Then, one of the wedges initiates a single macrofracture that travels across the specimen parallel with the direction of loading. The macrofracture is comprised of tensile microcracks and is driven by extensile deformation across the crack path (similar to an LEFM mode-I crack). This observation suggests that the Brazilian strength measured in these tests is related to the material fracture toughness, as is shown in Section 5.2, and that it may not correspond with the Brazilian strength measured in a valid Brazilian test on a real rock specimen. Such difficulties are inherent in tests that attempt to measure the tensile strength of a rock, because it is rarely possible to force the failure to 1341 occur simultaneously across the entire final fracture plane. 3.3. Macroproperties of Lac du Bonnet granite The following short-term laboratory properties of Lac du Bonnet granite [55] obtained from 63-mm diameter specimens with a length-to-diameter ratio of 2.5 are used to calibrate the short-term response of the granite model (see Table 2 for mean, standard deviation and number of tests): (a) the elastic constants of Young’s modulus, E; and Poisson’s ratio, u; measured from unconfined compression tests; (b) the unconfined compressive strength, qu ; (c) the crack-initiation stress, sci ; (d) the strength envelope for confining pressures in the range 0.1–10 MPa approximated as an equivalent Mohr–Coulomb material with a secant slope, N f ; friction angle, f; and cohesion, c; and (e) the Brazilian strength, st : The crack-initiation stress is defined as the axial stress at which non-elastic dilation just begins, identified as the point of deviation from linear elasticity on a plot of axial stress versus volumetric strain. The strength envelope, which is the relation between peak strength and confining pressure, is obtained by fitting the Hoek–Brown strength criterion [56] to results of triaxial tests at confining pressures up to 70 MPa. The Hoek–Brown strength criterion is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sf ¼ s3  ss2c  s3 sc m (21) where s and m are dimensionless material parameters and sc is equal to the unconfined compressive strength when s ¼ 1: The Hoek–Brown parameters for the strength envelope of Lac du Bonnet granite [57] are sc ¼ 213 MPa; m ¼ 31 and s ¼ 1: In general, the strength envelope will exhibit a decreasing slope for increasing confinement; however, it can be linearized using a secant approximation defined by the strength, sf ; at two confinements, P0 and P1 : If the slope is defined by Nf ¼ sf ðP1 Þ  sf ðP0 Þ ; P1  P0 (22) Table 2 Macroproperties of the models and Lac du Bonnet granite Property Lac du Bonnet granite PFC2D model, Davg ¼ 0:72 mm PFC3D model, Davg ¼ 1:53 mm E ðGPaÞ u qu ðMPaÞ sci ðMPaÞ scd ðMPaÞ Nf f ðdegÞ c ðMPaÞ st ðMPaÞ 69  5:8 ðn ¼ 81Þ 0:26  0:04 ðn ¼ 81Þ 200  22 ðn ¼ 81Þ 90 þ s3 ; s3 o30 150; s3 ¼ 0 13 59 30 9:3  1:3 ðn ¼ 39Þ 70:9  0:9 ðn ¼ 10Þ 0:237  0:011 ðn ¼ 10Þ 199:1  13:0 ðn ¼ 10Þ 71:8  21:8 ðn ¼ 10; s3 ¼ 0:1Þ NA 3:01  0:60 ðn ¼ 10Þ 29:5  4:8 ðn ¼ 10Þ 58:5  8:5 ðn ¼ 10Þ 44:7  3:3 ðn ¼ 10Þ 69:2  0:8 ðn ¼ 10Þ 0:256  0:014 ðn ¼ 10Þ 198:8  7:2 ðn ¼ 10Þ 86:6  11:0 ðn ¼ 10; s3 ¼ 0:1Þ 190:3  7:5 ðn ¼ 10; s3 ¼ 0:1Þ 3:28  0:33 ðn ¼ 10Þ 32:1  2:4 ðn ¼ 10Þ 55:1  4:2 ðn ¼ 10Þ 27:8  3:8 ðn ¼ 10Þ ARTICLE IN PRESS 1342 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 then the friction angle and cohesion can be written as [58]   1 N f  1 f ¼ sin ; Nf þ 1 q c ¼ puffiffiffiffiffiffiffi ; ð23Þ 2 Nf where qu is the unconfined compressive strength. The Brazilian strength is computed via Ff ; (24) pRt where F f is the peak force acting on the platens, and R and t are the radius and thickness, respectively, of the Brazilian disk. The crack-damage stress, scd ; is defined as the axial stress at which the volumetric strain reverses direction and becomes dilational. The crack-damage stress for saturated Lac du Bonnet granite is related to the peak strength, sf ; by the following fourth-order polynomial [53]: scd ¼ 0:75  0:031x þ 0:00285x2 sf  0:000104x3 þ 0:00000138x4 ; ð25Þ st ¼ Fig. 9. Typical specimens of the granite model: (a) PFC2D specimen (4017 particles, 7764 parallel bonds) and (b) PFC3D specimen (19,749 particles, 55,042 parallel bonds). where x is the confining stress (in MPa). The crackdamage stress is measured for the PFC3D material, but it is not used in the present calibration. 3.4. Macroproperties of the PFC models The microproperties in Table 1 with Dmin ¼ 0:55 mm and Dmin ¼ 1:19 mm were used to produce PFC2D and PFC3D materials with average particle diameters of 0.72 and 1.53 mm, respectively. Then, 10 rectangular PFC2D specimens (63.4 31.7 mm2) and 10 rectangular parallelepiped PFC3D specimens (63.4 31.7 31.7 mm3) with different packing arrangements and microstrengths were created by varying the seed of the random-number generator during material genesis. Typical PFC2D/ PFC3D specimens (shown in Fig. 9) contain approximately 4000/20,000 particles and have 44/20 particles across the specimen width. Biaxial, triaxial and Brazilian tests were performed upon these specimens to obtain the model macroproperties shown in Table 2. The testing procedures are described in Appendix A. The strength envelopes were obtained by performing a set of biaxial and triaxial tests at confining pressures of 0.1, 1, 5, 10, 20, 30 and 70 MPa. The test results for the PFC2D material are shown in Fig. 10, which includes the peak strength and crack-initiation stress from each test with lines joining the average values at each confining pressure. The strength envelope for the PFC3D material is similar. The PFC models match the macroproperties of E; n and qu of the Lac du Bonnet granite, and the variability Fig 10. Strength envelopes for PFC2D granite model (results of all 10 tests) and Lac du Bonnet granite. of these macroproperties (expressed as a ratio of standard deviation to mean) is less than that of the physical material. The unconfined compressive strength of the PFC3D material exhibits less variability than that of the PFC2D material. The model materials also exhibit the onset of significant internal cracking (measured by sci ) before ultimate failure. However, the model strengths only match that of the Lac du Bonnet granite for stresses near the uniaxial state. The slope of the strength envelope is too low (3 for both PFC2D and PFC3D versus 13 for the granite), and the Brazilian strength is too high (the ratio ðqu =st Þ is 4.5 and 7.2 for PFC2D and PFC3D, respectively, versus 21.5 for the granite). This discrepancy may arise from the use of circular and spherical grains in the present model, and it could be reduced by using grain shapes that more closely resemble the complex-shaped and highly interlocked crystalline grains in granite. A preliminary investigation, described below, suggests that introducing finitestrength particle clusters of complex interlocking shapes into the PFC2D particle assembly will increase the slope ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 of the strength envelope and will lower the crackdamage stress. A similar effect is expected for the PFC3D material. If these grains are also allowed to break, then damage mechanisms that more closely resemble those in granite will be possible. One mechanism that contributes to the large difference in the compressive and tensile strengths of granite is suggested by Mosher et al. [59], who have studied the cracking within granite thin sections when stressed in tension or compression under the microscope and have found that intergranular fractures predominate in tension and intragranular in compression. In the preliminary investigation, unbreakable particle clusters of complex interlocking shapes were introduced into the PFC2D particle assembly so that cracking is forced to occur along cluster boundaries (the white dots in Fig. 11(a)). The clustering algorithm is controlled by S c ; the maximum number of particles in a cluster. A cluster is defined as a set of particles that are adjacent to 1343 one another, where adjacent means that a connected path can be constructed between any two particles in a cluster by traversing bonded particle–particle contacts. The algorithm begins with a densely packed bonded material and identifies particle clusters by traversing the list of particles. Each cluster is grown by identifying the current particle as a seed particle and then adding adjacent particles to the cluster until either all adjacent particles have been added or the cluster size has reached S c : The algorithm provides no control over cluster shape but does produce a collection of complex cluster shapes that are fully interlocking—similar to the interlocking grains in a crystalline rock. The strength-envelope slope for such a material can be made to exceed that of Lac du Bonnet granite, as shown in Fig. 11(b), and dilation begins at lower stresses relative to peak (scd is lowered) as cluster size is increased. However, the biaxial-test damage does not localize into discrete macroscopic failure planes, and the post-peak response is plastic-like and does not exhibit the abrupt stress drop of Lac du Bonnet granite. Experiments on Barre granite [4] indicate that, after loading to approximately 87% of peak strength, almost all grain boundaries were cracked along their entire length. Because this had occurred before the peak strength had been reached, it suggests that additional damage, in the form of transgranular fracture (or grain splitting) may be occurring in this loading regime. If a finite strength were assigned to the bonds within clusters, then transgranular fractures could form, and this might produce the desired localization and abrupt post-peak stress drop. This presents a promising avenue for further study [60]. 3.5. Effect of particle size on macroproperties Fig. 11. Introducing complex grain-like shapes into the PFC2D material using particle clusters: (a) clusters of size 7 and bonds (black: intra-cluster bonds; white: inter-cluster bonds) and (b) effect of cluster size on strength envelope for unbreakable clusters. If the BPM is used to predict damage formation surrounding an excavation, then the absolute size of the potential damage region is fixed, and successively finer resolution models can be produced by decreasing Dmin while keeping all other microproperties fixed. The investigation in this section demonstrates that particle size is an intrinsic part of the material characterization that affects the Brazilian strength (and the unconfined compressive strength for PFC3D material); thus, particle size cannot be regarded as a free parameter that only controls model resolution. If the damage mechanisms involve extensile conditions similar to those existing in a Brazilian test, then particle size should be chosen to match the Brazilian strength or another appropriate property of the rock, perhaps the fracture toughness. For the PFC3D material, it may not be possible to match the ratio ðqu =st Þ without introducing nonspherical grains. Four PFC2D and four PFC3D granite materials were produced using the microproperties in Table 1 with ARTICLE IN PRESS 1344 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 different values of Dmin : The PFC2D and PFC3D materials have average particle diameters ranging from 2.87 to 0.36 mm and 5.95 to 1.53 mm, respectively. For each PFC2D and PFC3D material, 10 rectangular specimens (63.4 31.7 mm2) or 10 rectangular parallelepiped specimens (63.4 31.7 31.7 mm3), respectively, with different packing arrangements and microstrengths were created by varying the seed of the random-number generator during material genesis. Each specimen was then subjected to a Brazilian test and to two biaxial or triaxial tests at confinements of 0.1 and 10 MPa. The PFC3D Brazilian disk thicknesses were chosen to produce approximately 2.6 particles across the thickness. The test results are presented in Tables 3 and 4 in terms of the mean and coefficient of variation (ratio of standard deviation to mean) of each macroproperty. It is expected that as particle size continues to decrease, the coefficients of variation will converge to specific values. These values should be a true measure of the effect of both packing and strength heterogeneities in the model materials, because the number of particles across the specimen has become large enough to obtain a representative sampling of the material response. As the number of particles across the specimen is reduced, the scatter in measured properties increases. Also, a sufficient number of particles must be present for the model to resolve and reproduce the failure mechanisms that influence the strength properties. 3.5.1. PFC2D material response The elastic constants appear to be independent of particle size. The mean values of Poisson’s ratio are approximately the same for all particle sizes, and the mean values of Young’s modulus increase slightly (by less than 5%) as particle size decreases from 2.87 to 0.36 mm. This slight increase in Young’s modulus may be an artifact of having sampled too few data points to obtain a true measure of the mean values. If the elastic constants truly are independent of particle size, then the characterization of the elastic grain-cement n s microproperties in terms of E c ; Ē c ; ðkn =ks Þ; and ðk̄ =k̄ Þ provides a size-independent means of specifying the elastic properties of the material. The size independence is achieved by scaling the parallel-bond stiffnesses as a function of particle size via Eq. (18). The unconfined compressive strength also appears to be independent of particle size. The mean values differ by less than 8% and exhibit no clear increasing or decreasing trend. This behavior may be an artifact of having sampled too few data points to obtain a true measure of the mean values; however, the decreasing coefficients of variation are reasonable and indicate a possible convergence to a specific value. It is also reasonable that the variability of qu is greater than the variability of the elastic constants, because qu measures a complex critical-state phenomenon involving extensive damage formation and interaction, whereas the elastic constants merely measure the deformability of the particle assembly before any significant damage has developed. No definitive statements about the effect of particle size on friction angle and cohesion can be made based on the results in Table 3. The mean values of f and c differ by 19% and 21%, respectively, and exhibit no Table 3 Effect of particle size on PFC2D macroproperties Particle size Davg (mm) 2.87 1.44 0.72 0.36 Macroproperty (63.4 31.7 mm2 specimens, n=10, mean and coefficient of variation) E (GPa, %) n (, %) 68.3 68.8 70.9 71.5 0.231 0.249 0.237 0.245 10.5 3.3 1.3 0.8 19.0 7.6 4.6 2.9 qu (MPa, %) f (deg, %) c (MPa, %) st (MPa, %) 186.8 184.4 199.1 194.8 34.9 30.3 29.5 33.0 48.5 53.4 58.5 53.3 65.5 54.3 44.7 35.4 12.7 7.9 6.5 4.1 18.1 23.8 16.3 17.3 12.8 18.9 14.5 14.4 26.1 12.2 7.4 7.6 Table 4 Effect of particle size on PFC3D macroproperties Particle size Davg (mm) 5.95 3.05 2.04 1.53 Macroproperty (63.4 31.7 E (GPa, %) n (, %) 57.3 64.0 67.6 69.2 0.231 0.254 0.255 0.256 10.0 3.9 1.8 1.2 31.7 mm3 specimens, n=10, mean and coefficient of variation) 21.2 8.1 5.8 5.5 qu (MPa, %) f (deg, %) c (MPa, %) st (MPa, %) 127.9 169.6 186.9 198.8 25.9 30.6 32.3 32.1 40.0 48.4 51.6 55.1 43.6 35.4 33.0 27.1 11.9 3.4 1.5 3.6 14.2 9.8 9.2 7.4 11.4 7.6 6.9 7.6 27.8 26.1 21.9 13.7 ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 clear increasing or decreasing trends. For the three smallest particle sizes, f and c exhibit the most variability of all measured properties. The Brazilian strength exhibits a clear dependence upon particle size, with st decreasing from 65.5 to 35.4 MPa as particle size decreases from 2.87 to 0.36 mm. The coefficients of variation have converged to approximately 7.5%. The variability of st is slightly greater than that of qu : This suggests that the criticalstate conditions in a Brazilian test are more sensitive to packing and strength heterogeneities than are the critical-state conditions in an unconfined compression test. The measured decrease in Brazilian strength with decreasing particle size is explained in Section 5.2, where it is shown that Brazilian strength is proportional to fracture toughness, which, in turn, is proportional to particle size. The lack of a similar size effect for qu suggests that the critical state in an unconfined compression test on the PFC2D material does not consist of the unstable growth of a single macrofracture subjected to extensile conditions. We speculate that a more complex critical state exists in which a failure plane and secondary macrofractures form under mixed compressive-shear conditions, and that such conditions are not sensitive to particle size. 3.5.2. PFC3D material response The Poisson’s ratio is independent of particle size. The Young’s modulus exhibits a clear dependence on particle size, with E increasing from 57.3 to 69.2 GPa as particle size decreases from 5.95 to 1.53 mm. The characterization of the elastic grain-cement micropron s perties in terms of E c ; Ē c ; ðkn =ks Þ; and ðk̄ =k̄ Þ provides a size-independent means of specifying the Poisson’s ratio of the material, but the Young’s modulus exhibits a minor size effect. The unconfined compressive strength exhibits a clear dependence on particle size, with qu increasing from 127.9 to 198.8 MPa as particle size decreases from 5.95 to 1.53 mm. The coefficients of variation have converged to approximately 3.5%. The reason for this size effect is unknown. It corresponds with general observations that finer-grained rock is stronger than coarser-grained rock (medium and fine-grained Lac du Bonnet granite [61,62]). A credible explanation must encompass the fact that qu of the PFC2D material is independent of particle size. No definitive statements about the effect of particle size on friction angle and cohesion can be made based on the results in Table 4, although both properties seem to increase as particle size decreases. The Brazilian strength exhibits a clear dependence upon particle size, with st decreasing from 43.6 to 27.1 MPa as particle size decreases from 5.95 to 1.53 mm. The coefficients of variation have not yet 1345 converged. The variability of st is much greater than that of qu : This suggests that the critical-state conditions in a Brazilian test are more sensitive to packing and strength heterogeneities than are the critical-state conditions in an unconfined compression test. 3.6. Effect of stress and damage on elastic constants The elastic constants of the BPM are affected by both stress and damage. These effects have been quantified using a ‘‘strain probe’’, which applies a specified strain path to the particles at the boundary of an extracted core and monitors the stress–strain response using measurement regions within the core. Preliminary measurements performed upon the granite model indicate a reasonable correspondence with the properties of Lac du Bonnet granite [53]. The elastic constants are affected by the stress state as follows. Compressive loading induces microtensile forces that act orthogonal to the compressive direction, and the particle motions that produce these microtensile forces modify the fabric (distribution of force chains) of the assembly. Increasing the mean stress increases the coordination number (average number of contacts per particle), which modifies the magnitudes of the elastic constants but does not produce anisotropy. Anisotropy is produced by deviatoric stresses that modify the directional distribution of the force chains—compressive loading increases the number of contacts oriented in the compressive direction and decreases the number of contacts oriented orthogonal to the compressive direction. These effects are enhanced by the presence of damage in the form of bond breakages. The elastic modulus dependence on mean and differential stress, as identified in PFC modeling, has been qualitatively supported by in situ measurements of dynamic modulus [63]. Ongoing work [64] aims to develop a quantitative link between AE/MS field measurements of dynamic moduli and rock-mass damage by comparing the evolution of stress- and damage-induced modulus anisotropy in the PFC3D model with the results of both short- and long-term laboratory tests. 4. Measured macroscopic properties (field scale) Three sets of boundary-value simulations were performed to demonstrate the ability of the twodimensional BPM to predict the extent and fabric of damage that forms adjacent to an excavation and the effect of this damage on rock strength and deformability [53]. The boundary-value simulations were of three experiments at Atomic Energy of Canada Limited’s (AECL) Underground Research Laboratory (URL) near Lac du Bonnet, Manitoba, Canada— namely, the Mine-by Experiment [65,66], the Tunnel ARTICLE IN PRESS 1346 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Sealing Experiment [67] and Stage 1 of the Heated Failure Test [68]. The Excavation Stability Study [69] also was simulated using slightly different microproperties. A summary of comparisons between simulations and in situ observations is provided in Ref. [63]. Extensive work has been done at the URL to characterize excavation-induced damage [70]. One common feature of this damage is the formation of breakout notches (apparently in regions where the elastic tangential stress at the excavation boundary exceeds some critical value). This feature is investigated in the simulations, which embed the BPM within a larger continuum model such that only the region of a single notch is represented by the BPM. Only the results of the Mine-by Experiment simulations are presented here. The Mine-by Experiment took place on the 420-m level of the URL in sparsely fractured Lac du Bonnet granite and consisted of carefully excavating a 3.5-m diameter circular tunnel subparallel with the direction of the intermediate principal stress. The tunnel was excavated in 1-m stages by drilling overlapping holes about the tunnel perimeter, and then removing the material using hydraulic rock splitters in closely spaced drill holes in the tunnel face. During tunnel excavation, a multistage process of progressive brittle failure resulted in the development of v-shaped notches, typical of borehole breakouts, in the regions of compressive stress concentration in the roof and floor (as shown in Fig. 12). The major and minor in situ stresses are approximately 60 and 11 MPa, which produce an elastic tangential stress of 169 MPa in the roof and floor. This induced stress is lower than the 200 MPa unconfined compressive strength of the undamaged granite, and observation and analysis of other excavations at the 420-m level indicate that notches will form when the elastic tangential stress in horizontal tunnels exceeds approximately 120 MPa [69]. Various mechanisms have been suggested to account for the apparent degradation in material strength sufficient to cause breakout. Two such mechanisms include precracking, which results from the local increase and rotation of stresses ahead of the advancing tunnel face, and stress corrosion—time-dependent cracking in rock that depends on load, temperature and moisture [71,72]. Both of these mechanisms are investigated in Ref. [73] using a simplified form of the BPM that represents the cement using contact bonds. (A contact bond provides tensile and shear strengths only and behaves like a parallel bond with n s k̄ ¼ k̄ ¼ R̄ ¼ 0:) In that work, a time-dependent strength-degradation process was included in the simplified BPM by reducing the contact-bond strength at a specified rate if the contact force was above a specified microactivation force. It was found that (a) precracking coupled with uniform strength reduction (choosing microproperties to match the long-term instead of the short-term strength of the rock) produced spalling and not distinct notches, and (b) activation of stress corrosion (for microproperties that match the shortterm strength of the rock) produced notches that stabilized after approximately 56 h. Subsequent work, summarized in Ref. [53], utilized the full BPM that represents the cement using parallel bonds, which allow one to mimic more closely the micromechanisms operative during stress corrosion by reducing the parallel-bond radius at a specified rate if the maximum tensile stress in the parallel bond is above a specified microactivation stress. In such models, notches are formed either by uniformly reducing the material strength everywhere or by activating stress corrosion. The difference in model behavior, particularly the ability to replicate notch development using a uniform strength reduction for the parallel bonds, is believed to be related to the use of parallel versus contact bonds. The parallelbonded material becomes much more compliant in the damaged region than does the contact-bonded material, thereby shedding more load to the notch perimeter. 4.1. Coupling the BPM with a continuum model Fig. 12. Photograph of broken and buckled rock slabs comprising the notch tip in the floor of the Mine-by Experiment tunnel. The Mine-by Experiment was simulated using a coupled PFC2D-FLAC [74] model in which only the potential damage region is represented by the PFC2D material. The rectangular sample produced by the material-genesis procedure is installed into a corresponding rectangular void created in a FLAC grid adjacent to the tunnel surface, and a circular arc is removed from the particle assembly, corresponding to the tunnel boundary (see Fig. 13). The basis of the coupling scheme is that FLAC controls the velocities of a layer of particles on three sides of the PFC2D region. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Fig. 13. PFC2D portion of coupled PFC2D-FLAC Mine-by model (PFC2D granite model with Davg ¼ 32 mm). PFC2D returns the unbalanced forces of the same set of particles, and these are applied as boundary forces to the FLAC grid. A symmetry line, in the direction of the major principal stress, is used based on the assumption that cracking is similar in the roof and floor of the tunnel. The FLAC grid boundaries are placed at a distance of 10 times the tunnel radius from the center. The elastic constants to be used in the FLAC grid are specified, and the FLAC model is run in large-strain mode under conditions of plane strain. Coupling occurs during cycling such that each cycle in PFC2D corresponds with a cycle in FLAC. Initially, the system is regarded as stress free. There are no initial stresses in the FLAC model, and the initial stresses in the PFC2D material are low relative to the final stresses that will develop as a result of the excavation. In situ stresses are applied to the outer boundaries of the FLAC grid, and cycling is performed until equilibrium is established, indicating that the excavation-induced stress redistribution is complete. An alternative procedure would be to initialize stresses in the FLAC grid and install the same mean stresses in the PFC2D sample (using the stress-installation procedure described in Appendix A) prior to its connection to the FLAC grid. The resulting relaxation would correspond to the tunnel being created in a stressed material. Three PFC2D granite materials were produced using the microproperties in Table 1 with different values of Dmin : The PFC2D materials have average particle diameters of 32, 16 and 8 mm, respectively, and this corresponds with model resolutions of 33, 66 and 132 particles across the 1.05-m extent of the potential damage region as shown in Fig. 13. The elastic constants assigned to the FLAC grid were taken as the average values of Young’s modulus and Poisson’s ratio from Table 2. Although the PFC2D particle size that generated the data in Table 2 was much smaller ðDavg ¼ 0:72 mmÞ; the results in Table 3 indicate that particle size 1347 has only a minimal effect on the macroscopic elastic constants. The modeling sequence mimics the rock genesis followed by the tunnel excavation. Two forms of strength reduction: uniform strength reduction, or activation of stress corrosion, were then applied to replicate the lower strength long-term response of the granite. The uniform strength reduction occurs by multiplying both the tensile and shear strengths of all parallel bonds by a specified factor. The activation of stress corrosion introduces a time scale into the simulations such that the time evolution of the damage-formation process can be investigated—and, in cases where all microtensile stresses fall below the microactivation stress, a time to full damage stability can be determined. In this paper, the results from fieldscale simulations using a uniform strength-reduction factor are presented, whereas the stress-corrosion model is described elsewhere [53] and is only briefly discussed in the following sections. 4.2. PFC2D model behavior during excavation-damage studies In BPM simulations of the Mine-by Experiment, the notch-formation process is driven by either (a) monotonically increasing the far-field loading, (b) uniformly reducing the microstrengths (as shown in Fig. 14, which includes the outline of the excavation, the outline of the PFC2D material, and a set of dashed circles spaced at intervals of R=5; where R is excavation radius), or (c) activating stress corrosion. The notch forms by a progressive failure process that starts at the excavation boundary in the region of maximum compression and proceeds inward. First, a small notch (group of Fig. 14. Effect of strength-reduction factor on damage patterns in the potential damage region of the Mine-by model for PFC2D granite model with Davg ¼ 8 mm: (a) Strength-reduction factor of 0.75 (177 cracks, tensile/shear=129/48); (b) Strength-reduction factor of 0.60 (3428 cracks, tensile/shear=2865/563); (c) Strength-reduction factor of 0.50 (8076 cracks, tensile/shear=6704/1372). ARTICLE IN PRESS 1348 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 microcracks) forms. Wing-cracks form near the notch tip and extend parallel with the current notch boundary, eventually curving toward and intersecting the boundary, forming slab-like pieces of material. These slabs do not detach fully, still being joined to the rock mass by some unbroken bonds, but do become much softer than the rock mass and, thus, shed load deeper into the rock mass. This increased load initiates additional wing cracks, which combine to extend the notch boundary. During this process, the notch is in a meta-stable state and only grows if the far-field load is increased, the microstrengths are reduced, or stress corrosion is active. The notch is meta-stable because the notch shape induces a compressive zone to form at the notch tip, which effectively strengthens the rock mass by reducing the microtensions that are driving the wing-crack growth. If the slabs are removed, the notch grows again to re-establish a new stable shape. This meta-stable growth process continues until the notch reaches a critical depth relative to the tunnel radius, at which time the wing-cracks do not curve back toward the notch boundary; instead, they curve away from the boundary, forming what looks like a macroscopic fracture (see Fig. 14(c)). This secondary fracture is described in Ref. [39] as a ‘‘shear band’’—an elongated cluster of microcracks that emanates from the apex of one or both notches and extends approximately parallel with the excavation surface. Although there is no clear evidence for, or against, the formation of a shear band in the Mine-by Experiment, such features are often observed in laboratory tests. The rupture zone that formed in both the numerical (using a simplified BPM that represents the cement using contact bonds instead of parallel bonds) and laboratory test of a rectangular prism of Berea sandstone containing a circular hole and loaded in a plane strain (biaxial) test [75] may be such a feature. Potyondy and Autio [39] offer the following hypothesis regarding this feature: It evolves from a band or ‘‘cloud’’ of microcracks. The formation of this band seems to be related to a punching-type failure occurring at the top and/or bottom of the [excavation], in which the applied load tries to drive an intact rectangular region of material into the unloaded region surrounding the [excavation]. (This unloaded region encompasses the notch tips and results from the presence of the notches, which are much more compliant than the surrounding rock and therefore shed load deeper into the rock away from the [excavation].) The microfailure mode within the shear band need not be of a shearing type; many tensile microcracks may combine in an en-echelon fashion to accommodate the macroscopic shearing motion and thereby produce the shear band. 4.3. Discussion of simulation results for the Mine-by models No significant damage forms as a result of excavation in any of the BPM simulations of the Mine-by Experiment; however, application of a strength-reduction factor of 0.6 to the material with Davg ¼ 8 mm produces a stable breakout notch as seen in Fig. 14(b). Activation of stress corrosion to the material with Davg ¼ 8 mm after excavation produces a notch that grows over time and does not fully stabilize before reaching the PFC2D-FLAC boundary, which occurs at a time of approximately 14 years. A notch of similar extent (but different internal damage fabric) as that produced by a uniform strength-reduction factor of 0.6 has formed after approximately 2 months of stress corrosion. The notch-formation process is sensitive to the strength-reduction factor, as seen in Fig. 14. It is not clear what procedure should be used to calibrate the strength-reduction factor, other than comparing the final extent of any notches that form in boundary-value models of excavations. The use of a uniform strengthreduction procedure to predict excavation damage may be limited to qualitative assessments of the nature of the possible damage, and use of a properly calibrated stresscorrosion model may be required to make quantitative assessments. The stress-corrosion model was calibrated to match the static-fatigue behavior of the granite [76]. The fact that similar notches form for both strengthreduction procedures supports such use of a uniform strength-reduction approach. The notch-formation process is sensitive to the particle size, as seen in Fig. 15. Notches form more readily in material with smaller particles. One hypothesis to explain this behavior is that the material with Fig. 15. Effect of particle size on damage patterns in the potential damage region of the Mine-by model for PFC2D granite models with strength-reduction factors of 0.6. (a) Davg=32 mm (75 cracks, tensile/ shear=59/16); (b) Davg=16 mm (412 cracks, tensile/shear=325/87); (c) Davg=8 mm (3428 cracks, tensile/shear=2865/563). ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 smaller particles has a lower fracture toughness (as discussed below). The extent of the notch is controlled by macroscopic fracture formation in the form of the wing-cracks described above. Stress concentrations occur at the tips of these wing-cracks, giving rise to the conditions for which the principles of fracture mechanics apply, such that the material resistance to such macroscopic fracture formation is measured by the fracture toughness rather than the unconfined compressive strength. 5. Emergent properties of the BPM Systems composed of many simple objects commonly exhibit behavior that is much more complicated than that of the constituents [77]. The particles and contacts comprising the BPM are described by simple equations, with few parameters, but an assembly of such particles displays a rich spectrum of behaviors that closely resemble those of rock. The following behaviors are observed both in rock samples and in synthetic samples composed of bonded particles: 1. Continuously non-linear stress–strain response, with ultimate yield, followed by softening or hardening. 2. Behavior that changes in character, according to stress state; for example, crack patterns are quite different in the tensile, unconfined- and confinedcompressive regimes. 3. Memory of previous stress or strain excursions, in both magnitude and direction. This behavior is commonly expressed in terms of moving yield surfaces, or evolving anisotropic damage tensors. 4. Dilatancy that depends on history, mean stress and initial state. 5. Hysteresis at all levels of cyclic loading/unloading; cyclic energy dissipation is strongly dependent on cyclic amplitude. 6. Transition from brittle to ductile shear response as the mean stress is increased. 7. Dependence of incremental stiffness on mean stress and history. 8. Induced anisotropy of stiffness and strength with stress and strain path. 9. Non-linear envelope of strength. 10. Spontaneous appearance of microcracks and localized macrofractures. 11. Spontaneous emission of acoustic energy. It may be noted that there is no existing continuum constitutive model that reproduces all of these behaviors. Some models can reproduce selected items, but, usually, many ad hoc parameters are needed. An assembly of bonded particles exhibits all of the 1349 behaviors, using few microparameters (but at the expense of often quite lengthy calibration procedures). Further, continuum-based models (typically implemented via finite or boundary element techniques) have difficulty reproducing the spontaneous development of many microcracks and macrofractures. BPMs explicitly represent the evolution of damage (in terms of the density and orientation of microcracks), in contrast to continuum models, which commonly contain evolution laws that are phenomenological, and not mechanistically based. Many of the noted behaviors arise from the existence of many sites of non-linear response, each of which may be activated at a different stress level. Thus, the memory of a previous stress level, at which unloading occurred, is encoded by a set of contacts that are either on the point of sliding or opening/closing. In either case, a change in response occurs at a particular stress level when a contact starts to slide or it opens or closes. Each such contact site contributes to the stress-dependent, nonlinear response of the entire assembly—not only in magnitude but also in direction, because the contact activation depends critically on its orientation in relation to the direction of applied strain. Similar arguments may be made for the micromechanical bases of the other phenomena noted above. The following focus is on a single aspect of the emergent behavior of the BPM—namely, the ability to reproduce the fracturing behavior of a brittle material and its application to explain the size effect observed in Brazilian tests. 5.1. The reproduction of fracture mechanics behavior BPMs produce patterns of bond-breaks that are qualitatively similar to microcracks and extended cracks seen in real rock samples. Although bond-breaks in a particle assembly might appear conceptually different from monolithic cracks that propagate in a brittle continuum, it is possible to establish a formal equivalence between the mechanisms and parameters of the BPM and the concepts and equations of LEFM. Consider the bonded-disk assembly shown in Fig. 16, in which a ‘‘crack’’ (an unbonded line of contacts) is introduced into a cubic packing of unit-thickness (t ¼ 1) disks of radius R joined by contact bonds with kn =ks ¼ 2:5 and subjected to an extension strain normal to the crack. Large tensile contact forces occur at contacts in front of the crack tip, as shown by the red lines. The force magnitudes conform well with those that would be induced over finite line segments in front of a crack tip in an isotropic linear elastic continuum; see Fig. 17, in which the LEFM values are derived as follows. For a through-thickness crack of half-width a in an infinitely wide plate of isotropic linear elastic material subjected to a remote tensile stress, sf ; acting normal to ARTICLE IN PRESS 1350 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 versus disk sequence number. The BPM values are from a symmetric model containing 64.5 disks along the crack width and model boundaries at a distance of 3a and 2a above and in front of the crack, respectively. Having established that forces in the BPM conform well to those induced on line segments in a continuum, the next step is to examine the maximum contact force existing at m ¼ 1; pffiffiffiffiffiffi ¼ 2sf t aR: (29) F max n Fig. 16. Force distribution and displacement field near the crack tip in a cracked cubically packed contact-bonded disk assembly. Noting that the mode-I stress pffiffiffiffiffiffi intensity factor for the system is defined as K I ¼ sf pa; Eq. (29) can be written as rffiffiffiffi R max : (30) F n ¼ 2tK I p The true tensile strength of the BPM (i.e., the strength from a test in which no force concentrations exist) is denoted by s0t : For a cubically packed BPM loaded parallel with the packing direction, s0t ¼ fn =2Rt; where fn is the contact-bond tensile strength (in force units). At the condition of incipient failure, or crack extension, F max ¼ fn and K I ¼ K Ic ; where K Ic is the mode-I n fracture toughness. Substituting into Eq. (30), pffiffiffiffiffiffiffi (31) K Ic ¼ s0t pR: Fig. 17. Normalized tensile force of LEFM and BPM at the five contacts ahead of the crack. the crack, the induced stress, sn ; acting on the crack plane near the crack tip ðr aÞ is [78] rffiffiffiffiffi a sn ¼ sf ; (26) 2r where r is the distance from the tip. The force, F n ; acting over a line segment equal to a disk diameter (2R) at a mean distance r ¼ ð2m  1ÞR is given by rffiffiffi Z rþR Z rþR a Fn ¼ sn t dr ¼ sf t r1=2 dr 2 rR pffiffiffiffiffiffi pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffirR ¼ 2sf t aRð m  m  1Þ; ð27Þ where m is a positive integer denoting the disk sequence number. The tensile forces in the BPM are compared with those of LEFM in Fig. 17 by plotting normalized force F~ n ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi Fn pffiffiffiffiffiffi ¼ m  m  1 2sf t aR (28) The fracture toughness is thus related to the properties of the BPM. Notably, because particle radius enters into the expression, particle sizes cannot be chosen arbitrarily if it is required to match fracture toughness. This is not surprising, as the concept of fracture toughness implies an internal length scale, whereby the ratio of fracture toughness to material strength has the dimension of square root of length. The particle size in the BPM supplies this internal length scale. The foregoing analysis is based on incipient failure and does not address the condition for propagation of a crack. Considering Eq. (29), the new induced force, F 0n ; following breakage of the bond nearest the crack tip, is found by substituting the new crack length for the original crack length: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F 0n ¼ 2sf t ða þ 2RÞR ¼ F max 1 þ 2R=a: (32) n Thus, the new maximum contact force is always greater than the previous maximum and greater than the bond strength. This condition ensures that the crack will propagate in an unstable fashion. The analysis assumes that the crack is remote from other boundaries and that the far-field load remains constant, but proximate boundaries or fixed-grip loading might cause F 0n to be less than F max if the breaking of the critical bond leads n to a force on the next bond that is less than the bond strength. In such a case, the crack would not propagate; thus, the condition for propagation is not necessarily a local one, but depends on the influence of the complete model on the contact forces near the crack tip. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 The derivation for fracture toughness in Eq. (31) is based on the assumption of a cubically packed contactbonded assembly with perfectly brittle bonds. The toughness for more general cases is obtained by testing a set of self-similar Center Cracked Tension specimens, as shown in the inset of Fig. 18. This system is characterized by two length scales f¼ a 1 ¼ ; W 3 c¼ a ; 2R (33) where W is the specimen half-width and c is the crack resolution. A set of self-similar specimens is generated by varying c while keeping f and particle size, R; fixed. For LEFM conditions, the applied far-field tensile stress at incipient failure, s0 ; is [78] K Ic pffiffiffi a1=2 ; CðfÞ p   1=2 pf ½1  0:025f2 þ 0:06f4 : CðfÞ ¼ sec 2 s0 ¼ ð34Þ This equation is normalized by dividing by the true tensile strength and setting a ¼ c2R to obtain " # s0 K Ic pffiffiffiffiffiffiffiffiffi c1=2 : (35) ¼ 0 s0t st CðfÞ 2pR For each specimen, s0 is measured and normalized strength versus crack resolution is plotted in log–log space. If the slope equals minus one-half, then LEFM conditions apply, and pffiffiffi pffiffiffiffiffiffiffi K Ic ¼ ½10b CðfÞ 2s0t pR; (36) where b is the y-intercept. The results of a set of self-similar tests on the cubically packed contact-bonded material with brittle bonds and bonds that have softening slopes one-fifth and one-tenth 1351 of the loading slope are shown in Fig. 18. Results similar to those for the softening slope of k=10 were obtained for an arbitrarily packed contact-bonded assembly with brittle bonds assigned both uniform and normally distributed strengths. These results support the postulate that pffiffiffiffiffiffiffiffiffi ðcontact-bonded materialÞ; (37) K Ic ¼ s0t paR where aX1 is a non-dimensional factor that increases with packing irregularity, strength heterogeneity and bond ductility. The value of aR represents the apparent internal length scale of the BPM. The a values are obtained from the y-intercepts of the tests by comparing Eqs. (36) and (37) to obtain  2 K Ic pffiffiffiffiffiffiffi ¼ 2ð10b Cð1=3ÞÞ2 : a¼ (38) s0t pR For the cubically packed contact-bonded material with brittle bonds, softening slope k=5 and softening slope k=10; using the last five data points, the slopes are –0.47, and the y-intercepts are –0.1590, 0.0980 and 0.0034, giving a values of 1.11, 1.46 and 2.26, respectively. Therefore, the effect of decreasing the softening slope is to increase the apparent internal length scale of the material, making it possible to match a given K Ic for several different particle sizes. The asymptotic slope of all curves (at large crack resolutions) approaches minus one-half, but they are shifted to the right for increasing ductility. This is consistent with an increase in apparent length scale. The expression for fracture toughness in Eq. (37) only applies to a contact-bonded material. The results of a set of self-similar tests on the cubically packed paralleln s bonded material (with l̄ ¼ 1; k̄ =k̄ ¼ 2:5 and s0t ¼ l̄ s̄c ) are also shown in Fig. 18. The parallel-bonded material Fig. 18. Normalized strength of self-similar numerical specimens. ARTICLE IN PRESS 1352 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 is weaker than the contact-bonded material for the same true tensile strength, because a non-zero bending s moment, M̄ ; develops at the crack tip and increases the maximum tensile stress acting on the bond periphery (see Eq. (16)). This leads us to postulate that pffiffiffiffiffiffiffiffiffi K Ic ¼ bs0t paR ðparallel-bonded materialÞ; (39) where aX1 is a non-dimensional factor that increases with packing irregularity, strength heterogeneity and bond ductility, and bo1 is a non-dimensional factor that accounts for the weakening effect of the bending moment. For the results shown in Fig. 18, a ¼ 1; and b is obtained from the y-intercept of the tests (0.3623) by comparing Eqs. (36) and (39) to obtain pffiffiffi K Ic b ¼ 0 pffiffiffiffiffiffiffi ¼ 10b Cð1=3Þ 2 st pR pffiffiffi ¼ 100:3623 Cð1=3Þ 2 ¼ 0:66: ð40Þ Fig. 18 also provides evidence for a size effect whereby the behavior becomes plastic (strength independent of crack resolution) at small crack resolutions and brittle (conforming to LEFM with a slope of minus one-half) at large crack resolutions. This behavior, in terms of specimen size, is reported extensively for laboratory tests [79]. In the present analysis, increasing crack resolution corresponds with increasing specimen size. The BPM has been shown to be equivalent to a material described by LEFM in terms of observable behavior and mathematical description. However, the results cast fracture mechanics in a new light. The singularities discussed extensively in the literature of fracture mechanics simply do not arise in a discontinuous medium. Therefore, there is no need to propose devices, such as a process zone or tip plasticity, to eliminate the crack-tip singularity. The analysis of the initiation and stability of a crack in a BPM may be done in terms of forces, not global energy balance, as commonly done for fractures in a continuum [80,81]. Finally, the BPM shows naturally how size effect arises, in terms of the ratio between specimen size and particle size. A discontinuum interpretation of fracture mechanics has been given previously, notably by Eringen [82], who derives an equation almost identical to Eq. (31) and remarks: ‘‘No poorly known constant such as surface energy appears. . .’’ (in reference to the existence of surface energy in the Griffith criterion). An extensive comparison between BPMs and the results of fracture mechanics is given in Ref. [83]. single macrofracture that travels across the specimen parallel with the direction of loading. The macrofracture is comprised of tensile microcracks and is driven by extensile deformation across the crack path similar to an LEFM mode-I crack as shown in Fig. 19. The conditions at peak load can be idealized as shown in Fig. 20, where the wedge-shaped damage regions are replaced by edge cracks of length a: If the two cracks are not interacting ðD  2a4aÞ; then the mode-I stressintensity factor at each crack tip is pffiffiffi K I / s a; (41) where s is the tensile stress acting across the uncracked ligament. At peak load, K I ¼ K Ic ; s ¼ st and a ffi 0:2D (observed for all particle sizes of the granite models), Fig. 19. Damage in Brazilian specimen just past peak load (148 cracks, tensile/shear=125/23). 5.2. Relating Brazilian strength to fracture toughness During Brazilian tests on BPMs, a wedge-shaped region of cracks forms inside of the specimen below each platen. Then, at peak load, one of the wedges initiates a Fig. 20. Idealized conditions in Brazilian specimen at peak load showing two cracks subjected to tensile loading across the uncracked ligament. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 which can be substituted into Eq. (41) to give K Ic st / pffiffiffiffi ; D (42) where D is the diameter of the Brazilian disk. The form of Eq. (42) is the same as the empirical relation between fracture toughness and tensile strength of st ¼ 6:88K Ic suggested by Zhang [84]. Eq. (42) suggests that the Brazilian strength should increase as specimen size decreases, and this trend is exhibited by Lac du Bonnet granite, for which the Brazilian strength is about 10 MPa for specimens greater than 50-mm diameter and increases to about 15 MPa for specimens smaller than 50-mm diameter [55]. The analysis leading to Eq. (39) supports the postulate that pffiffiffiffi K Ic / s̄c R; (43) where s̄c ¼ t̄c is the mean parallel-bond strength with l̄ ¼ 1; and R is the mean particle radius. Combining Eqs. (42) and (43) yields rffiffiffiffi R st / s̄c (44) D which suggests that the Brazilian strength is affected by the ratio of particle size to Brazilian disk diameter. If model resolution, C; is defined as the average number of particles across the Brazilian disk ðC ¼ D=2RÞ; then models with constant C will have the same Brazilian strength. For example, if all other microproperties are kept fixed, then the same Brazilian strength will be measured by either doubling the size of the Brazilian disk for a fixed particle size, or halving the particle size for a fixed Brazilian disk size. The BPM for granite exhibits such dual-model similarity as shown in Fig. 21. Fig. 21. Dual-model similarity: If all microproperties are kept fixed, then the exact same six short-term macroproperties in Table 3 are measured for both models B and C. 1353 The test results in Table 3 indicate that the Brazilian strength decreases as particle size decreases. The same trend is suggested by Eq. (44). The test results are compared with Eq. (44) by plotting ðst =s̄c Þ versus ðR=DÞ on a log–log scale. The test data are well fit by a straight line with slope of approximately 0.3, whereas Eq. (44) suggests a slope of one-half. The slope will only be onehalf if LEFM conditions apply, but such conditions do not apply here because (1) the microcracks at peak load do not form sharp macrofractures, and (2) the crack resolution ðc ¼ a=2R ffi 0:2D=2RÞ is relatively small ranging from 2.2 to 17.6. Eq. (44) can be expressed as st / c1=2 (45) s̄c and the test data compared with Fig. 18 to see that the test data lies in the range of c for which the slope is less than minus one-half (indicating a size effect whereby the behavior is more plastic than brittle) for the two reasons cited above. 6. Conclusions The BPM approximates the mechanical behavior of rock by representing it as a cemented granular material. The model has been used to predict damage formation adjacent to excavations in Lac du Bonnet granite. The damage-producing mechanisms activated in rock loaded in compression are complex and not understood fully. There is a complex interplay occurring between the macroscale compression and the induced microscale tensions, which drives the damage processes. This interplay produces non-linear, anisotropic elastic behavior under low loads before any significant damage forms, and it subsequently controls the damage that results under higher loads. In addition, it contributes to the long-term strength degradation processes, such as stress corrosion, that are activated by the increased stresses adjacent to an excavation but which appear to be inactive (at least on engineering time scales) under in situ conditions. Particle size controls model resolution but is not a free parameter; instead, particle size is related directly to the material fracture toughness. When modeling damage processes for which macroscopic fractures form, the particle size and model properties should be chosen to match the material fracture toughness as well as the unconfined compressive strength. Such processes occur at the boundary of an excavation and contribute to the formation of breakout notches. This poses a severe limitation on the size of a region that can be represented with a BPM, because the present microproperty characterization is such that the particle size must be chosen to be of the same order as the grain size. ARTICLE IN PRESS 1354 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 The BPM for Lac du Bonnet granite, as presented here, has the following limitations. The material strength matches that of Lac du Bonnet granite only for stresses near the uniaxial state—i.e., the tensile strength is too high, and the slope of the strength envelope as a function of confining stress is too low. However, there is evidence to support the postulate that including complex-shaped, breakable grains of a size comparable to that of the rock will remove this limitation. Also, computations are time-consuming and are limited to small volumes of rock, unless the BPM is embedded within a larger continuum model. Despite the limitations noted above, BPMs have a number of advantages when compared with conventional continuum approaches. Damage and its evolution are represented explicitly in the BPM as broken bonds; no empirical relations are needed to define damage or to quantify its effect on material behavior. Localized microcracks form and coalesce into macroscopic fractures without the need for re-meshing or grid reformulation. Complex non-linear behaviors arise as emergent features, given simple behavior at the particle level; thus, there is no need to develop constitutive laws to represent these effects. BPMs provide a rational means of incorporating additional physical processes that occur at the grain scale, such as fluid flow [85], stress corrosion [53] and thermal loading [86]. In addition, secondary phenomena, such as AE, occur in the BPM without additional assumptions. In principle, the BPM is capable of representing all of the significant physical behavior mechanisms in rock. The BPM is complete in the sense of Linkov [87], who describes a hypothetical model that encompasses all aspects in a causal chain of events that occur in a geomechanical medium. The links in his chain of processes are elastic distortion, creep deformation, dynamic instability, bifurcation and localization and transient equilibrium. Linkov emphasizes that a model embracing all of these features would be able to produce synthetic seismograms during a simulation [88], for example, of underground construction. The BPM provides such a complete model and reproduces qualitatively all of the mechanical mechanisms and phenomena that occur in rock, although adjustments and modifications may be necessary to obtain quantitative matches in particular cases. Ongoing work [64] is addressing three-dimensional modeling of excavation damage. There are two thrusts to this work. The first addresses computational limitations inherent in performing such simulations. A scheme called AC/DC (adaptive continuum/discontinuum) is being developed that replaces elastic regions (remote from cracks) by a continuum formulation, represented by a linear matrix, thus achieving economies in computation time. The switch between continuum and discontinuum is done automatically, such that a BPM is in place in time for cracking to occur. The second thrust focuses on rock physics and aims to develop a quantitative link between passive monitoring (of AE/ MS events and velocity surveys) and rock mass damage by developing a BPM that can reproduce the evolution of the stress- and damage-induced modulus anisotropy in the rock. If it can be shown that such a BPM exhibits the same modulus anisotropy as the rock, then it is likely that the damage in such a BPM is similar to the damage in the rock. Also, the relations between BPM properties and fracture toughness are being explored to determine if one can develop a BPM material with a fracture toughness that is independent of particle size. Acknowledgements Much of the development of the bonded-particle model has been funded by Atomic Energy of Canada Limited and Ontario Power Generation during the years 1995–2001 as part of its Thermal-Mechanical Stability Study, with a goal of producing a mechanistically based numerical model that can predict excavation-induced rock-mass damage and long-term strength of Lac du Bonnet granite. We also thank Fabian Dedecker (Itasca Consultants S.A.) for performing the triaxial and Brazilian tests on the PFC3D material as part of a project co-funded by the European Commission and performed as part of the fifth EURATOM framework program, Nuclear Fission (1998–2002). Appendix A All vector and tensor quantities are expressed using indicial notation with respect to a fixed right-handed rectangular Cartesian coordinate system: the position vector is denoted by xi and the stress tensor is denoted by sij : The Einstein summation convention is employed; thus, the repetition of an index in a term denotes a summation with respect to that index over its range. For example, the traction vector ti ; acting in the direction ni ; is given by ti ¼ sij nj ¼ si1 n1 þ si2 n2 þ si3 n3 : Vertical braces denote the magnitude of a vector or the absolute value of a scalar. The notation ½;i denotes differentiation with respect to the coordinate xi and a dot over a variable indicates a derivative with respect to time (e.g., x_ i ¼ @xi =@t). The Kronecker delta, dij ; and the permutation symbol, eijk ; are employed. A.1. Stress-measurement procedure The average stress, s̄ij ; in a volume, V ; of material is defined by Z 1 sij dV (46) s̄ij ¼ V V ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 where sij is the stress acting throughout the volume. For a particulate material, stresses exist only in the particles; thus, the integral can be replaced by a sum over the N p particles contained within V as 1 X ðpÞ ðpÞ s̄ij V ; s̄ij ¼ V N (47) p where s̄ðpÞ ij is the average stress in particle ðpÞ: In the same way, the average stress in a particle ðpÞ can be written using Eq. (46) as Z 1 ðpÞ ðpÞ s̄ij ¼ ðpÞ sðpÞ (48) ij dV : ðpÞ V V The identity S ij ¼ dik S kj ¼ xi;k Skj ¼ ðxi Skj Þ;k  xi Skj;k (49) holds for any tensor S ij : Applying this identity to the stress in a particle, one can write  Z 
 1 ðpÞ ðpÞ ðpÞ s̄ij ¼ ðpÞ (50) xi skj  xi skj;k dV ðpÞ : ;k V V ðpÞ The stress in each particle is assumed to be continuous and in equilibrium. In the absence of body forces, the equilibrium condition is sij;i ¼ 0: The volume integral in Eq. (50) is re-written as a surface integral by applying the Gauss divergence theorem to the first term and noting that the second term vanishes in the absence of body forces such that Z
 1 ðpÞ xi sðpÞ nk dSðpÞ s̄ij ¼ ðpÞ kj V S ðpÞ Z 1 ðpÞ xi tðpÞ ð51Þ ¼ ðpÞ j dS ; ðpÞ V S where S ðpÞ is the particle surface, nk is the unit outward normal to the surface and tjðpÞ is the traction vector. If the moment carried by each parallel bond is neglected, then each particle is loaded by point forces acting at discrete contact locations, and the above integral can be replaced by a sum over the N c contacts as s̄ðpÞ ij 1 X ðcÞ ðcÞ ¼  ðpÞ xi F j ; V Nc (52) ðcÞ where xðcÞ i is the location and F j is the force acting at ðcÞ contact ðcÞ: F j includes both F i and F̄ i from Eqs. (4) and (13), and the negative sign is introduced to ensure that compressive/tensile forces produce negative/positive stresses. The contact location can be expressed as    ðcÞ ðpÞ ðpÞ  ðc;pÞ xðcÞ ; (53) i ¼ xi þ xi  xi ni is the location of the particle centroid and where xðpÞ i nðc;pÞ is the unit-normal vector directed from the particle i 1355 centroid to the contact location. By substituting Eqs. (53) into (52) and noting that X ðcÞ Fj ¼ 0 (54) Nc for a particle in equilibrium, one obtains  1 X ðcÞ ðpÞ  ðc;pÞ ðcÞ s̄ðpÞ xi  xi ni F j : ij ¼  ðpÞ V Nc (55) An expression for the average stress in a volume, V ; is obtained by substituting Eq. (55) into (47); however, definition of the volume in the resulting expression is problematic because of the particles that intersect the measurement region. The problem is overcome by noting that, in a statistically uniform assembly, the volume associated with each particle is V ðpÞ =ð1  nÞ; such that V ðpÞ ; (56) 1n where n is the porosity of the region. The average stress in a measurement region is found by combining Eqs. (47), (55) and (56) to yield 0 1   B 1  n C X X ðcÞ ðpÞ  ðc;pÞ ðcÞ (57) s̄ij ¼ @P ðpÞ A xi  xi ni F j ; V Np Nc V¼ Np where the summations are taken over the N p particles with centroids contained within the measurement region and the N c contacts of these particles, n is the porosity within the measurement region, V ðpÞ is the volume of particle ðpÞ; xðpÞ and xðcÞ i i are the locations of a particle centroid and its contact, respectively, nðc;pÞ is the uniti normal vector directed from a particle centroid to its contact location, and F ðcÞ j is the force acting at contact ðcÞ; which includes both F i and F̄ i from Eqs. (4) and (13) but neglects the parallel-bond moment. A.2. Strain-rate measurement procedure The procedure employed to measure local strain rate within a particle assembly differs from that used to measure local stress. In determining local stress, the average stress within a volume of material is expressed directly in terms of the discrete contact forces, as the forces in the voids are zero. It is not correct, however, to use the velocities in a similar way to express the average strain rate, as the velocities in the voids are non-zero. Instead of assuming a form for the velocity field in the voids, a velocity-gradient tensor, based on a best-fit procedure that minimizes the error between the predicted and measured velocities of all particles with centroids contained within the measurement region, is computed. The strain-rate tensor is the symmetric portion of the velocity-gradient tensor. ARTICLE IN PRESS 1356 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Before describing the least-squares procedure, the relation between the strain and strain rate will be reviewed, and the corresponding terms defined. The relation between the displacements ui at two neighboring points is given by the displacement gradient tensor, aij : Let points P and P0 be located instantaneously at xi and xi þ dxi ; respectively. The difference in displacement between these two points is dui ¼ ui;j dxj ¼ aij dxj : (58) The displacement-gradient tensor can be decomposed into a symmetric and an anti-symmetric tensor as aij ¼ eij  oij where eij ¼ 12ðui;j þ uj;i Þ; infinitesimal-strain tensor and oij ¼ 12 ðuj;i  ui;j Þ; rotation tensor: (59) In a similar fashion, the relation between velocities vi at two neighboring points is given by the velocitygradient tensor, a_ ij : (The velocity field is related to the displacement field by ui ¼ vi dt; in which vi is the velocity and dt is an infinitesimal interval of time.) Let points P and P0 be located instantaneously at xi and xi þ dxi ; respectively. The difference in velocity between these two points is dvi ¼ vi;j dxj ¼ a_ ij dxj : (60) The velocity-gradient tensor can be decomposed into a symmetric and an anti-symmetric tensor as _ ij a_ ij ¼ e_ij  o where e_ij ¼ 12 ðvi;j þ vj;i Þ; strain-rate tensor and   _ ij ¼ 12 vj;i  vi;j ; o spin tensor: (61) a_ ij is computed using the following least-squares procedure. The velocity-gradient tensor computed over a measurement region represents the best fit to the N p ðpÞ measured relative velocity values V~ i of the particles with centroids contained within the measurement region. The mean velocity and position of the N p particles is given by P ðpÞ P ðpÞ Vi xi V̄ i ¼ Np and Np x̄i ¼ Np Np ; (62) where V iðpÞ and xðpÞ are the translational velocity and i centroid location, respectively, of particle ðpÞ: The measured relative velocity values are given by ðpÞ V~ i ¼ V ðpÞ i  V̄ i : (63) For a given a_ ij ; the predicted relative velocity values v~ðpÞ i can be written, using Eq. (60), as _ ij x~ ðpÞ _ ij ðxðpÞ v~ðpÞ i ¼a j ¼a j  x̄j Þ: (64) A measure of the error in these predicted values is given by   
X
ðpÞ X  ðpÞ ðpÞ 2 ðpÞ ~ ðpÞ ;  V v~ðpÞ v~i  V~ i z¼ v~i  V~ i  ¼ i i Np Np (65) where z is the sum of the squares of the deviations between predicted and measured velocities. The condition for minimum z is that @z ¼0 : @_aij (66) Substituting Eq. (64) into Eq. (65) and differentiating, the following set of nine equations is obtained: 2 P ðpÞ ðpÞ P ðpÞ ðpÞ P ðpÞ ðpÞ 3 x~ 3 x~ 1 x~ 2 x~ 1 x~ 1 x~ 1 Np Np 78 _ 9 6 Np 7> ai1 > 6 = < > 6 P ðpÞ ðpÞ P ðpÞ ðpÞ P ðpÞ ðpÞ 7> 7 6 x~ 1 x~ 2 ~ ~ ~ ~ x x x x 2 2 3 2 a_ i2 7 6N Np Np > 7> 6 p ; : > 7> 6 4 P ðpÞ ðpÞ P ðpÞ ðpÞ P ðpÞ ðpÞ 5 a_ i3 x~ 1 x~ 3 x~ 2 x~ 3 x~ 3 x~ 3 Np Np Np 9 8 P ðpÞ > V~ i x~ 1ðpÞ > > > > > > > Np > > > > > > > = < P ðpÞ ðpÞ > ~ V i x~ 2 : ¼ Np > > > > > > > > > > > > > P V~ ðpÞ x~ ðpÞ > > ; : i 3 > ð67Þ Np These nine equations are solved by performing a single LU decomposition upon the 3 3 coefficient matrix and performing three back substitutions for the three different right-hand sides obtained by setting i ¼ 1; 2 and then 3. In this way, all nine components of the velocity-gradient tensor are obtained. (For the PFC2D model, four equations are used to obtain the four non-zero values of the velocity-gradient tensor.) A.3. Stress-installation procedure The stress-installation procedure installs a general state of uniform stress (i.e., stress that does not vary with position) within an arbitrarily shaped particle ensemble. The intended use is to initialize internal stresses—not to apply stress boundary conditions. The procedure is based upon the following relations. The displacementgradient tensor, aij ; can be decomposed into a symmetric and an anti-symmetric tensor [89] as @ui ¼ ui;j ¼ eij  oij @xj     ¼ 12 ui;j þ uj;i  12 uj;i  ui;j ; aij ¼ ð68Þ where eij is the infinitesimal-strain tensor and oij is the rotation tensor. The velocity-gradient tensor, a_ ij ; relates the velocities vi ¼ u_ i at two neighboring points. Let the ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 points P and P0 be located instantaneously at xi and xi þ dxi ; respectively. The difference in velocity between these two points is dvi ¼ @vi dxj ¼ vi;j dxj ¼ a_ ij dxj : @xj (69) The velocity at an arbitrary point xj (expressed as vi ðxj Þ) is found by integrating Eq. (69) with respect to position. Let x̄j be a fixed reference position, then Z xj Z xj Z xj @ _ ðeij  oij Þ dxj ; dvi ¼ aij dxj ¼ @t x̄j x̄j x̄j Z xj vi ðxj Þ  v̄i ¼ ð70Þ e_ij dxj ; x̄j where a condition of no rotation has been enforced by setting oij ¼ 0 in Eq. (68). If the time derivative of the strain tensor is approximated by Deij ; (71) NDt where Deij is a strain increment applied over N time steps, each of size Dt; then the velocity field can be written as   Deij (72) vi ðxj Þ ¼ v̄i þ ðxj  x̄j Þ NDt e_ij ffi by assuming that the strain field does not vary with position. (If the strain field did vary with position, one could express the variation explicitly before performing the integration and thus obtain a final expression for the velocity that would depend upon the parameters used to characterize the strain-field variation.) The stress-installation procedure utilizes an iterative approach whereby a set of boundary (and possibly also interior) particles is moved, then the boundary particles are fixed, the interior particles are freed and staticequilibrium conditions are allowed to develop. During each iteration, the applied particle displacements are computed from the strain increment that is related by linear elasticity to the stress increment needed to reach the target stress. The iterations continue until the target stress is obtained. The stress field is measured by averaging the stresses within a set of measurement regions that cover the ensemble. The stress increment is defined as the difference between the target stress and the current stress Dsij ¼ stij  sij (73) and the corresponding strain increment is found by assuming isotropic, elastic behavior:   1þn n (74) Dsij  Dsaa dij : Deij ¼ E E The applied particle velocities are given by Eq. (72), where x̄j and v̄i are chosen as the initial position and velocity of the ensemble centroid, respectively. Good 1357 performance is achieved by setting n ¼ 0; N ¼ 10 and E ¼ 2E 0 ; where E 0 is an estimate of the ensemble modulus. (Note that when E is greater than the true ensemble modulus, the strain required to reach the target stress will be underestimated; thus, the stresses will converge to the target stress from below in a smooth fashion. If E is too small, the stresses at the end of each iteration will overshoot the target stress; the stresses will then oscillate about the target stress and not converge.) A stress field can also be installed within an ensemble containing holes by including the hole boundaries in the set of boundary particles. For such a system, the installed stress field will not sense the presence of the holes. Such a scheme has been used to model a tunnel excavation without having to include particles within the tunnel itself. After installing the desired stress field, the tunnel boundary particles are freed to allow stress redistribution to occur as the tunnel is now sensed by the system. During the procedure, displacements can be applied to designated boundary particles or to all particles. For cases in which the material remains primarily elastic, applying displacements to all particles will: (1) speed convergence (as it will not be necessary to allow the entire system to respond to boundary displacements in order to reach final equilibrium; instead, the entire system is made to correspond with the required strain field instantaneously during each iteration); and (2) produce a stress field that is nearly uniform throughout the ensemble. Such would not be the case if either the boundary particles alone or a set of confining walls were moved, because then non-uniform force chains would develop throughout the ensemble, producing a locally non-uniform stress field—for compressive conditions, a compressive cage that shields the internal region from the full compressive stresses often forms. Note, however, that for cases in which localized failure occurs, the strain field is locally non-uniform. Thus, for such cases, it may be more appropriate to move only the boundary particles. A.4. Ball-generation procedure A procedure is described to generate a collection of particles of a given uniform size distribution with radii in the range ½Rmin ; Rmax  that fill a given area, A; or volume, V ; at a given porosity, n: The basic algorithm is presented first. It makes use of two relations that are derived at the end of this section. The number of particles, N; that satisfies the above constraints is given by 8 Að1  nÞ > > ; PFC2D < 2 Rmin þ Rmax pR̄ with R̄ ¼ : N¼ > 3V ð1  nÞ 2 > : ; PFC3D 3 4pR̄ (75) ARTICLE IN PRESS 1358 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 These particles are chosen from a uniform size distribution with radii in the range 12 ½Rmin ; Rmax  and placed randomly within the specified region such that they do not overlap one another or the region boundary. The particles are generated at half their target sizes to ensure that the no-overlap criterion can be satisfied. Next, the porosity of the generated assembly, n0 ; is computed, and the radii of all particles are multiplied by the factor   1  n 1=l m¼ where 1  n0 ( 2; PFC2D; ð76Þ l¼ 3; PFC3D: The radius multiplier must be computed (i.e., it will not equal 2) because Eq. (75) is an approximation. Eq. (76) is derived as follows. The porosity, n; of a collection of particles contained within an area, A; or volume, V ; is defined as 8 A  Ap > < ; PFC2D; A n¼ (77) > : V  V p ; PFC3D; V where Ap and V p are the total area and volume, respectively, of all particles given by X Ap ¼ pR2 ; X4 Vp ¼ pR3 ; ð78Þ 3 P where is over all particle radii, R: By combining Eqs. (77) and (78), one can write X Að1  nÞ ; R2 ¼ p X 3V ð1  nÞ R3 ¼ : ð79Þ 4p Denoting the old porosity and radii by n0 and R0 ; the new porosity and radii by n and R; and using Eq. (79), one can write P l R 1n (80) ¼P l: 1  n0 R0 If the same radius multiplier, m; is used for all particles, then R ¼ mR0 ; which can be substituted into Eq. (80) to obtain Eq. (76). Thus, given a collection of particles of porosity, n0 ; within an area, A; or volume, V ; Eq. (76) is an exact expression of the radius multiplier for all particles necessary to obtain a porosity,Pn: Eq. (75) arises from the approximation that Rl of a P l uniform distribution is equal to R̄ of a collection of same-size particles of size R̄: This approximation can be expressed as X X l l Rl ¼ (81) R̄ ¼ N R̄ ; which can be substituted into Eq. (79) to obtain Eq. (75). A.5. Isotropic stress installation procedure The isotropic stress, s0 ; of a particle assembly is defined as the average of the direct stresses s0 ¼ s̄kk l where l ¼ 2; PFC2D; 3; PFC3D: (82) In the scheme described here, the radii of all particles are scaled iteratively to modify the isotropic stress of the assembly. Combining Eqs. (55) and (47), an expression is obtained for the average stress, s̄ij ; in a volume, V ; of material 1 X X ~ ðc;pÞ ðc;pÞ ðcÞ R ni F j ; (83) s̄ij ¼  V N N p c   ðc;pÞ  ðpÞ  where R~ ¼ xðcÞ  x i i : The isotropic stress can then be evaluated by 1 X X ~ ðc;pÞ nðcÞ s̄kk s0 ¼ ¼ R F (84) lV N N l p c nðcÞ is the normal component of the force acting where F at contact ðcÞ: If the radii of all particles is scaled by the same factor, a; such that the change in radius is DRðpÞ ¼ aRðpÞ (85) and it is assumed that all particles remain fixed, then the change in the normal component of force acting at each contact will be DF nðcÞ ¼ aK nðcÞ fðcÞ ; ( ðAÞ R þ RðBÞ ; fðcÞ ¼ RðpÞ ; particle2particle contact; particle2wall contact; ð86Þ where K nðcÞ and fðcÞ are the normal stiffness and particle radii, respectively, at contact ðcÞ: Substituting in the incremental form of Eq. (84) gives a X X ~ ðc;pÞ nðcÞ ðcÞ Ds0 ¼  R K f : (87) lV N N p c Rearranging this expression gives a ¼ PP lV Ds0 : ðc;pÞ R~ K nðcÞ fðcÞ (88) Np Nc Ideally, this formula enables one to determine the change in radii that will produce a given change in isotropic stress. However, some particle rearrangement occurs when the radii are changed, so formula (88) is applied several times, until the measured isotropic stress is within some tolerance of the target isotropic stress. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 A.6. Floater-elimination procedure When a particle assembly is compacted, even under the condition of zero friction, typically 10% to 15% of the particles are found to have no contacts. These particles are termed ‘‘floaters’’, because they are detached from the material matrix and appear to be floating in space. In physical specimens of an unbonded, granular material, such as sand, floaters are believed to exist, and it is reasonable, therefore, to accept their presence in a numerical specimen. However, if a particle model is used to represent a solid material, such as rock, each floater is equivalent to a void in the material. The effect of such voids on material response can be minimized by adjusting microproperties to obtain correct ensemble properties—nevertheless, the voids introduce a pattern of inhomogeneity that has no physical counterpart. The effect of this inhomogeneity is unknown and may be small. To be safe, steps can be taken to eliminate its source; an algorithm to eliminate floaters is described here. It may be argued that a ‘‘solid’’ consisting of bonded circular or spherical particles, even with no floaters, contains many voids. However, in a dense assembly, these voids are inaccessible during deformation, because they are too small to allow particles to move into them. In this sense, they are invisible as far as potential mechanisms are concerned. The starting point for the algorithm is a stable, compacted assembly with no bonds in which all contacts carry similar levels of force (i.e., the force distribution is fairly uniform). Further, the average stress level is low compared to the final stress level to be carried by the material. The algorithm expands and moves floaters until every particle has the specified minimum number of contacts. The key to success (i.e., elimination of all floaters) is the recognition that it is not necessary to enforce local equilibrium. If the assembly were destined to become an unbonded granular medium, then equilibrium of the final state would be required (and floater elimination made more difficult). However, a bonded material always carries locked-in forces of magnitudes that are related to the levels of contact forces existing at the instant of bonding. These locked-in forces typically are made low compared to the forces carried by the material in its operating range. Extra forces of low magnitude (e.g., one-tenth of pre-bonding contact forces) introduced by the floater-elimination algorithm add little to the existing ‘‘noise’’ of locked-in forces, even though the forces are not in local equilibrium. After bonding, equilibrium is again established, with final forces at levels comparable to their prebonding levels. In summary, the algorithm is described as follows. All particles except floaters are fixed and given zero velocities. Floaters are then expanded by a large amount (30%), sufficient to ensure contact with all of their 1359 immediate neighbors. After being cycled to local equilibrium, the floaters are contracted by an amount (see Eq. (94)) that is calculated to reduce their mean contact normal force below a target force (one-tenth of the mean contact normal force for the assembly). If the mean contact normal force for a particle is below the target force, then it is declared ‘‘inactive’’ and is not contracted further at this stage. The contraction step is applied repeatedly until all floaters become inactive. After executing this iteration, the number of active floaters is reduced considerably. The complete procedure is repeated several times until all floaters have been removed. The algorithm works well, because the adjustments to floater radii (and positions, through the law of motion) are done for small groups within a ‘‘cage’’ of fixed particles. As mentioned, it is only necessary to establish equilibrium within this trapped group. Then, as floaters themselves become part of the cage, the scheme becomes increasingly more efficient as the number of particles in each floater group decreases. The following parameters control the algorithm: N f (minimum number of contacts to be a non-floater, set to 3 in our granite models); M r (initial radius multiplier, set to 1.3); f m (target fraction of F̄ a ; set to 0.1); F r (relaxation factor, set to 1.5); and H (hysteresis factor, set to 0.9). The following variables are used in the algorithm: F̄ a (mean contact normal force for entire assembly); F̄ p (mean contact normal force for single particle); R (particle radius); and kn (particle normal stiffness). The algorithm is defined by the following pseudo-code. Perform this outer loop up to 10 times. Get current mean contact normal force for the assembly, F̄ a : Scan all particles to identify floaters: a floater is a particle with less than N f contacts. If no floaters, then EXIT outer loop. Fix all non-floaters. For all floaters: –expand by M r , –set friction to zero. Cycle 200 steps to get floater groups to equilibrium. Perform this inner loop up to 100 times. Define ‘‘active’’ floaters. An active floater is a particle with more than one contact and mean contact normal force greater than f m F̄ a : If no active floaters, then EXIT inner loop. Shrink each active floater according to the formula: R :¼ R  F r ðF̄ p  H F̄ a f m Þ=kn : ARTICLE IN PRESS 1360 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 Cycle 100 steps to obtain approximate equilibrium. End of inner loop. End of outer loop. In order to speed up the algorithm, contact calculations are inhibited whenever the two particles comprising any contact are fixed. The inhibit flags are removed upon exit from the algorithm. For a large assembly, it is found that one or two floaters may be positioned almost exactly between two non-floaters such that the two contact directions are nearly coaxial. The algorithm cannot establish more than two contacts in this case. To allow convergence, the parameter N f is reduced by one if the same number of overall floaters is observed between iterations of the outer loop. The objective of the shrinkage procedure is to reduce the mean contact normal force on a particle below some target force, F target ¼ F̄ a f m ; by setting R R  F r ðF̄ p  H F̄ a f m Þ=kn : (89) This expression is now derived. The mean normal force on a single particle is given by 1 X nðcÞ F̄ p ¼ F (90) Nc N c nðcÞ is the normal component of the force acting where F at contact ðcÞ: The normal force at each contact is changed by the following amount, if the radius is changed by DR; assuming that particles remain fixed during the radius change: DF n ¼ 12 kn DR; (91) where kn is the particle normal stiffness. Thus, the new 0 mean normal force, F̄ p ; is  1 X  nðcÞ 0 (92) F̄ p ¼ F þ DF n : Nc N c Then, for the new mean normal force to be below the 0 target normal force, set F̄ p ¼ HF target ; where H is a factor less than unity. Thus,  kn DR 1 X  nðcÞ 1 X nðcÞ þ F þ DF n ¼ F HF target ¼ Nc N 2 Nc N c c (93) and DR ¼ 2ðF̄ p  H F̄ a f m Þ=kn : (94) This expression is identical to Eq. (89), except that F r replaces the factor of two. It is found that setting F r equal to 1.5 gives smoother convergence, because the radii of several neighboring particles are likely to be adjusted simultaneously. A value of F r less than 2 provides some damping in the algorithm, preventing excessive oscillation and possibly instability. A value of H ¼ 0:9 also prevents ‘‘noise’’, caused when particles satisfy the force criterion by only a small margin (thus possibly requiring multiple adjustments). A.7. Biaxial, triaxial and Brazilian testing procedures The macroscopic properties of the model material are obtained by performing biaxial, triaxial and Brazilian tests, in which each specimen is confined and loaded by pairs of opposing frictionless walls. These are the walls of the material vessel used during material genesis. For the biaxial and triaxial tests, the top and bottom walls act as loading platens, and the velocities of the lateral walls are controlled by a servomechanism that maintains a specified confining stress. Strictly speaking, the biaxial and triaxial tests simulate a polyaxial loading test—not a triaxial test, in which a specimen is encased in a membrane and confined by fluid pressure. Such a triaxial test could be simulated by replacing the side walls with a sheet of particles joined by parallel bonds to mimic the elastic membrane. The current biaxial and triaxial tests inhibit specimen bulging, whereby the specimen sides deform into a barrel shape, because the side walls remain straight. Such bulging should be minimal for granite, but may be important for softer rocks or cohesive soils. The lateral wall normal stiffnesses are set equal to a fraction (0.001–0.02 in PFC2D and 0.001–0.10 in PFC3D for confinements of 0.1–70 MPa) of the average particle normal stiffness to simulate a soft confinement. For the Brazilian test, the specimen is trimmed into a circular (in PFC2D) and cylindrical (in PFC3D) shape that is in contact with the lateral walls, and the lateral walls act as loading platens. For all tests, the loadingplaten normal stiffnesses are set equal to the average particle normal stiffness. The tests are run under displacement control such that the platens move toward one another at a constant rate. Hazzard et al. [49] ran similar tests with a servo-control to move the platens, so as to maintain a constant cracking rate, and observed a snap-back in the stress–strain response. Platen velocities of 0.05 m/s are chosen to ensure quasi-static test conditions, which are confirmed by demonstrating that reducing the platen velocities does not alter the measured macroproperties. The advantage of using the material-vessel walls to perform these tests is that the boundary particle alignment produces a rather uniform force transmission from the walls to the specimen—i.e., there are no particles protruding from the specimen and producing localized loading. However, the boundary particles tend to be aligned with these walls (see Fig. 5(d)) and, thus, are not fully representative of the internal microstructure. The model material mimics a sandstone created by first filling a vessel with sand, compacting the sand, and then cementing the sand at grain-grain contacts. A rock specimen that has been cored and polished is also not ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 fully representative of the internal microstructure, but for a different reason—namely, there is grain-scale damage at the specimen boundaries. Both of these specimen-creation processes will affect measured properties, but if the grains in the real rock and the particles in the modeled rock are small relative to the specimen dimensions, then these effects should be minimal. A more representative measure of the internal microstructural properties of the model material can be obtained by removing all particles outside of a nominal core region, identifying a set of boundary particles and specifying their motion while monitoring stress and strain within the core using a stress-installation procedure. Such measures are used to measure the evolving properties within selected regions of the model surrounding an excavation. Stresses and strains are computed using the specimen dimensions at the start of the test. Stresses are computed by dividing the average of the total force acting on opposing walls by the area of the corresponding specimen cross section. Stresses in the PFC2D models are computed assuming that each particle is a disk of unit thickness. Strains are computed by monitoring the motion of pairs of gauge particles that lie along the global coordinate axes at the centers of the specimen faces. (For both PFC2D and PFC3D models, the specimen axis is parallel with the global y-axis.) Computing strains via gauge particles removes the error introduced by particle–wall overlap when using the distances between opposing walls to compute strain. The elastic constants are computed using the stress and strain increments occurring between the start of the test and the point at which one-half of the peak stress has been obtained. For the PFC3D material, the Young’s modulus is computed by E¼ Dsy Dy ðPFC3D materialÞ (95) and the Poisson’s ratio is computed by 1 ðDx þ Dz Þ Dx ¼ 2 u¼  Dy Dy   1 Dv 1 ¼ ðPFC3D materialÞ; 2 Dy ð96Þ where the average of the two lateral strains is used to approximate the lateral strain, and volumetric strain Dv ¼ Dx þ Dy þ Dz : These relations are valid for constant lateral stress during the test. For the PFC2D material, first the Poisson’s ratio and Young’s modulus corresponding with a state of plane stress are computed as Dx Dy Ds y E0 ¼ : Dy n0 ¼  ðPFC2D material; plane stressÞ; ð97Þ 1361 Eq. (97) is valid for the case of plane stress ðsz ¼ 0Þ and constant lateral stress ðDsx ¼ 0Þ during the stress and strain increments. Then, the elastic constants corresponding with a state of plane strain are obtained via n0 ðPFC2D material; plane strainÞ; 1 þ n0 E ¼ E 0 ð1  n2 Þ: n¼ ð98Þ Eq. (98) is the general relation between plane-stress and plane-strain conditions [90]. When calibrating the PFC2D model, E and n (not E 0 and n0 ) are compared with the elastic constants measured during triaxial tests on real rock specimens. Recall that the PFC2D model behaves as an assembly of unit-thickness disks. The conditions are neither plane strain nor plane stress, because the corresponding stress–strain constitutive relations are not employed— i.e., the out-of-plane forces and displacements do not enter into the force–displacement law. Therefore, the elastic constants of the PFC2D model are measured by interpreting the force–displacement response. For the same PFC2D model, one force–displacement response is measured, but two sets of elastic constants, corresponding with plane-stress and plane-strain conditions, are computed. If a region of a 2D isotropic elastic continuum were extracted and replaced with this model material, and the corresponding set of elastic constants were assigned to the remaining continuum, then the deformation state of the model material would match that of the continuum—i.e., there would be full displacement compatibility along the extraction boundary. This procedure is used in the construction of boundary-value models of excavation damage in which the model material is inserted into a pre-defined region within a larger continuum model. For a tunnel simulation, the continuum model is run in planestrain mode and assigned the plane-strain elastic constants E and n: The initiation of cracking in the PFC models is controlled by the ratio of standard deviation to mean material strength. Increasing this ratio lowers the stress at which the first crack initiates. However, the axial stress at the time of the initiation of the first crack in a synthetic specimen will underestimate the crack-initiation stress, sci ; measured in typical laboratory experiments. Therefore, sci measured during a biaxial or triaxial test of a PFC model is defined as the axial stress at which there exists a specified fraction (chosen as 1% in the present simulations) of the total number of cracks existing at peak stress. The choice of 1% is rather arbitrary, and sci serves as a qualitative calibration property providing a means to match the onset of significant internal cracking before ultimate failure. (For Lac du Bonnet granite, sci ¼ 0:45qu :) ARTICLE IN PRESS 1362 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 References [1] Okui Y, Horii H. A micromechanics-based continuum theory for microcracking localization of rocks under compression. In: Mühlhaus HB, editor. Continuum models for materials with microstructure. New York: Wiley; 1995. p. 27–68. [2] Bieniawski ZT. Mechanism of brittle fracture of rock, part II— experimental studies. Int J Rock Mech Min Sci Geomech Abstr 1967;4:407–23. [3] Hallbauer DK, Wagner H, Cook NGW. Some observations concerning the microscopic and mechanical behaviour of quartzite specimens in stiff, triaxial compression tests. Int J Rock Mech Min Sci Geomech Abstr 1973;10:713–26. [4] Kranz RL. Crack growth and development during creep of Barre granite. Int J Rock Mech Min Sci Geomech Abstr 1979;16:23–35. [5] Kranz RL. Crack–crack and crack–pore interactions in stressed granite. Int J Rock Mech Min Sci Geomech Abstr 1979;16:37–47. [6] Lockner D. The role of acoustic emission in the study of rock fracture. Int J Rock Mech Min Sci Geomech Abstr 1993;30(7): 883–99. [7] Aubertin M, Julien MR, Li L. The semi-brittle behavior of low porosity rocks. In: Proceedings of the Third North American Rock Mechanics Symposium, vol. 2. Mexico: International Society for Rock Mechanics; 1998. p. 65–90. [8] Blair SC, Cook NGW. Analysis of compressive fracture in rock using statistical techniques: part I. A non-linear rule-based model. Int J Rock Mech Min Sci Geomech Abstr 1998;35(7):837–48. [9] Blair SC, Cook NGW. Analysis of compressive fracture in rock using statistical techniques: part II. Effect of microscale heterogeneity on macroscopic deformation. Int J Rock Mech Min Sci Geomech Abstr 1998;35(7):849–61. [10] Kemeny JM. A model for non-linear rock deformation under compression due to sub-critical crack growth. Int J Rock Mech Min Sci Geomech Abstr 1991;28(6):459–67. [11] Kemeny JM, Cook NGW. Micromechanics of deformation in rocks. In: Shah SP, editor. Toughening mechanisms in quasibrittle materials. Dordrecht: Kluwer Academic Publishers; 1991. p. 155–88. [12] Lockner DA, Madden TR. A multiple-crack model of brittle fracture, 1, non-time-dependent simulations. J Geophys Res 1991; 96:19623–42. [13] Lockner DA, Madden TR. A multiple-crack model of brittle fracture, 2, time-dependent simulations. J Geophys Res 1991;96:19643–54. [14] Lockner D. A generalized law for brittle deformation of westerly granite. J Geophys Res 1998;103(B3):5107–23. [15] Okui Y, Horii H. Stress and time-dependent failure of brittle rocks under compression: a theoretical prediction. J Geophys Res 1997;102(B7):14869–81. [16] Cundall PA, Potyondy DO, Lee CA. Micromechanics-based models for fracture and breakout around the mine-by tunnel. In: Martino JB, Martin CD, editors. Proceedings of the Excavation Disturbed Zone Workshop, Designing the Excavation Disturbed Zone for a Nuclear Repository in Hard Rock. Toronto: Canadian Nuclear Society; 1996. p. 113–22. [17] Cundall PA, Drescher A, Strack ODL. Numerical experiments on granular assemblies, measurements and observations. In: Vermeer PA, Luger HJ, editors. Deformation and failure of granular materials. Rotterdam: Balkema; 1982. p. 355–70. [18] Krajcinovic D. Damage mechanics. Mech Mater 1989;8:117–97. [19] Herrmann HJ, Roux S. Statistical models for the fracture of disordered media. New York: North-Holland; 1990. [20] Charmet JC, Roux S, Guyon E. Disorder and fracture. New York: Plenum Press; 1989. [21] Schlangen E, Garboczi EJ. Fracture simulations of concrete using lattice models: computational aspects. Eng Frac Mech 1997; 57:319–32. [22] D’Addetta GA, Kun F, Ramm E. On the application of a discrete model to the fracture process of cohesive granular materials. Granular Matter 2002;4:77–90. [23] Donzé F, Magnier SA. Formulation of a 3-D numerical model of brittle behaviour. Geophys J Int 1995;122:790–802. [24] Handley MF. An investigation into the constitutive behaviour of brittle granular media by numerical experiment. Ph.D. thesis, University of Minnesota, USA; 1995. [25] Jirasek M, Bazant ZP. Discrete element modeling of fracture and size effect in quasibrittle materials: analysis of sea ice. In: Proceedings of the Second International Conference on Discrete Element Methods. Cambridge, MA: IESL Publications; 1993. p. 357–68. [26] Bazant ZP, Tabbara MR, Kazemi MT, Cabot GP. Random particle model for fracture of aggregate or fiber composites. J Eng Mech 1990;116(8):1686–705. [27] Trent BC, Margolin LG, Cundall PA, Gaffney ES. The micromechanics of cemented granular material. In: Constitutive laws for engineering materials: theory and applications. Amsterdam: Elsevier; 1987. p. 795–802. [28] Lorig LJ, Cundall PA. Modeling of reinforced concrete using the distinct element method. In: Proceedings of the SEM/RILEM International Conference on Fracture of Concrete and Rock. Bethel, CT: SEM; 1987. p. 459–71. [29] Plesha ME, Aifantis EC. On the modeling of rocks with microstructure. In: Rock mechanics—theory–experiment–practice. New York: Association of Engineers and Geologists; 1983. p. 27–35. [30] Zubelewicz A, Mroz Z. Numerical simulation of rockburst processes treated as problems of dynamic instability. Rock Mech 1983;16:253–74. [31] Napier JAL, Malan DF. A viscoplastic discontinuum model of time-dependent fracture and seismicity effects in brittle rock. Int J Rock Mech Min Sci Geomech Abstr 1997;34(7):1075–89. [32] Steen B, Vervoort A, Napier JAL. Numerical modelling of fracture initiation and propagation in biaxial tests on rock samples. Int J Frac 2001;108:165–91. [33] Fang Z. Harrison JP. Development of a local degradation approach to the modelling of brittle fracture in heterogeneous rocks. Int J Rock Mech Min Sci Geomech Abstr 2002;39: 443–57. [34] Tang CA. Numerical simulation of progressive rock failure and associated seismicity. Int J Rock Mech Min Sci Geomech Abstr 1997;34(2):249–61. [35] Jing L. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. Int J Rock Mech Min Sci Geomech Abstr 2003;40:283–353. [36] Ingraffea AR, Wawrzynek PA. Crack propagation. In: Encyclopedia of materials: science and technology. New York: Elsevier Science; 2001. [37] Potyondy DO, Cundall PA, Lee C. Modeling rock using bonded assemblies of circular particles. In: Aubertin M, Hassani F, Mitri H, editors. Rock mechanics tools and techniques. Rotterdam: Balkema; 1996. p. 1937–44. [38] Robertson D, Bolton MD. DEM simulations of crushable grains and soils. In: Kishino Y, editor. Powders and grains. Lisse. The Netherlands: Balkema; 2001. p. 623–6. [39] Potyondy D, Autio J. Bonded-particle simulations of the in-situ failure test at Olkiluoto. In: Elsworth D, Tinucci JP, Heasley KA, editors. Rock mechanics in the national interest, vol. 2. Lisse, The Netherlands: Balkema; 2001. p. 1553–60. [40] Autio J, Wanne T, Potyondy D. Particle mechanical simulation of the effect of schistosity on strength and deformation of hard rock. In: Hammah R, Bawden W, Curran J, Telesnicki M, editors. NARMS-TAC 2002. Toronto: University of Toronto Press; 2002. p. 275–82. ARTICLE IN PRESS D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 [41] Itasca Consulting Group Inc. PFC2D/3D (Particle Flow Code in 2/3 Dimensions), Version 2.0. Minneapolis, MN: ICG; 1999. [42] Cundall PA. A computer model for simulating progressive large scale movements in blocky rock systems. In: Proceedings of the Symposium of International Society of Rock Mechanics, vol. 1, Nancy: France; 1971. Paper No. II-8. [43] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979;29(1):47–65. [44] Cundall PA. Formulation of a three-dimensional distinct element model—part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25(3):107–16. [45] Hart R, Cundall PA, Lemos J. Formulation of a threedimensional distinct element model—part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25(3):117–25. [46] Cundall PA, Hart R. Numerical modeling of discontinua. J Eng Comp 1992;9:101–13. [47] Cook BK, Jensen RP. Discrete element methods: numerical modeling of discontinua (Proceedings of the Third International Conference on Discrete Element Methods), Geotechnical Special Publication No. 117. Reston, VA: American Society of Civil Engineers; 2002. [48] Cundall PA. Distinct element models of rock and soil structure. In: Brown ET, editor. Analytical and computational methods in engineering rock mechanics. London: Allen and Unwin; 1987. p. 129–63. [49] Hazzard JF, Young RP, Maxwell SC. Micromechanical modeling of cracking and failure in brittle rocks. J Geophys Res 2000;105(B7):16683–97. [50] Feustel AJ. Seismic attenuation in underground mines: measurement techniques and applications to site characterization. Ph.D. thesis, Queen’s University, Kingston, Ontario, Canada; 1995. [51] Hazzard JF, Young RP. Simulating acoustic emissions in bondedparticle models of rock. Int J Rock Mech Min Sci Geomech Abstr 2000;37:867–72. [52] Holt RM. Particle vs. laboratory modelling of in situ compaction. Phys Chem Earth 2001;26(1–2):89–93. [53] Potyondy D, Cundall P. The PFC model for rock: predicting rock-mass damage at the underground research laboratory. Ontario Power Generation, Nuclear Waste Management Division Report No. 06819-REP-01200-10061-R00, May 2001. [54] Oreskes N, Shrader-Frechette K, Belitz K. Verification, validation, and confirmation of numerical models in the earth sciences. Science 1994;263:641–6. [55] Martin CD. The strength of massive Lac Du Bonnet granite around underground openings. Ph.D. thesis, University of Manitoba, Winnipeg, Canada, June; 1993. [56] Hoek E, Brown ET. Underground excavations in rock. London: Institute of Mining and Metallurgy; 1980. [57] Read RS. Interpreting excavation-induced displacements around a tunnel in highly stressed granite. Ph.D. thesis, University of Manitoba, Winnipeg, Canada; 1994. [58] Brady BHG, Brown ET. Rock mechanics for underground mining. London: Allen and Unwin; 1985. [59] Mosher S, Berger RL, Anderson DE. Fracture characteristics of two granites. Rock Mech 1975;7:167–76. [60] Boutt DF, McPherson BJOL. Simulation of sedimentary rock deformation: lab-scale model calibration and parameterization. Geophys Res Lett 2002;29(4). Available from Ontario Power Generation Inc., Nuclear Waste Management Division (16th Floor), 700 University Avenue, Toronto, Ontario, M5G 1X6. 1363 [61] Everitt, RA, Lajtai EZ. The influence of rock fabric on excavation damage in Lac du Bonnet granite. Int J Rock Mech Min Sci Geomech Abstr 2004;41(8). [62] Paterson MS. Experimental rock deformation—the brittle field. Berlin: Springer; 1978. p. 35. [63] Read RS, Chandler NA. An approach to excavation design for a nuclear fuel waste repository. Ontario Power Generation, Nuclear Waste Management Division Report No. 06819-REP-0120010086-R00, 2002. [64] SAFETI. Seismic validation of 3D thermo-mechanical models for the prediction of rock damage around radioactive waste packages in geological repositories. Project for European Atomic Energy Community; 2001–2004. [65] Read RS, Martin D. Technical summary of AECL’s mine-by experiment, phase 1: excavation response. Atomic Energy of Canada Limited, Report AECL-11311, 1996. [66] Martin CD, Read RS. AECL’s mine-by experiment: a test tunnel in brittle rock. In: Aubertin M, Hassani F, Mitri H, editors. Rock mechanics tools and techniques. Rotterdam: Balkema; 1996. p. 13–24. [67] Chandler N, Dixon D, Gray M, Hara K, Cournut A, Tillerson J. The tunnel sealing experiment: an in-situ demonstration of technologies for vault sealing. In: Proceedings of the 19th Annual Conference of Canadian Nuclear Society. Toronto: Canadian Nuclear Society; 1998. [68] Read RS, Martino JB, Dzik EJ, Oliver S, Falls S, Young RP. Analysis and interpretation of AECL’s heated failure tests. Ontario Power Generation, Nuclear Waste Management Division Report No. 06819-REP-01200-0070-R00, 1997. [69] Read RS, Martino JB, Dzik EJ, Chandler NA. Excavation stability study—analysis and interpretation of results. Ontario Power Generation, Nuclear Waste Management Division Report No. 06819-REP-01200-0028-R00, 1998. [70] Martino JB. A review of excavation damage studies at the underground research laboratory and the results of the excavation damage zone study in the tunnel sealing experiment. Ontario Power Generation, Nuclear Waste Management Division Report No. 06819-REP-01200-10018-R00, 2000. [71] Anderson OL, Grew PC. Stress corrosion theory of crack propagation with applications to geophysics. Rev Geophys Space Phys 1977;15(1):77–104. [72] Atkinson BK, Meredith RG. The theory of subcritical crack growth with applications to minerals and rocks. In: Atkinson BK, editor. Fracture mechanics of rock. London: Academic Press; 1987. p. 111–66. [73] Potyondy DO, Cundall PA. Modeling notch-formation mechanisms in the URL mine-by test tunnel using bonded assemblies of circular particles. Int J Rock Mech Min Sci Geomech Abstr 1998;35(4–5) Paper No. 067. [74] Itasca Consulting Group Inc. FLAC (Fast Lagrangian Analysis of Continua), Version 4.0. Minneapolis, MN: ICG; 2000. [75] Fakhimi A, Carvalho F, Ishida T, Labuz JF. Simulation of failure around a circular opening in rock. Int J Rock Mech Min Sci Geomech Abstr 2002;39:507–15. [76] Lau JSO, Chandler NA. Innovative laboratory testing. Int J Rock Mech Min Sci Geomech Abstr 2004;41(8). [77] Johnson S. Emergence: the connected lives of ants, brains, cities, and software. New York: Simon and Schuster; 2001. [78] Anderson TL. Fracture mechanics: fundamentals and applications. Boca Raton, Florida: CRC Press; 1991. [79] Bazant ZP, Planas J. Fracture and size effect in concrete and other quasibrittle materials. Boca Raton, Florida: CRC Press; 1998. Available from AECL, Scientific Document Distribution Office (SDDO), Chalk River Laboratories, Chalk River, Ontario, K0J 1J0. ARTICLE IN PRESS 1364 D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 [80] Griffith AA. The phenomena of rupture and flow in solids. Philos Trans Roy Soc Ser A 1920;221:163–98. [81] Griffith AA. The theory of rupture. In: Biezeno CB, Burgers JM, editors. Proceedings of the First International Congress Applied Mechanics. Delft: Tech. Boekhandel en Drukkerij J Walter Jr; 1924. p. 54–63. [82] Eringen AC. Continuum mechanics at the atomic scale. Crystal Lattice Defects 1977;7:109–30. [83] Huang H. Discrete element modeling of tool-rock interaction. Ph.D. thesis, University of Minnesota, USA; 1999. [84] Zhang ZX. An empirical relation between mode I fracture toughness and the tensile strength of rock. Int J Rock Mech Min Sci Geomech Abstr 2002;39:401–6. [85] Li L, Holt RM. Simulation of flow in sandstone with fluid coupled particle model. In: Elsworth D, Tinucci JP, Heasley KA, [86] [87] [88] [89] [90] editors. Rock mechanics in the national interest, vol. 1. Lisse, The Netherlands: Balkema; 2001. p. 165–72. Itasca Consulting Group Inc. PFC2D (Particle Flow Code in 2 Dimensions), Version 3.0. In: Optional features volume, thermal option. Minneapolis, MN: ICG; 2002. Linkov AM. Keynote lecture: new geomechanical approaches to develop quantitative seismicity. In: Gibowicz SJ, Lasocki S, editors. Rockbursts and seismicity in mines. Rotterdam: Balkema; 1997. p. 151–66. Hazzard JF, Young RP. Dynamic modelling of induced seismicity. Int J Rock Mech Min Sci Geomech Abstr 2004;41(8). Fung YC. A first course in continuum mechanics. Englewood Cliffs, NJ: Prentice-Hall; 1977. Ugural AC, Fenster SK. Advanced strength and applied elasticity. 2nd SI ed. New York: Elsevier Science; 1987. p. 71.