3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing

E.Ghzvinin,M.S.Diederichs,R.Quey

aGeoEngineering Centre,Queen’s University,Kingston,ON,Canada

bEcole Nationale Supérieure des Mines de Saint-Etienne,CNRS UMR 5307,France

3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing

E.Ghazviniana,*,M.S.Diederichsa,R.Queyb

aGeoEngineering Centre,Queen’s University,Kingston,ON,Canada

bEcole Nationale Supérieure des Mines de Saint-Etienne,CNRS UMR 5307,France

A R T I C L E I N F O

Article history:

Received 17 July 2014

Received in revised form

16 September 2014

Accepted 16 September 2014

Available online 6 October 2014

Numerical modelling

3D Voronoi tessellation

Discrete element method

Grain-based model

Crack damage thresholds

Fabric-guided micro-fracturing

Anisotropy

A grain-based distinct element model featuring three-dimensional(3D)Voronoi tessellations(random poly-crystals)is proposed for simulation of crack damage development in brittle rocks.The grain boundaries in poly-crystal structure produced by Voronoi tessellations can represent f l aws in intact rock and allow for numerical replication of crack damage progression through initiation and propagation of micro-fractures along grain boundaries.The Voronoi modelling scheme has been used widely in the past for brittle fracture simulation of rock materials.However the dif f i culty of generating 3D Voronoi models has limited its application to two-dimensional(2D)codes.The proposed approach is implemented in Neper,an open-source engine for generation of 3D Voronoi grains,to generate block geometry f i les that can be read directly into 3DEC.A series of Uncon f i ned Compressive Strength(UCS)tests are simulated in 3DEC to verify the proposed methodology for 3D simulation of brittle fractures and to investigate the relationship between each micro-parameter and the model’s macro-response.The possibility of numerical replication of the classical U-shape strength curve for anisotropic rocks is also investigated in numerical UCS tests by using complex-shaped(elongated)grains that are cemented to one another along their adjoining sides.A micro-parameter calibration procedure is established for 3D Voronoi models for accurate replication of the mechanical behaviour of isotropic and anisotropic(containing a fabric)rocks. ?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

The simulation of micro-fracturing and crack damage progression for brittle rocks can be performed implicitly in continuum or explicitly with the aid of discontinuum numerical approaches.In implicit simulation of micro-crack formation,the weakness caused by the formation of cracks is smeared within the material by means of constitutive relationships.The Damage Initiation-Spalling Limit (DISL)approach is an example of implicit brittle fracture simulation with continuum codes(Diederichs,2007).The explicit simulation of micro-fracturing in rock-like material can be accomplished by direct representation of cracks in models formulated based on the Discrete Element Method(DEM)(Cundall and Hart,1985;Cundall, 1988)or by using the hybrid Finite-Discrete Element Method (FDEM)codes(Mahabadi et al.,2012).

In DEM,the rock-like material can be simulated as a dense assembly of rigid or deformable particles that interact at their contact points.Also,the discrete bodies can detach and new contacts can automatically be detected(Itasca,2013a)and therefore fractures can be simulated at the grain(or block)boundary opening.In the Bonded-Particle Model(BPM),which falls within the DEM formulation,the particles are represented by rigid disks in 2D and rigid spheres in 3D(Potyondy and Cundall,2004).The breakage of bonds between the spheres in the BPM is comparable to fracturing in rock.The BPM which is implemented in Particle Flow Code(PFC)(Itasca,2008)has been used widely for crack damage simulation in rock.The concerns and possible solutions for the accuracy of PFC simulations have been documented in the past by many researchers,including Diederichs(2003),Potyondy and Cundall(2004),Cho et al.(2007),Yoon et al.(2008),and Ghazvinian(2010).

The random polygonal blocks(grains)in DEM are an alternative geometry to the disks and spheres employed in BPM approach for simulation of fracturing.Voronoi tessellation is one of the available techniques for generating the random polygonal grains within a domain.In this technique,a region is populated with random seed points.Lines or planar surfaces are generated so that the bounded region surrounding each seed point includes all the space that is closer to that seed point than any other.The grain boundaries in thepoly-crystal structure produced by Voronoi tessellation can be used to represent f l aws in intact rock and therefore allow for simulation of crack damage development through initiation and propagation of fractures along grain boundaries.

In the Voronoi modelling approach for numerical simulation of fracturing in intact rock,the grains(blocks)can be rigid,deformable or inelastic.The behaviour of contacts between the grains is commonly governed by the available constitutive formulations for rock mass joints and discontinuities.Some of the concerns associated with the PFC results were demonstrated to be no longer applicable to the Voronoi method(Kazerani and Zhao,2010;Lan et al.,2010).The Voronoi scheme has been used widely in the past for brittle fracture modelling of rock materials,however the complexity of generating 3D Voronoi domains has limited its application to 2D models.In this paper,a new methodology for generation of 3D Voronoi Grain-Based Models(GBMs)is proposed which would provide a means for simulation of rock brittle fracturing in 3DEC(Itasca,2013a).

Damage process in laminated intact rocks or rocks that have a preferred fabric orientation,which commonly represent an anisotropic mechanical behaviour,can also be simulated by means of an equivalent continuum such as the Ubiquitous Joint Model in FLAC (Itasca,2012),or with DEM and hybrid FDEM discontinuum approaches.The discontinuum methods either apply implicit anisotropic contact and material constitutive behaviour or introduce an anisotropic defect or pre-damage to the model(e.g.Diederichs et al.,2004;You et al.,2011;Kazerani,2013;Lisjak et al.,2014a). By using the proposed approach in this paper for generation of 3D Voronoi GBMs,an analogue for intact rocks with inherent anisotropy is suggested by scaling the models and therefore achieving elongated grains.The fabric-guided micro-fracturing in the anisotropic models is veri f i ed by laboratory test results of the Cobourg Limestone and the calibration procedure for the laminated GBMs is established from sensitivity analyses of the contact microparameters.

2.Numerical approach

The extensile“crack opening forces”at the micro-scale can be generated through different mechanisms.These mechanisms amplify the localized shear strains and subsequently facilitate the nucleation of inter-granular and intra-granular cracks at the Crack Initiation(CI)threshold(Diederichs,2003).Some of the mechanisms of extension crack generation under deviatoric stress conditions are illustrated in Fig.1.

To simulate the extensile crack opening forces at the micro-scale in DEM models,the rigid or deformable particles are approximated by simple geometries.The three commonly used geometries are disks(spheres in 3D)(Potyondy and Cundall,2004;Scholtès and Donzé,2013),Voronoi grains(Damjanac et al.,2007),and triangles(trigons in 3D)(Gao and Stead,2014).The geometry of the particles plays a key role in the generation of the crack extension and shearing forces acting on the sliding fractures.A schematic of an identical micro-mechanism simulated with Voronoi grains and circular disks is illustrated in Fig.2.The applied compressive forces in the Voronoi model can be resolved into extensile and shearing forces acting on the grain boundaries(sliding f l aw in this case).In the circular disk model,the applied forces will translate into tensile stresses at the contact between the disks and rotational moments acting upon the disks.The geometrically imposed restriction on the generation of shear stresses at the contacts between the disks in the BPM has a minor effect on the CI threshold since most of the microcracks in brittle rocks form in tension(Diederichs,2007).However, with increasing deviatoric stress and on reaching the Critical Damage(CD)threshold,micro-cracks start to interact and coalesce, mostly through shearing.Therefore,the BPM with circular disks has limited application for simulation of crack interaction(Diederichs, 2007).

The simplicity of generating triangular grains(trigons)in DEM for simulation of fracturing in rock is assisting the increasing rate of using this approach.The micro-properties of the trigons models can be calibrated to simulate the exact macro-properties of rock tested in the laboratory.At the micro-scale,in contrary to the BPM models,most of the failure occurs in shear(Gao and Stead,2014).In Fig.3a,the potential sliding path in a trigon model is shown(red lines).In Fig.3b,a Voronoi tessellated model with similar edge length to the trigon model is also illustrated. While the trigon model includes smooth pathways that encourage shear sliding,the grain boundaries of the Voronoi model provide asperities and rough failure paths which lead to extensile opening of some grain boundaries and simultaneous shear sliding along other grain edges.This f i gure shows that the existence of these smooth pathways for cracking encourages shear sliding.

Fig.1.Schematic illustration of crack opening forces generated through different mechanisms at the micro-scale under deviatoric stress condition.High(Hi)and low (Lo)E and v are the Young’s modulus and Poisson’s ratio,respectively(after Diederichs (2000)).

Fig.2.Micro-mechanisms for compression-induced tensile stress at the contacts between the particles in DEM.(a)The compressive force can be resolved into tensile and shear stresses acting on the boundaries between the grains.(b)The compressive forces are resolved into tensile stress at the contacts and moments trying to rotate the disks.

Fig.3.Smooth versus rough sliding paths for shearing in the(a)trigon and(b)Voronoi models,respectively.The potential sliding paths are shown as red lines.

3.Microstructure modelling for DEM

3.1.3D Voronoi tessellations

A 3D Voronoi tessellation is a partition of a domain of 3D space, D ?3,into a collection of cells.Given a number of seed points in D, {Si(xi)}for i={1,…,N},every seed is assigned a Voronoi cell,Ci,as follows:

where d(?,?)is the Euclidean distance.The seed locations are randomly chosen from a uniform distribution in the domain.A Voronoi cell can be seen as the region of in f l uence of a seed,namely the region of space closer to the seed than to any other seed.Voronoi cells are convex polyhedra intersecting along f l at faces, straight edges and vertices,for 2,3 and 4 cells,respectively.An example of a Voronoi tessellation containing 2000 cells is provided in Fig.4a.Different algorithms have been proposed for constructing Voronoi tessellations.In the present work,a cell-by-cell construction algorithm is used.A cell(Ci,of seed Si)is f i rst set as the whole domain.Ci is then modi f i ed by an iterative process,where other seeds(Sj)are considered by increasing distance from Si.At each iteration(seed Sj),Ciis reduced to the intersection of the previously computed cell and the half-space closer to Sithan to Sj.The iterative process can be stopped when the distance between seed Siand seed Sjbecomes high enough for the half-space closer to Sithan to Sjto necessarily include the whole cell.A simple,isotropic criterion is

where{Vi}is the set of vertices of cell Ci.The resulting collection of cells fully de f i nes the Voronoi tessellation.However,cell interfaces are de f i ned several times:a cell face,edge or vertex is de f i ned in 2, 3 or 4 cells,respectively.Therefore,for a proper description of a tessellation,duplicate features are merged.In general,this is also needed for further processing,for example regularization(see Section 3.2)and meshing(see Quey et al.(2011)).

The interest of Voronoi tessellations for microstructure modelling is twofold.First,similarly to regular 3D tessellations,based on cubes,rhombic dodecahedra or truncated octahedra,Voronoi tessellations can be described in a vectorial format(contrary to a raster format as provided by tomography techniques or different modelling approaches).Such a compact description is particularly wellsuited for DEM and can be meshed for other methods.Second, Voronoi tessellations show a wide variety of cell shapes and volumes,and spatial orientations of the contact surfaces between cells.This guarantees that the microstructure behaviour is not biased by particular cell arrangements or interface spatial orientations.This is why 3D Voronoi tessellations have been often used for f i nite element method(FEM)or boundary element method (BEM)simulations,e.g.Barbe et al.(2009),Barbe and Quey(2011), Quey et al.(2012),Benedetti and Aliabadi(2013).

Fig.4.Effect of regularization illustrated on a 2000-cell tessellation.(a)Standard Voronoi tessellation and(b)regularized Voronoi tessellation.Only half of the cells are shown and the small edges are plotted in red.8000 out of 25,000 edges are removed by regularization.

3.2.Regularization

As can be seen in Fig.4a,Voronoi tessellations contain a high proportion of edges and faces of dimension signi f i cantly lower than the average cell size.Such features are negligible from a cell morphology point of view and therefore should not play a signi f icant role in the results of the DEM simulations either.They also are highly detrimental to computation time in DEM simulations.

These drawbacks can be circumvented by“regularization”,a technique consisting of removing the“small edges”of a tessellation,of length typically lower than a tenth of the cell size(Quey et al.,2011).Regularization removes small faces indirectly, because the latter are composed of one or several small edges. During regularization,the small edges are removed by increasing length.An edge(and its two vertices)is collapsed to a single vertex on the condition that the resulting geometrical distortion in the vicinity of the edge is suf f i ciently low(typically,the neighbouring faces must remain planar within a 200 angular tolerance).Constraints are imposed on the vertices,edges and faces at the boundary of the domain so that the domain shape is not changed. As illustrated in Fig.4b,regularization removes the small edges and faces while retaining the global cell shapes.In most cases,face distortions remain lower than 1?-2?(Quey et al.,2011).In 3DEC simulations,such distorted faces can be treated as standard,planar faces.

Fig.5.Two-scale Voronoi tessellations applied to modelling(a)a two-scale material (e.g.sedimentary rocks)and(b)intra-grain cracking paths.

3.3.Two-scale Voronoi tessellations

A“two-scale Voronoi tessellation”is a tessellation for which every cell of a f i rst,“primary”Voronoi tessellation,is partitioned by another,“secondary”Voronoi tessellation.Fig.5 provides two possible applications of two-scale Voronoi tessellations.First,such tessellations enable to model microstructures that clearly exhibit two scales such as sedimentary rocks,consisting of 45?-tilted bands decomposed into individual blocks(Fig.5a).The primary Voronoi tessellation was created from seeds regularly spaced along a line (hence the band structure)and the secondary Voronoi tessellations were obtained using random seed locations.Second,two-scale tessellations enable inclusion of potential intra-grain crack paths into a regular,single-scale material(Fig.5b).Intra-grain interfaces can be assigned micro-properties different from grain-to-grain interfaces.The algorithm for constructing a two-scale Voronoi tessellation consists of sequentially constructing standard,singlescale Voronoi tessellations.Secondary Voronoi tessellations can be constructed independently from each other.The Voronoi tessellation algorithm described in Section 4.1 is used,using the primary cell(instead of the whole domain)as a starting point for the secondary cells.This results in a collection of secondary tessellations.Similarly to single-scale Voronoi tessellations,duplicated,or overlapping features are found at the intersections between secondary tessellations.These features must be merged to get a proper description of the tessellation.The procedure involves merging,creation and decomposition of vertices,edges and faces located at the intersection between secondary tessellations.Details will be provided separately.Similar to single-scale tessellations,the resulting two-scale tessellations contain small edges and faces, which can be removed by regularization.

3.4.The Neper software

All algorithms described previously have been implemented in a free(open-source)software:Neper(Quey et al.,2011;Quey,2014). The software offers support for both single-scale and two-scale Voronoi tessellations,in 2D and 3D,however the focus of this paper is on single-scale 3D Voronoi tessellated models.Thanks to ef f i cient and robust algorithms,very large tessellations can be easily created(typically 105 cells).Any convex domain shape can be used;cubes(or any parallelepiped)and cylinders are directly supported.The tessellations can also be stretched to get elongated cells.Rotation and cropping are also straightforward operations. Finally,tessellations can be written at various formats,including 3DEC for DEM simulations,and meshed for FEM simulations.

4.Grain-based model for 3DEC

The Voronoi blocks can be used for simulating fractures at different scales,from micro-cracks within a laboratory specimen (Lan et al.,2010)to modelling spalling and slabbing in the walls of underground excavations(Itasca,2011).The focus of this paper is on GBMs.With the currently available computational power,GBMs are commonly restricted to the laboratory specimen scale.To investigate the application of the 3D Voronoi tessellation approach for simulation of brittle micro-fracturing in rock,a series of Uniaxial Compressive Strength(UCS)tests are simulated in 3DEC from the tessellations generated in Neper.

4.1.Specimen generation

Cylindrical specimens were decided to be used for the simulation of numerical UCS tests.Any 3D convex domain geometry can be f i lled with 3D Voronoi diagrams with Neper.To avoid grains whose geometry is affected by the domain boundary in a cylindrical specimen,a cubical domain was generated which was then carved into a cylindrical specimen in 3DEC(Fig.6).For generation of the specimens,a cylinder with the height of 140 mm was cut from the centre of the tessellated block and then two 10 mm slices were cut from each end of the specimen.The grains in these 10 mm thick slices were joined to form the platens.The f i nal specimen had a diameter of 55 mm and a length of 120 mm.The cubical domains in Neper were generated with 40,000 grains while the f i nal cylindrical specimens contained 5300 grains on average.

4.2.Numerical test set-up

Fig.6.Generation of a Voronoi tessellated cube in Neper and cutting the cube into a cylindrical specimen in 3DEC.

Loading of the UCS models was performed by applying a constant velocity to the top and bottom platens until a moderate stateof monotonically increasing axial stress was reached in the specimens.After this stage,a servo-controlled system took over the loading control.This servo mechanism kept the average unbalanced mechanical forces in the model within a speci f i ed range by adjusting the applied velocity to the platens.The unbalanced force servo system is particularly important at the failure stage to keep the loading condition within the static range.The loading velocity at the beginning of the test and the speci f i ed range for unbalanced forces for the servo-mechanically controlled loading stage was kept small enough to ensure the uniform distribution of stress within the specimen.

The axial strain for the specimens was measured by averaging the displacements between f i ve pairs of grid points,located 10 mm away from each end of the specimens in the vicinity of the models’axial axis,in principle similar to laboratory extensometers.The locations of these grid points are illustrated in Fig.7.The specimens’lateral strain was calculated from the average relative displacement of three pairs of grid points in the middle of the specimens as shown in Fig.7,in the x and y directions.

The axial stress in the specimens was measured by averaging the axial stress in the zones of blocks whose centroid falls within a cube in the centre of the cylindrical specimen.A trial test was run to determine the suitable size for the cube which was representative of the axial stress within the specimen.For this purpose,the axial stress within four boxes with dimensions shown in Fig.8a was measured during a UCS test.The results are shown in Fig.8b.It can be seen that the difference between the measured stresses tends to disappear with increasing box size,consistent with the emergence of a representative volume element(Kanit et al.,2003).Although the stress measured in the two largest cubes appear to be very similar,it was decided to use the biggest box size with dimensions of 35 mm?35 mm?40 mm to measure axial stress within the specimens.An example for the zones used for measurement of axial stress in a specimen is shown in Fig.8c.

4.3.Estimation of crack damage thresholds

Estimation of the crack damage thresholds for the models was performed according to the Acoustic Emission(AE)method described in Ghazvinian et al.(2012).The comparable phenomenon to AE activity in laboratory specimens can be easily correlated to the bond breakage in PFC models.However,the same is not true for GBMs in 3DEC.Contact deletion,sub-contact deletion,sub-contact failure,area of the failed sub-contact,etc.,can all be associated with the AE activity.In this paper,the number of failed sub-contacts was used,as shown in Fig.9,to identify the damage thresholds. Nevertheless further investigation is required on this subject.

4.4.GBM calibration

The behaviours of 3D Voronoi tessellated models are controlled by two main components:the grains(3DEC blocks)and the contacts between grains.The grains in a 3DEC model can be simulated as rigid bodies,either elastic or plastic,where continuum intragrain failure is allowed to occur.Common practice with Voronoi GBMs is to have rigid(Kazerani and Zhao,2010)or elastic grains (Itasca,2011)for simplicity and to allow for failure to occur only at contacts between the grains.In this case,the parameters requiring calibration,when a simple Coulomb slip joint model is used for the grain contacts in 3DEC,are as follows: ?

contact normal stiffness(kn) ?

contact normal to shear stiffness ratio(kn/ks) ?

contact peak cohesion(cp) ?

contact residual cohesion(cr) ?

contact peak tensile strength(Tp) ?

contact residual tensile strength(Tr) ?

contact peak friction angle(φp) ?

contact residual friction angle(φr) ?

grain(block)Young’s modulus(E) ?

grain(block)Poisson’s ratio(ν)

For simplicity in the calibration process,the residual tensile strength(Tr)and cohesion(cr)are commonly assumed to be zero (Kazerani and Zhao,2010;Itasca,2011;Gao and Stead,2014)for Voronoi GBMs.

Martin(1997)showed that the mobilization of the frictional and cohesional components of shearing in brittle rocks is not a simultaneous process.The friction gradually mobilizes as a function of damage to the brittle rock,while the cohesion is degrading. Hajiabdolmajid et al.(2002)successfully implemented the Cohesional Weakening-Frictional Strengthening(CWFS)behaviour in their continuum modelling to simulate notch formation around the tunnels of deep underground openings in massive brittle rocks.To adopt the CWFS concept for the contact behaviour between the grains,the peak cohesion(cp)and residual friction angle(φr)were de f i ned for the contacts while the residual cohesion(cr)and peak friction angle(φp)were set to zero.Therefore for simplicity from now on in this paper,the peak cohesion,peak tensile strength and residual friction angle will be referred to as c,T andφ,respectively. The schematic representation of the contact behaviour with the implemented CWFS approach in the generated GBMs is shown in Fig.10.

For the sensitivity analysis of the GBM response to the grain contact properties,a set of micro-parameters were chosen to closely replicate the mechanical properties of the Lac du Bonnet granite(Diederichs,2000).The micro-properties were adjusted accordingly in an iterative process suggested by Potyondy and Cundall(2004)and Itasca(2008)until the desired macroresponse was reached in terms of UCS and other macroproperties.The f i nal micro-properties,model properties and stress-strain curve for the model are shown in Table 1 and Fig.11, respectively.The micro-properties listed in Table 1 will be used in Section 5 for the parametric study.

Fig.7.Location of grid points used to measure axial and lateral strains.

Fig.8.Measurement of axial stress for the specimens.(a)The four box sizes investigated for axial stress measurement.(b)Comparison between the measured axial stress within the four boxes of different dimensions.(c)Axial stress for the specimen is measured as the mean axial stress recorded in the red zones.

5.Sensitivity analysis for GBM micro-properties

The primary requirement prior to any GBM modelling is to calibrate the model micro-properties.The calibration process can be challenging and laborious if done in a random fashion.Therefore,a series of sensitivity analyses were performed on the GBM contact properties and grain size to better understand the role of these micro-properties on the overall response of the model and to establish the calibration procedure for a 3D Voronoi GBM.The sensitivity of the model response to the contact micro-properties was investigated by varying one micro-property at a time while the rest of the micro-parameters were kept constant,equal to the values listed in Table 1.The grain elastic constants(Young’s modulus and Poisson’s ratio)were considered real and were obtained from reported laboratory tests on the composing minerals of the Lac du Bonnet granite(Lama and Vutukuri,1978)as opposed to rigid(Kazerani and Zhao,2010),or semi-rigid(by assuming very high stiffnesses for them(Itasca,2011)).

5.1.Contact normal stiffness(kn)

The earlier work by Kazerani and Zhao(2010)on 2D rigid Voronoi GBMs and similar work by Gao and Stead(2014)on UDEC (Itasca,2013b)trigon models,suggested that the simulated Young’s modulus of the specimen is strongly correlated with the normal and shear stiffnesses of the contacts between the grains.They showed that by increasing the normal and shear stiffnesses of the contacts,the model’s Young’s modulus increases.Similarly,the contact stiffness in the 3D Voronoi models directly controls the specimen’s Young’s modulus as shown in Fig.12.The increase in the normal and shear stiffnesses of the contacts stiffens the model’s contact lattice and consequently enhances the stiffness of the model.Stiffening the models contact lattice also increases the role of grain deformation in controlling the model’s macro-behaviour. Therefore,increasing the contact normal and shear stiffnessvalues changes the models Poisson’s ratio to within a closer range to the grain’s Poisson’s ratio.For the models shown in Fig.12,the models overall Poisson’s ratio approaches 0.20(grains Poisson’s ratio)with increasing contact stiffness values.The Poisson’s ratio of the model with constant kn/ksratio is independent of the change in the contact stiffness when grains are assumed rigid(Kazerani and Zhao,2010).

Fig.9.Estimation of crack damage thresholds from the number of failed sub-contacts in 3DEC.

Fig.10.3DEC simulation of fracturing along grain boundaries.(a)Sequential mobilization of cohesion and frictional component of contact shear strength.(b)Constitutive behaviour of grain contacts with implemented CWFS approach(σnTis the normal stress at which the contact properties are switched from peak to residual,σsmaxis the maximum shear stress at a contact with normal stress greater thanσnT).

Table 1Micro-parameters used for calibration of the 3D Voronoi GBM and the resulting model macro-response.

Fig.11.Stress-strain curve for the UCS model with the micro-parameters listed in Table 1.

5.2.Contact normal to shear stiffness ratio(kn/ks)

The Poisson’s ratio is directly governed by the normal to shear stiffness ratio of the grain bonds in the BPM formulation as demonstrated by Diederichs(2000),Potyondy and Cundall(2004) and Cho et al.(2007).Similar works in UDEC with rigid Voronoi grains(Kazerani and Zhao,2010)and elastic trigon grains(Gao and Stead,2014)support this.

The stress-strain curves from the sensitivity analysis for the contact kn/ksratio for the 3D Voronoi GBM is shown in Fig.13.In this series of models,the contact normal stiffness(kn)was kept constant,at 68,000 GPa/m,while different kn/ksratios(2,5,8 and 16)were examined.The change in the model’s elastic constants as a function of kn/ksratio is plotted in Fig.14.This f i gure shows that increasing the kn/ksratio increases the Poisson’s ratio,even beyond the realistic threshold of 0.5,before the crack initiation threshold is reached.Conversely,increasing the kn/ksratio decreases the Young’s modulus and reduces the brittleness(increases ductility)of the model.By decreasing the shear stiffness and keeping the normal stiffness constant(increasing kn/ksratio),while this makes the modelsofter and consequently lowers the Young’s modulus,the contacts are encouraged to dilate in shear as oppose to extensile opening,which results in a larger lateral deformation and thus a larger Poisson’s ratio.Therefore in 3D Voronoi GBMs,the kn/ksratio is not only the major controlling factor for the Poisson’s ratio,but also controls the Young’s modulus and the brittle to ductile transition of the model(Fig.13).

Fig.12.Effect of contact normal stiffness on the model’s elastic constants.

5.3.Contact c/T

The contact or bond strength properties in UDEC and PFC directly control the strength of the models in compression or tension.Increasing the contact strength-related micro-properties improves the specimen’s strength(Diederichs,2000;Potyondy and Cundall,2004;Gao and Stead,2014).In Diederichs(2000),the shear to normal strength ratio of the contact bonds in PFC controls the peak strength and brittleness of the model.The material response shifts towards a more brittle fashion with increasing shear to tensile bond strength ratio.This is because increasing this ratio for the contact bonds forces the micro-fractures to form mostly in tension,which is similar to the initiation of extensile fractures in laboratory brittle rock specimens(Diederichs,2007).

A suite of models were run with identical grain contact peak cohesion(c),but varying peak tensile strengths(T)for the grain contacts to achieve different ratios of cohesion to tensile strengths (c/T)for the 3D Voronoi GBMs.The UCS and UCS/CI results for this series of models are shown in Fig.15.As expected,the peak strength of the specimen increases with decreasing c/T ratio,since T is increasing,causing the strength of the grain contacts to improve.The other macro-property that is affected by the c/T ratio is the CI threshold in the models.The opening of fractures at the CI threshold is assumed to be in an extensile fashion and therefore governed by the tensile strength of the contacts(Diederichs,2003). This is supported by the models in Fig.15 which shows that by decreasing the c/T ratio,the UCS/CI ratio decreases drastically (increasing CI with increasing contact tensile strength).

Fig.13.Stress-strain curves for the models with constant normal stiffness (68,000 GPa/m)but different ratios of contact normal to shear stiffness.

Diederichs(2000,2003)reported limited success with PFC2Dfor simulation of crack propagation beyond a single disk at the CD threshold.This is due to stress redistribution amongst the bonds in the vicinity of a bond once it breaks,which impedes stress concentration at the crack tips and therefore suppresses crack extension.Gao and Stead(2014)argued that they have successfully simulated the crack propagation and coalescence within their trigon approach in UDEC and 3DEC.This was however solely based on the visual evaluation of contact openings.Martin and Chandler (1994)identi f i ed the deviatoric stress corresponding to the maximum volumetric strain of a laboratory specimen under compressive loading,as the onset of unstable cracking,or in other words,the CD threshold.Here,the volumetric strain is plotted versus the axial stress normalized to the peak strength of each specimen for models with different c/T ratios in Fig.16.If the reversal of volumetric strain can also be identi f i ed as the crack propagation threshold in the numerical models,it can be seen that unstable micro-cracking or crack propagation occurs in models with c/T ratios equal to 4,10 and 100.In Fig.16,the decreasing c/T ratio shifts the axial stress associated with the reversal point of the volumetric strain closer to the UCS.For c/T=2,the volumetric strain of the model never reaches the reversal point.Trans-granular cracks are major contributors to unstable cracking and the crack propagation threshold in brittle rocks.The absence of this factor can suppress the unstable cracking to a large extent in GBMs with unbreakable grains.Sub-tessellation of the Voronoi grains as described in Section 3.3 allows for trans-granular fracture propagation and can be a solution to this limitation.The decreasing c/T ratio as observed in Fig.16 increases the slope of the volumetric strain curve for the models.This can be attributed to the increase in the tensile strength of the contacts which allows for larger contact dilation before failure and therefore larger volumetric strains at the failure stress.

The cohesion,tensile strength and friction angle of the grain contacts are also shown in 2D Voronoi simulations of triaxial tests by Kazerani and Zhao(2010)and Gao and Stead(2014),to control the macro cohesion,tensile strength and friction angle of the model respectively.

Fig.14.Effect of contact normal to shear stiffness ratio(kn/ks)on the model’s elastic constants.

Fig.15.Effect of contact c/T ratio on the model’s peak strength and UCS/CI ratio.

Fig.16.Effect of contact c/T ratio on the model’s volumetric strain behaviour.

5.4.Contact friction angle(φ)

To more accurately simulate the behaviour of brittle rocks,the CWFS approach was employed in the models and was discussed in detail in Section 4.4.The peak friction angle of the grain contacts was set to zero,therefore the material friction is controlled by the residual friction angle.This approach is more common in PFC modelling(Diederichs,2003).According to Kazerani and Zhao (2010)and Gao and Stead(2014),the material friction angle can be controlled by the peak friction angle.In the CWFS approach,the material friction angle is governed by the residual friction angle. The results of a series of 3D Voronoi GBMs with different residual friction angles showed that,as expected,increasing the residual friction angle improves the strength of the models.

5.5.Effect of grain size

A comprehensive investigation on the effect of particle size on the PFC2Dand PFC3Dmodel macro-properties was done by Potyondy and Cundall(2004).They showed that in the PFC2Dmodels,the elastic constants and the peak strengths are independent of the particle size while in PFC3D,the peak strength and the Young’s modulus of the model are governed by the particle size. Similar to PFC2D,the Poisson’s ratios of PFC3Dmodels are independent of the particle size.

To investigate the effect of GBM grain size on the model response,a set of cylindrical specimens were carved from three cubes f i lled with 20,000,40,000 and 60,000 Voronoi grains.The resulting specimens contained approximately 3000,5300 and 7800 grains,with average grain volumes of 0.111 cm3,0.062 cm3and0.042 cm3,respectively.The stress-strain curves obtained from the UCS simulation for these models are shown in Fig.17.The results show that the peak strength of the models is independent of the grain size,within the grain sizes investigated in these models.A minor dependency of the model’s elastic response to the grain size is observed to occur after the initiation of damage in the models (deviation of lateral strain from linearity(Martin,1993)).As fracturing initiates,the macro-response of the models becomes slightly softer with decreasing grain size.This can be attributed to the relative magnitude of the grains elastic constants in comparison to the shear and normal stiffnesses of the grain contacts.

6.Calibration procedure for GBMs with regular 3D Voronoi grains

A better understanding of the role of each of the microparameter on the overall mechanical response of the GBMs was achieved from the sensitivity analyses.Therefore,the following procedure can be suggested for calibration of 3D Voronoi grainbased models from UCS data:

(1)Decide on the Voronoi grain size(number of Voronoi cells within the specimen)based on the rock mineralogy and available computational power.

(2)De f i ne the elastic constants(Young’s modulus and Poisson’s ratio)for the grains according to laboratory testing or available literature for the composing minerals of the simulated rock.

(3)Set the grain contact strength related components to high values or run the 3DEC models in“small strain”mode for calibration of the model elastic constants.

(4)Adjust the contact normal to shear stiffness ratio(kn/ks)for the correct Poisson’s ratio.

(5)Adjust the normal stiffness(kn)for the model to achieve a proper Young’s modulus.

(6)Cohesion,c/T and the residual friction angle of the grain contacts are adjusted in the model to achieve the correct peak strength and crack initiation threshold.In the case of matching the model response to triaxial and tensile strength testing data, the model cohesion,friction angle and tensile strength can be controlled by the grain contact cohesion,residual friction angle and tensile strength,respectively.

(7)A distribution of strengths can be assigned to the grain contacts to further re f i ne the crack initiation threshold.

(8)Steps 3-6(without changing the contact strength parameters to high values)can be reiterated for a more accurate calibration.

The suggested calibration procedure is for GBMs with elastic grains and grain contact behaviour as shown in Fig.10.

7.DEM modelling of damage progression in anisotropic rocks

Many rocks near the Earth’s surface represent anisotropic behaviour due to layering,foliation,f i ssuring,bedding,strati f i cation,jointing,etc.(Amadei,1996).The available fabric arrangement in rocks directly in f l uences their anisotropic behaviour(Gatelier et al.,2002).The term“fabric”is used in this paper as a general term to describe any planar feature in intact rock.The available fabric in intact rocks can be represented by one or a combination of the fabric elements shown in Fig.18.The difference in mineral composition(Fig.18a),preferred orientation of elongated grains (Fig.18b,c),spatial variation of grain sizes(Fig.18d),preferred orientation of platy minerals,aggregates or planar micro-fractures (Fig.18e-g)or the existence of any combination of these elements in intact rock(Fig.18h)can introduce anisotropic effects to its mechanical behaviour.

Previous attempts for explicit simulation of mechanical anisotropic behaviour for brittle intact rocks by different researchers have mainly concentrated on introducing anisotropic damage to the intact rock,similar to the fabric element as shown in Fig.18g. Diederichs et al.(2004)studied the effect of pre-existing damage in brittle rocks on their mechanical behaviour such as CI and CD by inducing different crack intensities(removing or breaking bondsbetween disks)at different orientations in their PFC2Dmodels.In a hybrid continuum/discontinuum approach,Lisjak et al.(2014a) successfully applied a small scale Discrete Fracture Network (DFN)to a combined FDEM model to study the anisotropic behaviour of Opalinus Clay.They later extended this approach to a larger scale to study the development of the Excavation Damaged Zone(EDZ)around tunnels excavated in this formation(Lisjak et al., 2014b).

Fig.17.Effect of grain size on the material response.

The Voronoi tessellation provides an opportunity to generate GBMs that resemble the fabric elements similar to the ones in Fig.18c-f at the micro-scale.In the simplest form,a GBM with different scaling factors can be generated from a regular Voronoi tessellation to mimic the intensity of the anisotropy in the intact rock.Fig.19a-c which host the fabric elements as shown in Fig.18c can be approximated by scaled Voronoi grains with different length/width ratios.

Fig.18.Fabric elements that can introduce anisotropic behaviour in intact rocks(after Passchier and Trouw(2005),modi f i ed from Hobbs et al.(1976)).

7.1.Laminated GBM for fabric-guided micro-fracturing simulation

Laminated GBMs with fabrics composed of elongated 3D Voronoi grains can be constructed with scaling.An example for generation of a cylindrical laminated specimen oriented 45?with respect to the loading direction is illustrated in Fig.20.A prismatic domain with height/width ratio equal to the desired length/width ratio of grains is generated,scaled and then rotated for proper fabric orientation.The generated model is then exported to 3DEC and cut for a cylindrical specimen.

7.2.Calibration

To investigate the application of the laminated GBMs for simulation of anisotropic rock behaviour,a series of UCS models were generated in 3DEC by using the elongated grains for different fabric orientations,following the procedure introduced in Section 7.1.The Cobourg Limestone,which is composed of alternating layers of fossil rich packstone and argillaceous wisps and blebs,was chosen for the model veri f i cation.The mechanical properties of the Cobourg Limestone were obtained from Ghazvinian et al.(2013a, b).In these studies,Cobourg Limestone specimens with different layering orientations with respect to the loading direction were drilled from large diameter cores obtained from a limestone quarry in Bowmanville,Ontario.For f i ve groups of specimens with apparent bedding oriented at 0?,30?,45?,60?and 90?,four to six specimens were tested for their mechanical properties includingthe peak strength,elastic constants and crack damage thresholds. An example of the Cobourg Limestone specimens and the associated GBMs are shown in Fig.21.

Fig.19.Crystalline minerals with various degrees of recrystallization.(a)Polygonal fabric of scapolite(width of view 4 mm).(b)Elongated grains in recrystallized quartz (width of view 1.8 mm).(c)Strongly orientated mineral grains de f i ned by parallel grains of biotite,muscovite and quartz(width of view 1.8 mm),(modi f i ed after Passchier and Trouw(2005)).

In addition to the micro-parameters discussed in Section 4.4 that need to be calibrated for GBMs with regular Voronoi tessellations,here the length to width(l/w)ratio of the grains for laminated GBMs also needs to be calibrated.For a speci f i c rock,the models with different fabric orientations have identical microparameters.It is the l/w ratio of the grains which de f i nes the mechanical anisotropy for the behaviour of the GBMs with different fabric orientations.

The calibration process was conducted by deciding on an initial value for the grains l/w ratio and adjusting the micro-parameters for the horizontally(90?)and vertically(0?)laminated models simultaneously with an identical set of micro-parameters to achieve the proper elastic constants.Then,the strength-related microparameters were adjusted for the horizontal model to calibrate its strength properties.Subsequently by applying those microparameters to the vertically laminated model and f i ne-tuning the grains l/w ratio to attain the proper strength macro-properties,the global set of micro-parameters for the entire suite of specimens with different fabric orientations was obtained.The stress-strain curves for the numerical and laboratory specimens with vertical and horizontal fabrics are shown in Fig.22.The calibrated microproperties,listed in Table 2,were applied to the suite of GBMs with f i ve different fabric orientations as shown in Fig.21.

Fig.20.Generation of layered GBMs for anisotropic simulations in 3DEC(an example of a 45?specimen).(a)Generation of 3D Voronoi tessellations with a prismatic domain.(b) Scaling the domain to achieve the layered model with elongated grains.(c)Rotation of the domain to achieve the desired layering orientation.(d,e)Cutting through the block to make the cylindrical specimen and platens.

7.3.Numerical results

The laminated GBM approach was successful in capturing the effect of fabric orientation on the failure mode of the models.In Fig.23,which shows the fractures in the models after failure on a plane cut through the centre of the model,the switch between the failure modes can be observed to occur as a function of the fabric orientation.Axial splitting and brittle failure can be seen in the 0?and 90?specimens respectively and the formation of shear bands is evident in the other orientations.Remarkably the formation of axial micro-fractures and their role towards formation of a shear band by kinking can be seen in the 45?model.

The comparison between the laboratory testing results of the Cobourg Limestone and the corresponding numerical macroproperties are shown in Fig.24.The UCS and CI thresholds between the two sets of data are in good agreement and the classical U-shape behaviour(Jaeger and Cook,1969)for anisotropic rocks is well-captured by the GBMs.The elastic response of the numerical models,which were calibrated before the strength properties, approximately captures the lower-bound Young’s moduli(except the vertically and horizontally laminated models)and the upperbound of the laboratory measured Poisson’s ratios.Since calibration of the strength-related components of the macro-behaviour for GBMs takes place after the calibration of the elastic response, the Young’s moduli and the Poisson’s ratio of the models can be slightly altered,as observed in Fig.24.Therefore,following the calibration procedure for a couple of iterations between adjustments for the elastic response and strength of the models can increase the accuracy of the calibrated micro-parameters.

7.4.Sensitivity analysis

The micro-properties that can be calibrated for a speci f i c set of macro-properties in general are non-unique.A series of sensitivity analyses were conducted for the most important micro-propertiesof a laminated GBM to better understand their role in de f i ning the classical U-shape behaviour that is anticipated from anisotropic rocks.The U-shape behaviour was f i rst introduced by Jaeger and Cook(1969)to show the effect of discontinuities on the rock mass strength.It is employed here to demonstrate the variation of mechanical properties with respect to the intact rocks fabric orientation.

Fig.21.The Cobourg Limestone specimens with different fabric orientations with respect to the loading direction(before testing)and the corresponding 3DEC instances.

Fig.22.Stress-strain curves for 0?and 90?specimens tested in the laboratory and numerically simulated in 3DEC.

7.4.1.Scaling(grain l/w ratio)

Numerous models with different grain l/w ratios were achieved by scaling prismatic domains of different heights,but all f i lled withan identical number of grains,down to a proper cube.This was discussed in detail in Section 7.1.Two new sets of laminated GBMs with l/w ratios of 1.5 and 2.5 were generated and numerically tested by using the micro-properties listed in Table 2.The results for these two new sets of models are compared in Fig.25 with the original suite of GBMs as discussed in Section 7.3.

Table 2Calibrated micro-parameters based on the horizontal and verticallaminated models.

The l/w ratio of grains appears to have a major impact on the strength of the horizontally laminated models,while its effect on the strength of models with other fabric orientations seems to be minute.Regarding the Young’s modulus and Poisson’s ratio,the l/w ratio appears to control the depth of the anisotropy U-shape curve, and the elevation of the curve.

Fig.23.Cross-sections through the laminated GBMs after failure.

7.4.2.Contact c/T ratio

The c/T ratio of contacts in GBM plays a key role in controlling the macro-response of the models as observed in Section 5.3.Two new sets of models with identical micro-parameters to the original set of laminated models(as listed in Table 2,c/T=10),but different contact tensile strengths(T),to achieve c/T ratios of 2 and 20 were numerically tested.The results are plotted in Fig.26.Similar to the GBM with regular Voronoi grains,reducing the c/T ratio by increasing the tensile strength(T)of contacts improves the peak strength of the models and subsequently shifts the U-shape anisotropy curve to a higher elevation(model with higher peak strengths).This is also true for the effect of the c/T ratio on theYoung’s moduli of the models.The c/T ratio controls the f l atness and the direction of the anisotropy curve for the CI thresholds and Poisson’s ratios as shown in Fig.26.The c/T ratio is particularly important for its role in controlling the anisotropy curve for the CI threshold.In some rocks,the CI threshold is entirely independent of the fabric orientation(Hakala et al.,2007)while it can change according to the orientation of the lamination in other rock types (Ghazvinian et al.,2013a).Therefore the c/T parameter can be used to adjust the sensitivity of the CI threshold to the orientation of the fabric present in a rock.

Fig.24.Comparison between the mechanical properties of the Cobourg Limestone measured in the laboratory and simulated in 3DEC(the dotted lines are 2nd order polynomial best f i t to the 3DEC results).

Fig.25.Effect of grain l/w ratios on the GBM macro-response(the dotted lines are 2nd order polynomial best f i t to the data points).

7.4.3.Contact friction angle

Increasing the friction angle of the contacts(i.e.the residual friction angle according to the contact model as shown in Fig.10) improves the peak strengths of the models.The increase in the strength of the models intensi f i es as the orientation of fabric in the rock changes from vertical to horizontal.As shown in Fig.27,the maximum strength increase of the laminated GBMs as a function of increasing the models’friction angle is observed in the horizontally laminated model.

7.5.Calibration procedure

Veri f i cation of the laminated GBM approach with the available test data from the Cobourg Limestone con f i rmed the application of this modelling method for numerical replication of brittle rock behaviour with inherent anisotropic mechanical properties.Based on the conducted sensitivity analyses for micro-properties,the calibration procedure for the GBMs with elongated 3D Voronoi grains can be summarised as follows:

(1)Decide on the grain size for the GBM and an initial grain l/w ratio.

(2)De f i ne the elastic constants(Young’s modulus and Poisson’s ratio)for the grains.It is nearly impossible to match the elastic constants for the horizontally and vertically laminated models during the calibration stage with an identical set of micro-parameters with rigid grains or unrealistically very stiff grains.

(3)Set the strength components of the contacts to a high value so that the deformational behaviour of the model is fully elastic (or set the“small strain”in 3DEC)and adjust the Poisson’s ratio for the models with the horizontal and vertical fabric simultaneously by modifying the kn/ksratio.

(4)Calibrate the Young’s moduli for the models with the horizontal and vertical fabric simultaneously by adjusting the knfor the contacts.

(5)Calibrate the peak strength and CI threshold for the horizontally laminated model by f i ne-tuning the c/T ratio and c.Adjust the friction angle based on the post-peak behaviour(if needed).

(6)Apply the micro-parameters calibrated thus far to the vertically laminated model and adjust the grains l/w ratio accordingly.

(7)Apply the micro-parameters calibrated thus far to the entire suite of laminated GBMs and correct the U-shape anisotropy curve by adjusting the micro-parameters according to the discussions in Section 7.4 for a few iterations.

8.Discussion

This paper showed that brittle rock damage through explicit simulation of micro-fracturing can be accurately captured in numerical models by means of GBMs built with 3D Voronoi grains. Similar work has been done previously in 2D(e.g.Kazerani and Zhao,2010;Lan et al.,2010;Itasca,2011),however the complexity of generating proper 3D Voronoi tessellations has encumbered its development in 3D DEM codes to this point.The mechanics of the Voronoi tessellated models at the micro-scale remains nearly the same in transition from 2D to 3D.However some complexities in regards to the calibration process can intensify when switching from UDEC to 3DEC for Voronoi GBMs.In UDEC and 3DEC,the contact models commonly used for grains are the constitutive slip models developed for rock mass joints and therefore,attempting to formulate the micro-scale grain contact behaviour by using an analogue established for large-scale discontinuities can be challenging.In this study,the CWFS approach was successfully implemented in the grain contact model to better address the peak and residual friction angles and gradual mobilization of the frictional strength component of the contacts.This approach can be further re f i ned with GBM simulations under triaxial conditions to capture the correlation between the material friction angle and the friction angle in the CWFS model for grain contacts.

Grain contact normal and shear stiffnesses are two of the most important micro-parameters in terms of de f i ning the model’s macro-response.The transition from 2D to 3D models drastically increases the importance of these two micro-parameters.The unrealistic representation of contact stiffnesses(particularly the normal stiffness)with a constant value leads to high stiffness parameters for the contacts and therefore an enormous stiffness contrast between the grain and the grain contact,which can be erroneous.It is anticipated that de f i ning the normal stiffness of thegrain contact as a function of the grains’overlapping area(in 2D)or volume(in 3D)can better represent the micro-mechanical behaviour of crystalline rocks.

3D Voronoi tessellated models can be used as the building block for the Synthetic Rock Mass(SRM)(Mas Ivars et al.,2011)to represent intact rock.Its combination with the in-situ joint network,commonly referred to as the DFN,can represent rock masses which can be numerically investigated.The regular and laminated GBMs can be used for SRM representation of rock masses,including intact rock with isotropic or inherent anisotropic mechanical behaviour.Success of the 3D Voronoi GBMs in simulating the micro-fracturing in laboratory scale specimens supports its ability in simulation of rock damage and fracturing in the walls of deep underground excavations in brittle rocks by means of proper up-scaling rules.

Fig.26.Effect of contact c/T ratio on the mechanical behaviours of laminated GBMs(the dotted lines are 2nd order polynomial best f i t to the data points).

Fig.27.Effect of contact friction angle on the strength of laminated GBMs(the dotted lines are 2nd order polynomial best f i t to the data points).

9.Conclusions

The concerns with the formulation of bonded particle models (i.e.PFC)lead to increasing application of Voronoi grain-based models for simulation of fracture formation in rocks despite its poor computational speed compared to the BPM.While clumping and clustering in the BPM(Cho et al.,2007;Ghazvinian,2010)and the modi f i cations to the BPM formulation(development of the f l atjoint contact model(Potyondy,2012))in the recent years have solved the concerns with the accuracy of the BPM results to a great extent,the Voronoi based GBMs are still popular due to simplicity in generating these models.To extend the application of the Voronoi tessellated DEM models from 2D to 3D,a new methodology was introduced in this paper for generation of models with any arbitrary 3D convex shape domains in 3DEC tessellated with Voronoi grains.Furthermore,a new model was introduced for simulation of inherent mechanical anisotropic behaviour of rocks containing a fabric based on the freshly developed 3D Voronoi GBM.The fabric-guided micro-fracturing with this approach was veri f i ed successfully with the available laboratory testing data of the Cobourg Limestone.The sensitivity analysis of the grain contact micro-properties proved that the dependency or independency of the CI threshold to the fabric orientation in the rock can be controlled for different rocks and different types of anisotropic behaviours and f i nally was employed to establish a calibration procedure for the 3D Voronoi laminated grain-based models.

Con f l ict of interest

The authors wish to con f i rm that there are no known con f l icts of interest associated with this publication and there has been no signi f i cant f i nancial support for this work that could have in f l uenced its outcome.

Acknowledgements

The authors would like to acknowledge the Nuclear Waste Management Organization of Canada(NWMO),the Swedish Nuclear Fuel and Waste Management Company(SKB)and the National Science and Engineering Research Council of Canada(NSERC)for supporting this research.The discussions with Dr.Jim Hazzard andDr.Branko Damjanac from the Itasca Consulting Group signi f i cantly helped in the preparation of this paper.Special thanks to Felipe Duran,Michelle van der Pouw Kraan and Dr.Derek Martin for their inputs.Assistance received from Mark Jensen and Tom Lam from NWMO,Rolf Christiansson from SKB and Denis Labrie from Natural Resources Canada is much appreciated.

Amadei B.Importance of anisotropy when estimating and measuring in situ stresses in rock.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1996;33(3):293-325.

Barbe F,Quey R,Musienko A,Cailletaud G.Three-dimensional characterization of strain localization bands in high-resolution elastoplastic polycrystals.Mechanics Research Communications 2009;36(7):762-8.

Barbe F,Quey R.A numerical modelling of 3D polycrystal-to-polycrystal diffusive phase transformations involving crystal plasticity.International Journal of Plasticity 2011;27(6):823-40.

Benedetti I,Aliabadi MH.A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials.Computer Methods in Applied Mechanics and Engineering 2013;265:36-62.

Cho N,Martin CD,Sego DC.A clumped particle model for rock.International Journal of Rock Mechanics and Mining Sciences 2007;44(7):997-1010.

Cundall PA,Hart RD.Development of generalized 2-D and 3-D distinct element programs for modeling jointed rock.Minneapolis,MN,USA:Itasca Consulting Group Inc.;1985.

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.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1988;25(3):107-16.

Damjanac B,Board M,Lin M,Kicker D,Leem J.Mechanical degradation of emplacement drifts at Yucca Mountain-a modeling case study:part II:lithophysal rock.International Journal of Rock Mechanics and Mining Sciences 2007;44(3):368-99.

Diederichs MS.Instability of hard rock masses:the role of tensile damage and relaxation.PhD Thesis.Waterloo,Ontario,USA:University of Waterloo;2000. Diederichs MS.Rock fracture and collapse under low con f i nement conditions.Rock Mechanics and Rock Engineering 2003;36(5):339-81.

Diederichs MS,Kaiser PK,Eberhardt E.Damage initiation and propagation in hard rock during tunnelling and the in f l uence of near-face stress rotation.International Journal of Rock Mechanics and Mining Sciences 2004;41(5):785-812.

Diederichs MS.The 2003 Canadian geotechnical colloquium:mechanistic interpretation and practical application of damage and spalling prediction criteria for deep tunnelling.Canadian Geotechnical Journal 2007;44(9):1082-116.

Gao FQ,Stead D.The application of a modi f i ed Voronoi logic to brittle fracture modelling at the laboratory and f i eld scale.International Journal of Rock Mechanics and Mining Sciences 2014;68:1-14.

Gatelier N,Pellet F,Loret B.Mechanical damage of an anisotropic porous rock in cyclic triaxial tests.International Journal of Rock Mechanics and Mining Sciences 2002;39(3):335-54.

Ghazvinian E,Diederichs MS,Martin CD.Identi f i cation of crack damage thresholds in crystalline rock.In:Proceedings of Eurock 2012,Stockholm,Sweden;2012.

Ghazvinian E,Perras M,Diederichs M,Labrie D.The effect of anisotropy on crack damage thresholds in brittle rocks.In:Proceedings of the 47th US Rock Mechanics/Geomechanics Symposium.San Francisco,California,USA:American Rock Mechanics Association;2013a.

Ghazvinian E,Perras M,Langford C,Diederichs MS,Labrie D.A comprehensive investigation of crack damage anisotropy in Cobourg limestone and its effect on the failure envelope.In:Proceedings of GeoMontreal 2013.Montreal,QC, Canada:Canadian Geotechnical Society;2013b.

Ghazvinian E.Modelling and testing strategies for brittle fracture simulation in crystalline rock samples.MS Thesis.Kingston,Ontario:Queen’s University; 2010.

Hajiabdolmajid V,Kaiser PK,Martin CD.Modeling brittle failure of rock.International Journal of Rock Mechanics and Mining Sciences 2002;39:731-41.

Hakala M,Kuula H,Hudson JA.Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction:a case study of the Olkiluoto mica gneiss,Finland.International Journal of Rock Mechanics and Mining Sciences 2007;44(1):14-46.

Hobbs BE,Means WD,Williams PF.An outline of structural geology,vol.570.New York:Wiley;1976.

Itasca.3DEC(3 dimensional distinct element code)version 5.0.Minneapolis,MN, USA:Itasca Consulting Group Inc.;2013a.

Itasca.UDEC(Universal distinct element code).Minneapolis,MN,USA:Itasca Consulting Group Inc.;2013b.

Itasca.FLAC(Fast Lagrangian analysis of Continua)version 7.0.Minneapolis,MN, USA:Itasca Consulting Group Inc.;2012.

Itasca.Long-term geomechanical stability analysis:OPG’s deep geological repository for low&intermediate level waste.2011.Technical Report,NWMO DGR-TR-2011-17.

Itasca.PFC3D(Particle f l ow code 3D)version 4.0.Minneapolis,MN,USA:Itasca Consulting Group Inc.;2008.

Jaeger JC,Cook NG.Fundamentals of rock mechanics.London,UK:Methuen&Co Ltd.;1969.

Kanit T,Forest S,Galliet I,Mounoury V,Jeulin D.Determination of the size of the representative volume element for random composites:statistical and numerical approach.International Journal of Solids and Structures 2003;40:3647-79.

Kazerani T,Zhao J.Micromechanical parameters in bonded particle method for modeling of brittle material failure.International Journal for Numerical and Analytical Methods in Geomechanics 2010;34(18):1877-95.

Kazerani T.A discontinuum-based model to simulate compressive and tensile failure in sedimentary rock.Journal of Rock Mechanics and Geotechnical Engineering 2013;5(5):378-88.

Lama RD,Vutukuri VS.Handbook on mechanical properties of rock,vol.II.Clausthal,Germany:Trans Tech Publications;1978.

Lan H,Martin CD,Hu B.Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading.Journal of Geophysical Research 2010;115:1-14.

Lisjak A,Tatone BS,Grasselli G,Vietor T.Numerical modelling of the anisotropic mechanical behaviour of Opalinus Clay at the laboratory-scale using FEM/DEM. Rock Mechanics and Rock Engineering 2014a;47(1):187-206.

Lisjak A,Grasselli G,Vietor T.Continuum-discontinuum analysis of failure mechanisms around unsupported circular excavations in anisotropic clay shales.International Journal of Rock Mechanics and Mining Sciences 2014b;65:96-115.

Mahabadi OK,Lisjak A,Grasselli G,Munjiza A.Y-Geo:a new combined f i nitediscrete element numerical code for geomechanical applications.International Journal of Geomechanics 2012;12(6):676-88.

Martin CD,Chandler NA.The progressive fracture of Lac du Bonnet granite.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1994;31(6):643-59.

Martin CD.Seventeenth Canadian geotechnical colloquium:the effect of cohesion loss and stress path on brittle rock strength.Canadian Geotechnical Journal 1997;34(5):698-725.

Martin CD.The strength of massive Lac du Bonnet granite around underground openings.PhD Thesis.Winnipeg,Manitoba:University of Manitoba;1993.

Mas Ivars D,Pierce ME,Darcel C,Reyes-Montes J,Potyondy DO,Young RP, Cundall PA.The synthetic rock mass approach for jointed rock mass modelling. International Journal of Rock Mechanics and Mining Sciences 2011;48(2):219-44.

Passchier CW,Trouw RAJ.Microtectonics.2nd ed.Berlin:Springer;2005.

Potyondy DO,Cundall PA.A bonded-particle model for rock.International Journal of Rock Mechanics and Mining Sciences 2004;41(8):1329-64.

Potyondy DO.A f l at-jointed bonded-particle material for hard rock.In:Proceedings of the 46th U.S.Rock Mechanics/Geomechanics Symposium.Chicago,USA: American Rock Mechanics Association;2012.

Quey R,Dawson PR,Barbe F.Large-scale 3D random polycrystals for the f i nite element method:generation,meshing and remeshing.Computer Methods in Applied Mechanics and Engineering 2011;200(17):1729-45.

Quey R,Dawson PR,Driver JH.Grain orientation fragmentation in hot-deformed aluminium:experiment and simulation.Journal of the Mechanics and Physics of Solids 2012;60(3):509-24.

Quey R.Neper:numerical descriptors of polycrystals(version 2.0).2014.Available at:http://neper.sourceforge.net/.

Scholtès L,Donzé FV.A DEM model for soft and hard rocks:role of grain interlocking on strength.Journal of the Mechanics and Physics of Solids 2013;61(2): 352-69.

Yoon JS,Jeon S,Stephansson O,Zang A,Dresen G.A new method of microparameter determination for PFC2Dsynthetic rock model generation.In:Proceedings of the 1st International FLAC/DEM Symposium on Numerical Modeling.Minneapolis,MN,USA:Itasca Consulting Group Inc.;2008.

You S,Zhao GF,Ji HG.Model for transversely isotropic materials based on distinct lattice spring model(DLSM).Journal of Computers 2011;6(6):1139-44.

Ehsan Ghazvinian graduated with a BSc in Civil Engineering from K.N.Toosi University of Technology in Tehran,Iran(2008).He pursued his graduate studies in the if eld of geomechanics at Queen’s University in Kingston, Canada where he received his MASc in Geological Engineering(2010)and is currently completing his PhD in the same fi eld(2015).His research interests include brittle fracturing of hard rocks and the application of discontinuum modelling to the design of underground excavations.He is currently involved in ISRM Spalling Commission for work on suggested methods for laboratory estimation of crack damage thresholds.

*Corresponding author.Tel.:+1 613 539 3330.

E-mail addresses:e.ghazvinian@queensu.ca,eghazvinian@gmail.com (E.Ghazvinian).

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

1674-7755?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

http://dx.doi.org/10.1016/j.jrmge.2014.09.001