当前位置:首页 >> 电力/水利 >>

A bonded-particle model for rock


ARTICLE IN PRESS

International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364

A bonded-particle model for rock
D.O. Potyondya,?,1, P.A. Cundall

b
a

www.elsevier.com/locate/ijrmms

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 con?nement. 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 in?uences 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 ?delity. 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: david.potyondy@utoronto.ca (D.O. Potyondy), itasca@itascacg.com (P.A. Cundall). 1 Formerly with Itasca.

thus exhibits a rich set of emergent behaviors that correspond very well with those of real rock. The BPM provides both a scienti?c 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

1365-1609/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2004.09.011

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 uncon?ned or con?ned 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 ?lled with cement). All of these items in?uence 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 ?lled with cement. The force chains are similar to those that form in a granular material [17]. The cement-?lled 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 1331

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 uncon?ned 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 classi?ed 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.

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 dif?cult 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 ?eld-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 identi?ed. 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 ?eld 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.

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 ?nite 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 de?ned in Ref. [46] as methods that allow ?nite 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

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 de?nition in the ?eld 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 ?nite amount of space. The rigid particles interact only at the soft contacts, which possess ?nite 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 ?nite mass. 2. The particles move independently of one another and can both translate and rotate.

ARTICLE IN PRESS
D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1333

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 ?nite-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 dif?culty, 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 coef?cient, 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)

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 de?ned 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 coef?cient 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 coef?cient corresponding with a quality factor of 100 and found that the energy released by crack events (i.e., bond breakages) may have a signi?cant in?uence 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 uncon?ned 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, ?rst 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 coef?cient, 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? in Fig. 3), and two i 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 s

increment of elastic shear force is given by DF s ? ?ks DU s ; where ks is the contact shear stiffness given by ks ? k?A? k?B? s s k?A? ? k?B? s s ; (8) (7)

(4)

where F and F denote the normal and shear force components, respectively, and ni and ti are the unit vectors that de?ne 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 ? k?A? k?B? n n
?B? k?A? ? kn n k?A? and k?B? n n

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     _i _i V i ? x?c? ? x?c? B  A  ?B? ?B? _ ? xi ? eijk oj x?c? ? x?B? k k    ?A? ?A? ?c? _ ? xi ? eijk oj xk ? x?A? ; ?9? k _ xi and oj are the particle translational and rotational velocities, respectively, and eijk is the permutation symbol. The relative shear–displacement increment vector is DU s ? V s Dt ? ?V i ? V n ?Dt ? ?V i ? V j nj ni ?Dt: i i i (10)

;

(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

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 coef?cient ? ? m ? min m?A? ; m?B? ; (11) with m?A? and m?B? being the particle friction coef?cients, 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 1335

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 ?ve parameters that de?ne a parallel bond (see Fig. 3): normal and shear stiffnesses per unit area, ?n ?s k and k ; tensile and shear strengths, sc and tc ; and ? ? ? bond-radius multiplier, l; such that parallel-bond radius ? ? ? ? R ? l min R?A? ; R?B? ; (12)

are given by ? n ?n DF ? k A DU n ; ?s ?s DF ? ?k A DU s ; ?s ?n DM ? ?k J Dyn ; ? ? DM ? ?k I Dys ;
s n

?14?

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 t i ;

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 ( ? 2Rt; t ? 1; PFC2D; A? ?2 pR ; PFC3D; 8 2?3 < R t; t ? 1; PFC2D; 3 I? :1 ?4 pR ; PFC3D; (4 NA; PFC2D; J? ?15? 1 ?4 PFC3D: 2 pR ; The maximum tensile and shear stresses acting on the parallel-bond periphery are calculated from beam theory to be ? ? ? ?F jM jR ; ? ? s ? I A s n ? ? ? jF j jM jR ? : ?16? tmax ? ? A J If the maximum tensile stress exceeds the tensile strength ?smax Xsc ? or the maximum shear stress exceeds ? ? the shear strength ?? max X? c ?; then the parallel bond t t breaks, and it is removed from the model along with its accompanying force, moment and stiffnesses. A simpli?ed 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 de?ned 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 n s

?13?

?n ?s ?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 de?ne 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 rotation   increment ?DU n ; DU s ; Dyn ;  ?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

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 in?uences 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 ?xed 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 identi?es Dmin as the controlling length scale of the material. The ?nal item that characterizes the BPM is the grain-cement microproperties: ? ? ? ? E c ; kn =ks ; m ; grain microproperties; ? ? n s? ? ? ? ? ? ? ? l; E c ; k =k ; sc ; tc ; cement microproperties; ?17? ? where E c and E 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 coef?cient; and sc and tc 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 ? ? Ec ?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 con?rmed 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 speci?ed directly. It is possible, however, to ? relate the grain and cement moduli, E c and E 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 ? k?A? ? k?B? in Eq. (6). For the n n cement-based behavior, ? ? AE c AE c ?n ? ?A? k A? L R ? R?B? ? n ? ?A? ? ? ) E c ? k R ? R?B? :

?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 ?xed 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 ?ve-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 ?ll 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) ?oating 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 Cement

? l?1 ? Ec ?

62 GPa; 72 GPa; ?n ?s ?k =k ? ? 2:5

&

PFC2D PFC3D & 157 ? 36 MPa; 175 ? 40 MPa; PFC2D PFC3D

sc ? tc ? ?mean ? std: dev:?? ? ?

vessel is 16% for PFC2D and 35% for PFC3D. The particles, at half their ?nal size, are placed randomly such that no two particles overlap. Then, the particle radii are increased to their ?nal 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 speci?ed isotropic stress: The radii of all particles are reduced uniformly to achieve a speci?ed isotropic stress, s0 ; de?ned 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 ‘‘?oating’’ 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 ‘‘?oating’’ 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 ?oating 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 ?oating particles

to zero, we obtain a bonded assembly for which nearly all particles away from the specimen boundaries have at least three contacts. The ?oaterelimination 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 10?6 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 sc and ? tc 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 signi?cant 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 1339

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 veri?able [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 ?eld 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 ?lls 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 (de?ned in Section 3.3), and the mean value of the cement strengths is chosen to match the uncon?ned compressive strength. The parti-

cle-friction coef?cient 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 sc ? tc ; both tensile and shear microfai? ? lures are possible. If microtensile failure is excluded (by setting sc to in?nity), 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 tc to in?nity) ? or allowed (by setting sc ? tc ) are investigated for the ? ? PFC2D material by Potyondy and Autio [39], who study damage formation adjacent to a circular hole subjected to far-?eld compressive loading. Both materials produce breakout notches in compressive regions and tensile macrofractures in tensile regions that are generally similar. The material with sc ? tc allows ? ? more well-de?ned 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, sc ? tc 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 ?gures 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 con?nements 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 1341

and do not seem to be interacting. The onset of this early cracking is quanti?ed 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 con?nement 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 con?nement 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 con?nements, 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 con?nement than is the pre-peak damage-formation process. This change in behavior is related in part to the fact that, as con?nement increases, the ratio of shear microcracks to tensile microcracks increases. The con?nement 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 dif?culties are inherent in tests that attempt to measure the tensile strength of a rock, because it is rarely possible to force the failure to

occur simultaneously across the entire ?nal 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 uncon?ned compression tests; (b) the uncon?ned compressive strength, qu ; (c) the crack-initiation stress, sci ; (d) the strength envelope for con?ning 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 de?ned as the axial stress at which non-elastic dilation just begins, identi?ed 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 con?ning pressure, is obtained by ?tting the Hoek–Brown strength criterion [56] to results of triaxial tests at con?ning pressures up to 70 MPa. The Hoek–Brown strength criterion is given by q????????????????????????? sf ? s3 ? ss2 ? s3 sc m (21) c where s and m are dimensionless material parameters and sc is equal to the uncon?ned 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 con?nement; however, it can be linearized using a secant approximation de?ned by the strength, sf ; at two con?nements, P0 and P1 : If the slope is de?ned by Nf ? sf ?P1 ? ? sf ?P0 ? ; P1 ? P0 (22)

Table 2 Macroproperties of the models and Lac du Bonnet granite Property E ?GPa? u qu ?MPa? sci ?MPa? scd ?MPa? Nf f ?deg? c ?MPa? st ?MPa? Lac du Bonnet granite 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? PFC2D model, Davg ? 0:72 mm 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? PFC3D model, Davg ? 1:53 mm 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 ? pu ; ?23? 2 Nf where qu is the uncon?ned 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 de?ned 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 ? where x is the con?ning 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 con?ning 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 con?ning 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.

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).

of these macroproperties (expressed as a ratio of standard deviation to mean) is less than that of the physical material. The uncon?ned compressive strength of the PFC3D material exhibits less variability than that of the PFC2D material. The model materials also exhibit the onset of signi?cant 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 ?nitestrength 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 1343

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 de?ned as a set of particles that are adjacent to

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 identi?es 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 ?nite 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 If the BPM is used to predict damage formation surrounding an excavation, then the absolute size of the potential damage region is ?xed, and successively ?ner resolution models can be produced by decreasing Dmin while keeping all other microproperties ?xed. The investigation in this section demonstrates that particle size is an intrinsic part of the material characterization that affects the Brazilian strength (and the uncon?ned 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

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.

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 con?nements 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 coef?cient of variation (ratio of standard deviation to mean) of each macroproperty. It is expected that as particle size continues to decrease, the coef?cients of variation will converge to speci?c 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 suf?cient number of particles must be present for the model to resolve and reproduce the failure mechanisms that in?uence 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 ; E 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 uncon?ned 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 coef?cients of variation are reasonable and indicate a possible convergence to a speci?c 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 signi?cant damage has developed. No de?nitive 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) Macroproperty (63.4 ? 31.7 mm2 specimens, n=10, mean and coef?cient of variation) E (GPa, %) 2.87 1.44 0.72 0.36 68.3 68.8 70.9 71.5 10.5 3.3 1.3 0.8 n (?, %) 0.231 0.249 0.237 0.245 19.0 7.6 4.6 2.9 qu (MPa, %) 186.8 184.4 199.1 194.8 12.7 7.9 6.5 4.1 f (deg, %) 34.9 30.3 29.5 33.0 18.1 23.8 16.3 17.3 c (MPa, %) 48.5 53.4 58.5 53.3 12.8 18.9 14.5 14.4 st (MPa, %) 65.5 54.3 44.7 35.4 26.1 12.2 7.4 7.6

Table 4 Effect of particle size on PFC3D macroproperties Particle size Davg (mm) Macroproperty (63.4 ? 31.7 ? 31.7 mm3 specimens, n=10, mean and coef?cient of variation) E (GPa, %) 5.95 3.05 2.04 1.53 57.3 64.0 67.6 69.2 10.0 3.9 1.8 1.2 n (?, %) 0.231 0.254 0.255 0.256 21.2 8.1 5.8 5.5 qu (MPa, %) 127.9 169.6 186.9 198.8 11.9 3.4 1.5 3.6 f (deg, %) 25.9 30.6 32.3 32.1 14.2 9.8 9.2 7.4 c (MPa, %) 40.0 48.4 51.6 55.1 11.4 7.6 6.9 7.6 st (MPa, %) 43.6 35.4 33.0 27.1 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 1345

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 coef?cients 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 uncon?ned 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 uncon?ned 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 micropro?n ?s ? perties in terms of E c ; E 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 uncon?ned 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 coef?cients of variation have converged to approximately 3.5%. The reason for this size effect is unknown. It corresponds with general observations that ?ner-grained rock is stronger than coarser-grained rock (medium and ?ne-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 de?nitive 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 coef?cients of variation have not yet

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 uncon?ned 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 quanti?ed using a ‘‘strain probe’’, which applies a speci?ed 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 modi?es 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 identi?ed 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 ?eld 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 (?eld 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 ?oor (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 ?oor. This induced stress is lower than the 200 MPa uncon?ned 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 suf?cient 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 simpli?ed 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 simpli?ed BPM by reducing the contact-bond strength at a speci?ed rate if the contact force was above a speci?ed 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 speci?ed rate if the maximum tensile stress in the parallel bond is above a speci?ed 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 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.

Fig. 12. Photograph of broken and buckled rock slabs comprising the notch tip in the ?oor of the Mine-by Experiment tunnel.

ARTICLE IN PRESS
D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1347

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 ?oor 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 speci?ed, 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 ?nal 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

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 speci?ed 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 ?eldscale simulations using a uniform strength-reduction factor are presented, whereas the stress-corrosion model is described elsewhere [53] and is only brie?y 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-?eld 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-?eld 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 simpli?ed 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 signi?cant 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 ?nal 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 1349

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 uncon?ned 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, uncon?ned- and con?nedcompressive 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

behaviors, using few microparameters (but at the expense of often quite lengthy calibration procedures). Further, continuum-based models (typically implemented via ?nite or boundary element techniques) have dif?culty 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 ?nite 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 in?nitely 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; p?????? (29) F max ? 2sf t aR: n Noting that the mode-I stress ? intensity factor for the p????? system is de?ned as K I ? sf pa; Eq. (29) can be written as r???? R max F n ? 2tK I : (30) 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), p??????? (31) K Ic ? s0t pR: 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: p??????????????????? p????????????????????? 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-?eld load remains constant, but proximate boundaries or ?xed-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 in?uence of the complete model on the contact forces near the crack tip.

Fig. 16. Force distribution and displacement ?eld near the crack tip in a cracked cubically packed contact-bonded disk assembly.

Fig. 17. Normalized tensile force of LEFM and BPM at the ?ve contacts ahead of the crack.

the crack, the induced stress, sn ; acting on the crack plane near the crack tip ?r ( a? is [78] r????? 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 r??? Z r?R Z r?R a Fn ? sn t dr ? sf t r?1=2 dr 2 r?R r?R p?????? p???? p???????????? ? 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 ~ Fn ? p???? p???????????? Fn p?????? ? m ? m ? 1 2sf t aR (28)

ARTICLE IN PRESS
D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1351

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)

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 p????????? ?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 p??????? ? 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 ?ve 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 parallel? ?n ?s ?? bonded material (with l ? 1; k =k ? 2:5 and s0t ? l sc ) are also shown in Fig. 18. The parallel-bonded material

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; ?xed. For LEFM conditions, the applied far-?eld tensile stress at incipient failure, s0 ; is [78] K Ic p??? a?1=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 p????????? c?1=2 : ? 0 (35) 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 p??? p??????? K Ic ? ?10b C?f? 2?s0t 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-?fth and one-tenth

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 p????????? 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 p??? K Ic b ? 0 p??????? ? 10b C?1=3? 2 st pR p??? ?40? ? 10?0:3623 C?1=3? 2 ? 0:66: 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 Grif?th criterion). An extensive comparison between BPMs and the results of fracture mechanics is given in Ref. [83]. 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

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 p??? 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 ? 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).

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 1353

which can be substituted into Eq. (41) to give K Ic st / p???? ; 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 p???? K Ic / sc R; (43) ? where sc ? tc is the mean parallel-bond strength with ? ? ? l ? 1; and R is the mean particle radius. Combining Eqs. (42) and (43) yields r???? R st / sc (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 de?ned 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 ?xed, then the same Brazilian strength will be measured by either doubling the size of the Brazilian disk for a ?xed particle size, or halving the particle size for a ?xed Brazilian disk size. The BPM for granite exhibits such dual-model similarity as shown in Fig. 21.

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 =sc ? versus ?R=D? ? on a log–log scale. The test data are well ?t 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 ? 0:2D=2R? is relatively small ranging from 2.2 to 17.6. Eq. (44) can be expressed as st / c?1=2 (45) sc ? 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 signi?cant 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 uncon?ned 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.

Fig. 21. Dual-model similarity: If all microproperties are kept ?xed, then the exact same six short-term macroproperties in Table 3 are measured for both models B and C.

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 con?ning 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 de?ne 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 ?uid ?ow [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 signi?cant 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 modi?cations 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 ?rst 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 ?fth EURATOM framework program, Nuclear Fission (1998–2002). Appendix A All vector and tensor quantities are expressed using indicial notation with respect to a ?xed 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., _ xi ? @xi =@t). The Kronecker delta, dij ; and the permutation symbol, eijk ; are employed. A.1. Stress-measurement procedure The average stress, sij ; in a volume, V ; of material is ? de?ned by Z 1 sij ? sij dV (46) ? V V

ARTICLE IN PRESS
D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1355

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? sij V ; sij ? ? ? V N
p

centroid to the contact location. By substituting Eqs. (53) into (52) and noting that X ?c? Fj ? 0 (54)
Nc

(47)

where s?p? is the average stress in particle ?p?: In the ? ij same way, the average stress in a particle ?p? can be written using Eq. (46) as Z 1 ?p? sij ? ?p? s?p? dV ?p? : (48) ? ij ?p? V V The identity S ij ? dik S kj ? xi;k Skj ? ?xi Skj ?;k ? xi Skj;k (49)

for a particle in equilibrium, one obtains  1 X ?c?  ?p?  ?c;p? ?c? s?p? ? ? ?p? ? ij xi ? xi ni F j : V Nc

(55)

An expression for the average stress in a volume, V ; is obtained by substituting Eq. (55) into (47); however, de?nition 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) 1?n 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? sij ? ?@P ?p? A (57) ? xi ? xi ni F j ; V Np Nc V?
Np

holds for any tensor S ij : Applying this identity to the stress in a particle, one can write ! Z   1 ?p? ?p? ?p? sij ? ?p? xi skj ? xi skj;k dV ?p? : (50) ? ;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 ?rst 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? sij ? ?p? ? kj ?p? V ZS 1 ? ?p? xi t?p? dS?p? ; ?51? j ?p? V S where S ?p? is the particle surface, nk is the unit outward normal to the surface and t?p? is the traction vector. If j 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)

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? are the locations of a particle i i 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? is the force acting at contact j ? ?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 ?eld in the voids, a velocity-gradient tensor, based on a best-?t 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.

where x?c? is the location and F ?c? is the force acting at i j ? contact ?c?: F ?c? includes both F i and F i from Eqs. (4) j 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     x?c? ? x?p? ? x?c? ? x?p? n?c;p? ; (53) i i i i i where x?p? is the location of the particle centroid and i n?c;p? is the unit-normal vector directed from the particle i

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 de?ned. 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 and oij ? 1 ?uj;i ? ui;j ?; rotation tensor: 2 (59) where eij ? 1?ui;j ? uj;i ?; 2 infinitesimal-strain tensor

A measure of the error in these predicted values is given by 2 X   X  ?p?  ~ ?p?  ~ ?p? v?p? ? V ?p? ; ~ ~ ~i ~i v?p? ? V i z?  vi ? V 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 : @_ ij a (66)

In a similar fashion, the relation between velocities vi at two neighboring points is given by the velocity_ gradient tensor, aij : (The velocity ?eld is related to the displacement ?eld by ui ? vi dt; in which vi is the velocity and dt is an in?nitesimal 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 ? aij dxj : (60) The velocity-gradient tensor can be decomposed into a symmetric and an anti-symmetric tensor as _ _ _ aij ? eij ? oij where _ eij ? 1 ?vi;j ? vj;i ?; strain-rate tensor 2 and ? ? _ oij ? 1 vj;i ? vi;j ; 2 spin tensor: (61)

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 ~ ~ ~ ~ ~ ~ x1 x1 x2 x1 x3 x1 Np Np 78 _ 9 6 Np 7> ai1 > 6 < = 6 P ?p? ?p? P ?p? ?p? P ?p? ?p? 7> > 6 x1 x2 ~ ~ ~ ~ ~ ~ 7 _ x2 x2 x3 x2 7 a 6N i2 Np Np 7> > 6 p : ; 7> > 6 _ 4 P ?p? ?p? P ?p? ?p? P ?p? ?p? 5 ai3 ~ ~ ~ ~ ~ ~ x1 x3 x2 x3 x3 x3
Np

9 8 P ?p? ~ ~ > > V i x?p? > > 1 > >N > > p > > > > > > > >P < ?p? ?p? = ~ x V i ~2 : ? > > Np > > > > > >P > > > > V ?p? x?p? > ~ ~ > > ; : i 3 >
Np

Np

Np

?67?

_ aij is computed using the following least-squares procedure. The velocity-gradient tensor computed over a measurement region represents the best ?t 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 Np Np ? Vi ? and xi ? ; (62) ? Np Np where V ?p? and x?p? are the translational velocity and i i centroid location, respectively, of particle ?p?: The measured relative velocity values are given by ~ ?p? Vi ? V ?p? i ? ? V i: (63)

These nine equations are solved by performing a single LU decomposition upon the 3 ? 3 coef?cient 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 aij ? @ui ? ui;j ? eij ? oij @xj ? ? ? ? ? 1 ui;j ? uj;i ? 1 uj;i ? ui;j ; 2 2

_ ~i For a given aij ; the predicted relative velocity values v?p? can be written, using Eq. (60), as _ ~j _ j ~i v?p? ? aij x?p? ? aij ?x?p? ? xj ?: ? (64)

?68?

where eij is the in?nitesimal-strain tensor and oij is the _ rotation tensor. The velocity-gradient tensor, aij ; relates _ the velocities vi ? ui 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 1357

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 ? aij 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 xj be a ?xed reference position, then ? Z xj Z xj Z xj @ _ ij dxj ? ?eij ? oij ? dxj ; dvi ? a @t ? ? ? xj xj xj Z xj _ ?70? vi ?xj ? ? vi ? eij dxj ; ?
? xj

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 _ eij ? Deij ; (71) NDt where Deij is a strain increment applied over N time steps, each of size Dt; then the velocity ?eld can be written as   Deij (72) vi ?xj ? ? vi ? ?xj ? xj ? ? ? NDt by assuming that the strain ?eld does not vary with position. (If the strain ?eld did vary with position, one could express the variation explicitly before performing the integration and thus obtain a ?nal expression for the velocity that would depend upon the parameters used to characterize the strain-?eld 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 ?xed, 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 ?eld is measured by averaging the stresses within a set of measurement regions that cover the ensemble. The stress increment is de?ned as the difference between the target stress and the current stress Dsij ? st ? sij ij (73)

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 ?eld 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 ?eld 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 ?eld, 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 ?nal equilibrium; instead, the entire system is made to correspond with the required strain ?eld instantaneously during each iteration); and (2) produce a stress ?eld that is nearly uniform throughout the ensemble. Such would not be the case if either the boundary particles alone or a set of con?ning walls were moved, because then non-uniform force chains would develop throughout the ensemble, producing a locally non-uniform stress ?eld—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 ?eld 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 ?ll a given area, A; or volume, V ; at a given porosity, n: The basic algorithm is presented ?rst. It makes use of two relations that are derived at the end of this section. The number of particles, N; that satis?es the above constraints is given by 8 > A?1 ? n? ; > PFC2D < Rmin ? Rmax ?2 pR ? with R ? N? : > 3V ?1 ? n? 2 > : ; PFC3D ?3 4pR (75)

and the corresponding strain increment is found by assuming isotropic, elastic behavior:   1?n n Deij ? (74) Dsij ? Dsaa dij : E E The applied particle velocities are given by Eq. (72), where xj and vi are chosen as the initial position and ? ? velocity of the ensemble centroid, respectively. Good

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 1 ?Rmin ; Rmax ? and 2 placed randomly within the speci?ed 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 satis?ed. 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; l? ?76? 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 de?ned as 8 > A ? Ap ; PFC2D; < A (77) n? > V ? Vp : ; PFC3D; V where Ap and V p are the total area and volume, respectively, of all particles given by X Ap ? pR2 ; X4 pR3 ; Vp ? ?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? : ?79? R3 ? 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 1?n R ?P l: (80) 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, n: P 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 ? NR ;

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 de?ned as the average of the direct stresses & 2; PFC2D; skk ? where l ? (82) s0 ? l 3; PFC3D: 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, sij ; in a volume, V ; of ? material 1 X X ~ ?c;p? ?c;p? ?c? sij ? ? (83) R ni F j ; ? V N N p c   ~ ?c;p? ? x?c? ? x?p?  The isotropic stress can then where R  i i : be evaluated by 1 X X ~ ?c;p? n?c? skk ? ?? s0 ? R F (84) lV N N l
p 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)

n?c?

and it is assumed that all particles remain ?xed, 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
Np Nc

lV Ds0 : ~ ?c;p? R K n?c? f?c?

(88)

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 1359

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 ‘‘?oaters’’, because they are detached from the material matrix and appear to be ?oating in space. In physical specimens of an unbonded, granular material, such as sand, ?oaters 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 ?oater 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 ?oaters is described here. It may be argued that a ‘‘solid’’ consisting of bonded circular or spherical particles, even with no ?oaters, 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 ?nal stress level to be carried by the material. The algorithm expands and moves ?oaters until every particle has the speci?ed minimum number of contacts. The key to success (i.e., elimination of all ?oaters) 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 ?nal state would be required (and ?oater elimination made more dif?cult). 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 ?oater-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 ?nal forces at levels comparable to their prebonding levels. In summary, the algorithm is described as follows. All particles except ?oaters are ?xed and given zero velocities. Floaters are then expanded by a large amount (30%), suf?cient to ensure contact with all of their

immediate neighbors. After being cycled to local equilibrium, the ?oaters 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 ?oaters become inactive. After executing this iteration, the number of active ?oaters is reduced considerably. The complete procedure is repeated several times until all ?oaters have been removed. The algorithm works well, because the adjustments to ?oater radii (and positions, through the law of motion) are done for small groups within a ‘‘cage’’ of ?xed particles. As mentioned, it is only necessary to establish equilibrium within this trapped group. Then, as ?oaters themselves become part of the cage, the scheme becomes increasingly more ef?cient as the number of particles in each ?oater group decreases. The following parameters control the algorithm: N f (minimum number of contacts to be a non-?oater, 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 ? p (mean contact normal force for single assembly); F particle); R (particle radius); and kn (particle normal stiffness). The algorithm is de?ned 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 ?oaters: a ?oater is a particle with less than N f contacts. If no ?oaters, then EXIT outer loop. Fix all non-?oaters. For all ?oaters: –expand by M r , –set friction to zero. Cycle 200 steps to get ?oater groups to equilibrium. Perform this inner loop up to 100 times. De?ne ‘‘active’’ ?oaters. An active ?oater is a particle with more than one contact and mean contact normal force ? greater than f m F a : If no active ?oaters, then EXIT inner loop. Shrink each active ?oater 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 ?xed. The inhibit ?ags are removed upon exit from the algorithm. For a large assembly, it is found that one or two ?oaters may be positioned almost exactly between two non-?oaters 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 ?oaters 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? ? Fp ? F (90) Nc N
c

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 con?ned 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 speci?ed con?ning 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 con?ned by ?uid 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 con?nements of 0.1–70 MPa) of the average particle normal stiffness to simulate a soft con?nement. 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 con?rmed 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 ?rst ?lling 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

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 ?xed during the radius change: DF n ? 1 kn DR; 2 (91)

n?c?

where kn is the particle normal stiffness. Thus, the new ?0 mean normal force, F p ; is ? 1 X ? n?c? ?0 (92) F ? DF n : Fp ? 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

ARTICLE IN PRESS
D.O. Potyondy, P.A. Cundall / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1329–1364 1361

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)

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 n? n0 ?PFC2D material; plane strain?; 1 ? n0 E ? E 0 ?1 ? n2 ?:

?98?

and the Poisson’s ratio is computed by
1 ?Dx ? Dz ? Dx u? ? ? ?2 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, ?rst the Poisson’s ratio and Young’s modulus corresponding with a state of plane stress are computed as Dx Dy Dsy E0 ? : Dy n0 ? ? ?PFC2D material; plane stress?; ?97?

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-de?ned 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 ?rst crack initiates. However, the axial stress at the time of the initiation of the ?rst 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 de?ned as the axial stress at which there exists a speci?ed 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 signi?cant 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 [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] Donze 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 ?ber 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.

References
[1] Okui Y, Horii H. A micromechanics-based continuum theory for microcracking localization of rocks under compression. In: Muhlhaus 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.

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. Geotechnique 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. Veri?cation, validation, and con?rmation 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 in?uence 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 ?eld. 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, Scienti?c Document Distribution Of?ce (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 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 ?rst 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. [80] Grif?th AA. The phenomena of rupture and ?ow in solids. Philos Trans Roy Soc Ser A 1920;221:163–98. [81] Grif?th 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 ?ow in sandstone with ?uid coupled particle model. In: Elsworth D, Tinucci JP, Heasley KA,

[86]

[87]

[88] [89] [90]


相关文章:
Contact Models
To use, enable the bonding model for both particle-particle and particle-geometry interactions. Although particles cannot be bonded to geometry elements, the...
美国药典色谱柱型号对照
L30 Ethyl silane chemically bonded to a totally porous silica particle, 3 to 10 ?m in diameter. L31 A strong anion-exchange resin-quaternary amine ...
新混凝土分析
(H2 PO4)3, a group IIA metal bonded to ...Aggregates, such as gravel or trap rock are ...However, other particle sizes, such as positive ...
专业英语A
2Gradation of particle sizes of the aggregate ...Multiple layers are bonded by adhesives to make ...Stones must be cut from the solid rock in a ...
国外标准索引 屋面 墙面工程
1998 Code of practice for the selection and application of particleboard, oriented strand board (OSB), cement bonded particleboard and wood fibreboards ...
HPLC色谱柱选择
units Ethyl silane chemically bonded to a totally porous silica particle, 3 to 10um in diameter. A strong anion-exchange resin-quaternary amine bonded on...
Contact Models
Although particles cannot be bonded to geometry elements, the bonded model ensures that the particles contact the geometry based on the particle’s physical...
Particle Replacement Tutorial
Particles do not become bonded to the geometry when using the Hertz-Mindlin with Bonding model: adding a particle-geometry interaction ensures the difference...
更多相关标签:
model for murder2016 | model for murder | v for v model | model for | vue v for v model | editorformodel | html.editorformodel | sas for mixed model |