Three-dimensional FDEM numerical simulation of failure processes observed in Opalinus Clay laboratory samples

Omid Mhdi,Ptrick Kifosh,Pul Mrschll,Tim Vietor

aGeomechanica Inc.,90 Adelaide St W,Suite 300,Toronto,ON,M5H 3V9,Canada

bNational Cooperative for the Disposal of Radioactive Waste(NAGRA),Wettingen,Switzerland

Three-dimensional FDEM numerical simulation of failure processes observed in Opalinus Clay laboratory samples

Omid Mahabadia,*,Patrick Kaifosha,Paul Marschallb,Tim Vietorb

aGeomechanica Inc.,90 Adelaide St W,Suite 300,Toronto,ON,M5H 3V9,Canada

bNational Cooperative for the Disposal of Radioactive Waste(NAGRA),Wettingen,Switzerland

A R T I C L E I N F O

Article history:

Received 4 September 2014

Received in revised form

16 October 2014

Accepted 20 October 2014

Available online 7 November 2014

Three-dimensional(3D)hybrid f i nitediscrete element method(FDEM)

Intermediate principal stress

Discrete element methods

True triaxial behaviour

Failure envelope

This study presents the fi rst step of a research project that aims at using a three-dimensional(3D)hybrid fi nite-discrete element method(FDEM)to investigate the development of an excavation damaged zone (EDZ)around tunnels in a clay shale formation known as Opalinus Clay.The 3D FDEM was fi rst calibrated against standard laboratory experiments,including Brazilian disc test and uniaxial compression test.The effect of increasing con fi ning pressure on the mechanical response and fracture propagation of the rock was quanti fi ed under triaxial compression tests.Polyaxial(or true triaxial)simulations highlighted the effect of the intermediate principal stress(σ2)on fracture directions in the model:as the intermediate principal stress increased,fractures tended to align in the direction parallel to the plane de fi ned by the major and intermediate principal stresses.The peak strength was also shown to vary with changingσ2. ?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

Clay shales possess favourable long-term isolation properties, and as a result,have been considered as a possible host rock for the geological disposal of radioactive waste.However,during the excavation,the isolation properties of the intact rock can be degraded as an excavation damaged zone(EDZ)forms around the underground openings.Numerical modelling has been extensively used to understand the failure mechanisms that lead to the formation of an EDZ in Opalinus Clay(a type of clay shale),and in particular,around the tunnel boundaries(e.g.Popp et al.,2008; Yong et al.,2010;Lisjak,2013;Lisjak et al.,2014a).Among others, these studies have provided valuable insights into the failure processes involved during the formation of EDZ in Opalinus Clay. Nevertheless,these studies were either performed using twodimensional(2D)methods or simplistic numerical tools.

The present study should be considered the f i rst step of a longer-term research project that aims at using three-dimensional (3D)hybrid f i nite-discrete element method(FDEM)to investigate the development of the EDZ around tunnels in Opalinus Clay.As an initial stage,the feasibility of using the FDEM tool to model 3D problems was assessed both quantitatively and qualitatively by modelling laboratory experiments in rocks,including Brazilian disc test,uniaxial,triaxial,and polyaxial or true triaxial compression tests.The 3D results were compared with 2D FDEM and experimental observations and data.

Particular emphasis was placed on analysing the in f l uence of the intermediate principal stress,σ2,on rock fracturing and strength near excavation boundaries.Using a prismatic sample,this simulation reproduced the loading conditions ofσ1≠0,σ2≠0,σ3=0 that exist at the tunnel boundary(Fig.1a).

It is anticipated that the use of 3D models will improve the quality and reliability of the simulations used to investigate the behaviour of Opalinus Clay during the excavation of tunnels and will allow to better understand how the advancement of the tunnel face in f l uences the shape and extension of the EDZ.

2.Previous modelling work

Rock mechanics laboratory experiments have been extensively modelled using continuum,discontinuum,and hybrid continuumdiscontinuum methods in 2D and 3D.A full literature review of these studies requires a dedicated article by itself,and thus,only selected publications are discussed here.

Potyondy and Cundall(2004)used PFC2D and PFC3D to model uniaxial and triaxial(biaxial in 2D)tests.Although acceptable fracture patterns were reported,the stress-strain response did not show any brittle to ductile transition even as the con f i ning pressure was increased toσ3=70 MPa.More recently,Zhang(2014)used a synthetic rock mass approach within PFC3D(Itasca,2012)andreproduced an appropriate brittle to ductile transition in triaxial and true triaxial tests.However,their results showed an arti f i cial increase of model stiffness(Young’s modulus)as the con f i ning pressure increased.Moreover,the peak strength of the true triaxial tests seemed to always increase as the con f i nement increased.

Fig.1.(a)Schematic of principal stresses around a tunnel.(b)Representation of the transversely isotropic constitutive law.

Kazerani and Zhao(2010)and Kazerani(2013)used bonded particle modelling(BPM)within UDEC(Itasca,2013a)to model Brazilian and uniaxial compression tests.Their results illustrated that as the con f i ning pressure increased,more isolated fractures developed and it became more dif f i cult to determine a plane of failure.Also,as the con f i nement increased so did the material stiffness.In addition,Gao and Stead(2014)used a modi f i ed Voronoi method within UDEC and 3DEC(Itasca,2013a,b)to model laboratory experiments.While their 2D uniaxial model exhibited a shear-dominated failure mechanism,the equivalent 3D model exhibited a more realistic,tensile-dominated mechanism.They attributed this phenomenon to the lack of porosity in the 2D sample,which was effective under a plane strain condition,thus limiting the out-of-plane deformation.Their Brazilian simulation showed a mixture of tensile splitting fractures and shear fractures closer to loading platens.

Jia et al.(2012)used RFPA3D to model polyaxial experiments and tunnel excavations in 3D.They showed that,with the increase of con f i ning pressure,the damaged elements align themselves parallel to the free surface with zero con f i nement.Scholtès and Donzé(2013)used 3D YADE discrete element method(DEM) code to model laboratory experiments and reproduced acceptable failure envelopes and stress-strain responses,including the brittle-ductile transition with increasing con f i nement.However,the bond failures did notseem to form any type of macroscopic features (fracture planes).

Elmo(2006)used 2D and 3D ELFEN FDEM code(Rock f i eld, 2003)to model cylindrical and prismatic uniaxial compression tests.The models did not appear to show a realistic fracturing process and any plane of failure.Cai(2008)used ELFEN to study the in f l uence of intermediate principal stress on the mechanical response of a rock.Their polyaxial simulations were performed under loading conditions similar to this study(σ1≠0,σ2≠0,σ3=0) and reproduced fractures parallel to maximum and intermediate principal stresses.Due to intrinsic assumptions of ELFEN,Mode II fracturing was not modelled.Also,the fracturing process,including coalescence of individual fractures into fractures planes was not reproduced in their results.More recently,Hamdi et al.(2014)used 3D ELFEN to model Brazilian and compression tests.The Brazilian disc simulations showed reasonable fracturing,while in the uniaxial compression tests it was dif f i cult to distinguish fracture planes.They captured the brittle to ductile transition as the con f i nement increased.However,the pre-peak stress-strain responses showed a perfectly linear behaviour.

3.Modelling software:FDEM

The hybrid FDEM is a numerical method which combines continuum mechanics principles with DEM algorithms to simulate multiple interacting deformable bodies(Munjiza,2004).In 3D FDEM,each solid is discretized as a mesh consisting of nodes and tetrahedral elements.An explicit time integration scheme is applied to solve the equations of motion for the discretized system and to update the nodal coordinates at each simulation time step.In general,the governing equation can be expressed as(Munjiza et al., 1995):

where M and C are the system mass and damping diagonal matrices,respectively;x is the vector of nodal displacements;Fint, Fext,and Fcare the vectors of internal resisting forces,applied external loads,and contact forces,respectively.

Contact forces,Fc,are calculated either between contacting discrete bodies or along internal discontinuities(i.e.pre-existing or newly created fractures)(Section 3.1).Internal resisting forces,Fint, include the contribution from the elastic forces,Fe,and the crack element bonding forces,Fb,which are used to simulate material elastic deformation and progressive failure,respectively,as further explained in Sections 3.2 and 3.3.

Numerical damping is introduced in the governing equation to account for energy dissipation due to non-linear material behaviour or to model quasi-static phenomena by dynamic relaxation (Munjiza,2004).The matrix C is equal to

whereμand I are the damping coef f i cient and the identity matrix, respectively.

3.1.Contact detection and interaction

An FDEM simulation can comprise a very large number of potentially interacting distinct elements.To correctly capture this behaviour,contacting couples(i.e.pairs of contacting discrete elements)must f i rst be detected.Subsequently,the interaction forces, Fc,resulting from such contacts can be de f i ned.Contact interaction forces are calculated between all pairs of elements that overlap in space.Two types of forces are applied to the elements of each contacting pair:repulsive forces and frictional forces.The repulsive forces between the elements of each contacting pair(i.e.couples) are calculated using a penalty function method(Munjiza,2004). Contacting couples tend to penetrate into each other,generating distributed contact forces,which depend on the shape and size of the overlap between the two bodies and the value of the penalty term,pn.As penalty values tend to in f i nity,a body impenetrability condition is approached.The frictional forces between contacting couples are calculated using a Coulomb-type friction law.These frictional forces are used to simulate the shear strength of intact material and of pre-existing and newly created fractures(Mahabadi et al.,2012).

3.2.Elastic behaviour

Since the material strain is expected to be localised in the crack elements(see Section 3.3),the bulk material is treated as linearelastic using constant-strain tetrahedral elements.Since Opalinus Clay shows strong transverse isotropy,the Hooke’s law based on Young’s modulus and Poisson’s ratio in the x-y symmetry plane (Ep,νp),Young’s modulus and Poisson’s ratio in the z-direction (Ez,νpz),and shear modulus in z-direction(Gzp)has been implemented in the code(Fig.1b).The stiffness matrix for transversely isotropic materials is given as(Singh,2007):

In general,the elastic response of an FDEM model depends not only on the constitutive relationship of the tetrahedral elements but also on the properties of the crack elements.Since the elastic deformation before the onset of fracturing takes place in the bulk material,no deformation should in theory occur in the crack elements before the intrinsic strength is exceeded.However,a f i nite stiffness is required for the crack elements by the explicit time integration scheme of FDEM.Such an arti f i cial compliance is represented by the normal,tangential,and fracture penalty values,pn, pt,and pf,for compressive,shear,and tensile loading conditions, respectively.For practical purposes,the contribution of crack elements to the overall model compliance can be largely limited by adopting very high penalty values(Mahabadi et al.,2012).Consequently,the elastic response is effectively controlled by the choice of the stress-strain relationship of the continuum tetrahedral elements.

3.3.Material failure

The progressive failure of rock material is modelled using a cohesive-zone approach,a technique f i rst introduced in the context of the elasto-plastic fracturing of ductile metals(Dugdale,1960; Barenblatt,1962).This approach aims at capturing the non-linear interdependence between stresses and strains that characterises the zone ahead of a macro-crack tip known as the fracture process zone(FPZ).As depicted in Fig.2,the FPZ in brittle rocks manifests itself in the form of micro-cracking and interlocking related to the presence of micro-scale inhomogeneities(e.g.mineral grains and pre-existing defects or voids)(Labuz et al.,1987).When using cohesive-zone models,the failure of the material progresses based solely on the strength degradation of dedicated interface elements (referred herein to as crack elements)and therefore emerges as a natural outcome of the deformation process without employing any additional macroscopic failure criterion.

In the FDEM code,the bonding stresses transferred by the material are decreasing functions of the displacement discontinuity across the crack elements according to the cohesive laws illustrated in Fig.3b.These constitutive relationships represent a modi f i ed version of the crack model response proposed by Munjiza et al. (1999).Mode I(i.e.opening)and Mode II(i.e.shearing)fracturing are simulated by a cohesive model based on the FPZ model originally developed for concrete by Hillerborg et al.(1976)under the name of f i ctitious crack model.A Mode I fracture is assumed to initiate when the crack tip opening,o,reaches a critical value,op, related to the intrinsic tensile strength of the material,ft.When the crack opens,the normal bonding stress,σ,is not assumed to fall to zero at once but to decrease with increasing crack opening until a residual opening value,or,is reached and a stress-free surface is created(i.e.σ=0).Mode II(i.e.shearing)fracturing is simulated by a slip-weakening model conceptually similar to that of Ida(1972).For shear loading conditions,a tangential bonding stress,τ,exists between the two fracture walls,which is a function of the amount of slip,s,and the normal stress on the fracture,σn.The critical slip,sp, corresponds to the intrinsic shear strength of the rock,fs,de f i ned by the Mohr-Coulomb failure criterion as

where c is the material cohesion,φis the internal friction angle of the material,andσnis the normal stress acting across the fracture surfaces.Upon undergoing the critical slip,sp,the tangential bonding stress is gradually reduced to a residual value,fr,which corresponds to a purely frictional resistance:

Fig.2.Cohesive-zone approach for material failure modelling in FDEM.(a)Conceptual model of a tensile crack in a heterogeneous rock material(modi f i ed after Labuz et al.(1987)). (b)Theoretical FPZ model of Hillerborg et al.(1976).(c)Cross-section of FDEM implementation of the FPZ using tetrahedral elastic elements and six-noded crack elements to represent the bulk material and the fracture,respectively.Triangles are shrunk for illustration purposes.

In the current crack element implementation,the unloading path in the softening branch coincides with the loading path. Therefore,the model is strictly valid for monotonic loading conditions only.

For mixed mode fracturing,the rupture of a crack element is de f i ned by the following coupling criterion between crack opening and slip:

The external energy required to fully break a unit surface area of cohesive crack corresponds to the input speci f i c fracture energy,Gc. Gcis de f i ned in terms of the material properties GIcand GIIcwhich correspond to the strain energy release rates for Mode I and Mode IIfracturing,respectively.The crack residual displacement values,orand sr,are such that:

Fig.3.Simulation of fracturing with the FDEM code.(a)Representation of a cross-section of a continuum using cohesive crack elements interspersed throughout a mesh of tetrahedral elastic elements.Triangles are shrunk for illustration purposes.(b)Constitutive behaviour of the crack elements de f i ned in terms of normal and tangential bonding stresses,σandτ,vs.crack relative displacements,o and s(i.e.opening and sliding).

Following an approach similar to that pioneered by Xu and Needleman(1996),the crack elements in FDEM are interspersed throughout the material(i.e.across the edges of all tetrahedral element pairs)from the very beginning of the simulation.Thus, cracks are allowed to nucleate and grow without any additional assumption or criterion other than the crack element constitutive response.Upon breakage of the cohesive surface,the crack element is removed from the simulation and therefore the model transition from a continuum to discontinuum is locally completed.The newly created discontinuity is treated by the contact algorithm through the contact forces,Fc,brie f l y described in the previous section.As the simulation progresses,f i nite displacements and rotations of discrete bodies are allowed and new contacts are automatically recognised.

4.Modelling methodology

4.1.Geometry and boundary conditions

As the model geometries in Fig.4 show,the 3D laboratory-scale models included a 17.4 mm(diameter)?35.8 mm(height)cylindrical specimen for the uniaxial and triaxial compression tests,a 29.85 mm(diameter)?14.85 mm(thickness)circular disc specimen for the Brazilian test,and a prismatic specimen with a square cross-section of 8.7 mm?8.7 mm and a height of 35.8 mm for the polyaxial(true triaxial)tests.The specimens were meshed with a uniform,unstructured tetrahedral grid of nominal element size of h=2.0 mm(Fig.4).The top and bottom of the Brazilian disc were fl attened over the width of an element to avoid point loading (Fig.4a).The equation of motion for the discretized system,Eq.(1), was integrated with a time step of 5?107ms;this value was the largest time step size that ensured numerical stability for the explicit time integration scheme of the code.Uniaxial loading conditions were obtained by means of two rigid platens moving in opposite directions.As shown in Fig.4d,the loading velocity linearly increased from zero over t1=30,000 time steps to a constant value of v1=1 m/s.This gradual application of the velocity minimised dynamic artefacts in the system.The FDEM graphical user interface Y-GUI(Mahabadi et al.,2010)was used to assign boundary conditions and material properties to the model.

4.2.Input parameters

As reported in Table 1,input parameters were based on material properties previously calibrated by Lisjak et al.(2014b).Since the current 3D formulation does not account for strength anisotropy,the calibrated strength parameters were averaged and used as inputs for the 3D simulations(Table 1).A friction coef f i cient equal to 0.1 was assumed at the platen-sample interfaces.The effect of the crack element compliance on the overall model stiffness was minimised by selecting appropriate values for the penalty coef f i cients.Based on the recommendations of Mahabadi et al.(2012),normal penalty,pn, tangential penalty,pt,and fracture penalty,pf,were set equal to 10, 10,and 5 times of the average Young’s modulus,E,respectively.As illustrated in Fig.1b,the x-y plane is assumed to be the plane of symmetry for the transversely isotropic constitutive law of Eq.(3), representing the P-sample of Lisjak(2013).

5.Results

In this section,the results of the Brazilian disc test,uncon f i ned uniaxial compression test,the triaxial and polyaxial(true triaxial) compression tests at different con f i ning pressures are discussed. Note that throughout this manuscript,compressive stresses are negative.

5.1.Brazilian disc test

The stress evolution(σzz)and fracturing process of the Brazilian disc are illustrated in Fig.5.As the specimen is loaded,tensile stresses form in the centre of the disc.However,high shear stresses also develop closer to the loading platens.Thus,the specimen experiences a mix of tensile failure along the loading direction and crushing failure closer to the loading points.Similar trends in 3D simulations have been reported by other researchers(Gao and Stead,2014).

Fig.4.3D meshes of(a)the Brazilian disc test,(b)uniaxial and triaxial compression tests,and(c)polyaxial tests.(d)Loading was applied by prescribed velocities(v)that followed a linear increase to v1over t1time steps.Dimensions are given in units of millimetres.

Table 1Input parameters of the FDEM model calibrated based on laboratory-scale rock mechanics tests on Opalinus Clay(Lisjak et al.,2014b).

It is argued that tensile failure would dominate over shear failure if a lower Mode I fracture energy(GIc)was used(Table 1). This phenomenon is clearly demonstrated when the value was reduced to GIc=0.02 J/m2(Fig.6).Note that GIcis related to Mode I fracture toughness(KIc)through Young’s modulus.However,if KIcis not available,GIcshould be calibrated via an iterative process as depicted by Tatone(2014).

The indirect tensile strength of the Brazilian disc,calculated as

is shown in Fig.7(with P as the reaction force,D as the disc diameter,and B as the disc thickness).As expected,the stresses increased linearly until reaching about 80%of the peak(strain hardening),where the linearity was lost.Upon reaching the peak stress,brittle failure was followed during which the stresses decreased in the post-peak(strain softening).The tensile strength was slightly overestimated compared to the laboratory experiments with a tensile strength ranging between 0.67 MPa and 1.5 MPa.

5.2.Uniaxial and triaxial compression tests

Triaxial simulations were performed using variableσ3values of 1 MPa,2.5 MPa,5 MPa,7.5 MPa,10 MPa,and 12.5 MPa.Similar to the loading velocity,the con f i ning pressure was applied gradually over 30,000 time steps(Fig.4d).

5.2.1.Fracturing process

Fig.5.Stress(σzz)and the fracture pattern of the Brazilian specimen along y-z plane during the simulation time at(a)0.05 ms,(b)0.1 ms,(c)0.15 ms,(d)0.2 ms,(e)0.25 ms,and(f) 0.3 ms after the start of simulation,respectively.Compressive stresses are negative.The bottom right inset shows where the cross-sections were taken.

Fig.8 demonstrates the major principal stress(σ1)and the fracture pattern of the uniaxial specimen along two orthogonalviews(x-y and y-z planes)during the simulation time.As the specimen was gradually loaded(Fig.8a-d),isolated fractures formed in various locations in the specimen(Fig.8e-h).These fractures later coalesced to form more macroscopic fractures (Fig.8i-l).Fracturing of the triaxial specimen atσ3=1 MPa followed a similar trend as the uniaxial model.

X-ray style 3D fracture volume evolution of the specimen during the simulation time is presented in Fig.9.The gray-scale intensity in these images represents fractures localised into planes,i.e. fracture density.As shown in Fig.9d,the major planes of failure were aligned approximately at 52?from the x-z plane.This value corresponds well with the theoretical value according to the Mohr-Coulomb criterion:π/4+φ/2=56?.

The application of more con f i ning pressure(σ3?2.5 MPa) limited the growth of tensile fractures parallel to the major principal stress direction(parallel to y-axis),thus resulting in a more shear-dominated fracturing process.Unlike the models at low con f i nement(σ3=0-1 MPa),in which one or two major macroscopic fractures developed,the fractures in triaxial models with σ3?2.5 MPa formed shear bands or faults(Figs.10 and 11).As the con f i ning pressure exceeded 5 MPa,for instance,in the case of σ3=12.5 MPa in Figs.12 and 13,shear failures formed a network of small fractures that continued to deform plastically.As a result,no noticeable plane of failure was formed and a more ductile postyield stress response was captured(see Section 5.2.2).These types of fracturing behaviour correspond well with the published literature of experimental observations(Mogi,2007;Mahabadi, 2012)and of previous 2D FDEM simulations(Mahabadi,2012; Lisjak et al.,2014b).

5.2.2.Mechanical response

The stress-strain curves of the uniaxial and triaxial simulations are presented in Fig.14(major principal stress,σ1,vs.major principal strain,ε1).At low or no con f i ning pressure(σ3=0-2.5 MPa), the specimen showed strain hardening until reaching the peak strength,which was then followed by brittle strain softening, during which the stresses decreased rapidly.As the con f i ning pressure increased,the peak strength increased and a more ductile post-peak behaviour was observed(σ3=5 MPa).At higher con f i ning pressures(σ3＞5 MPa),a completely ductile post-peak was simulated.Overall,this trend agrees with laboratory f i ndings and previous numerical simulations performed in 2D(Lisjak et al., 2014b).It should be noted that the strains in Fig.14 were calculated from platen displacements divided by specimen height.

The stress-strain curves seemed to depart from linearity at about 60%-70%of peak stress(Fig.14).This range coincides with the so-called crack damage stress,which marks the beginning of unstable cracking(Bieniawski,1967;Martin and Chandler,1994).In the case of an FDEM simulation,this stage is represented by softening of individual crack elements.

Fracturing in rocks usually occurs parallel to the major principal stress(σ1)direction.Therefore,applying any con f i ning pressure tends to arrest the initiation and development of fractures.For this reason,the rock exhibits higher strength,possibly strain hardening, and less fracturing occurs as the con f i ning pressure increases.

For comparison,2D uniaxial and biaxial models were built using similar element sizes and loading rates as the 3D models(Fig.15). Unlike the 3D models,the 2D models showed a perfectly linear elastic mechanical response up until yielding,followed by a sudden drop in the post-peak stress at low con f i ning pressures(σ3=0-5 MPa)or by strain hardening at higher pressures(σ3?7.5 MPa). Both the 2D and 3D models captured the increase in strength as the con f i nement increased.

Using the results of the Brazilian,uniaxial,and triaxial(biaxial in 2D)tests,the failure envelope of the simulations was derived (Fig.16).While both the 3D and 2D failure envelopes reproduced the increase in peak strength as a function of increasing con f i ning pressure,the 3D envelope signi f i cantly underestimated the slope of the curve,which corresponds to the emergent internal friction angle(Table 2).However,the simulated cohesion was relatively close to the input value.Since the model inputs,in particular, penalty terms were adopted from calibrated 2D simulations, further investigation is required to assess the validity of these parameters in 3D.

Fig.6.Stress(σzz)and the fracture pattern of the Brazilian specimen with reduced Mode I fracture energy along y-z plane during the simulation time at(a)0.05 ms,(b)0.0525 ms, (c)0.055 ms,and(d)0.0575 ms after the start of simulation,respectively.

Fig.7.Indirect tensile stress vs.platen displacement of the 3D Brazilian disc simulation.

Fig.8.Major principal stress(σ1)and the fracture pattern of the uniaxial specimen along two orthogonal views(x-y and y-z planes)during the simulation time at(a)0 ms,(b) 0.0625 ms,(c)0.125 ms,(d)0.1875 ms,(e)0.199 ms,(f)0.2015 ms,(g)0.2065 ms,(h)0.2115 ms,(i)0.2165 ms,(j)0.2215 ms,(k)0.2265 ms,and(l)0.2315 ms after the start of simulation,respectively.Compressive stresses are negative.

Fig.9.3D fracture volume of the uniaxial specimen during the simulation time at(a)0.2015 ms,(b)0.2115 ms,(c)0.2215 ms,and(d)0.2315 ms after the start of simulation, respectively.The dashed lines in(d)represent approximate failure plane orientations.

Fig.10.Major principal stress(σ1)and the fracture pattern of the triaxial specimen withσ3=5 MPa along two orthogonal views(x-y and y-z planes)during the simulation time at (a)0 ms,(b)0.0625 ms,(c)0.125 ms,(d)0.1875 ms,(e)0.2355 ms,(f)0.2505 ms,(g)0.2655 ms,(h)0.2805 ms,(i)0.2955 ms,(j)0.3105 ms,(k)0.3255 ms,and(l)0.3375 ms after the start of simulation,respectively.Compressive stresses are negative.

Finally,the 3D models captured the elastic response of the material very well:the emergent Young’s modulus showed a merely 5%divergence from the input value(Table 2).

5.3.Polyaxial compression tests

Polyaxial or true triaxial simulations were performed using σ3=0 MPa with variableσ2values of 0 MPa,1 MPa,2.5 MPa,5 MPa, 7.5 MPa,10 MPa,and 12.5 MPa.Initial simulation runs revealed that the models withσ2=10-12.5 MPa would fracture mainly due to the rapid application ofσ2.Therefore,unlike the triaxial models where the con f i ning pressure was applied gradually over 30,000 time steps(Fig.4d),the intermediate principal stress in the polyaxial models was applied gradually at a rate of 80 MPa/ms until reaching the prescribed,constant value.

5.3.1.Fracturing process

The fracture propagation of the polyaxial model withσ2=0 MPa is demonstrated in Figs.17 and 18.Similar to the uniaxial and triaxial models(Figs.8 and 9),the formation of two distinct fracture planes parallel toσ2can be observed.A similar behaviour was noticed forσ2=1 MPa and 2.5 MPa.

As the intermediate principal stress increased,σ2?5 MPa,the angle between the faults and the maximum principal stress,σ1,direction was signi f i cantly reduced(Figs.19 and 20).As a result,two distinct fracture planes can be easily distinguished in Fig.20 where the fractures are parallel toσ2and are mainly in the direction ofσ1. This f i nding agrees with experimental observations(Mogi,1967). Also,this kind of fracturing is similar to spalling-type fractures observed around underground openings.Similar to the triaxial simulations,at higher intermediate principal stresses(σ2=10-12.5 MPa),more isolated fractures formed and it became harder to distinguish well-de f i ned fracture surfaces(Figs.21 and 22).

5.3.2.Mechanical response

Fig.11.3D fracture volume of the triaxial specimen withσ3=5 MPa during the simulation time at(a)0.2505 ms,(b)0.2805 ms,(c)0.3105 ms,and(d)0.3375 ms after the start of simulation,respectively.

Fig.12.Major principal stress(σ1)and the fracture pattern of the triaxial specimen withσ3=12.5 MPa along two orthogonal views(x-y and y-z planes)during the simulation time at(a)0 ms,(b)0.125 ms,(c)0.25 ms,(d)0.2925 ms,(e)0.3175 ms,(f)0.3425 ms,(g)0.3675 ms,(h)0.3925 ms,(i)0.4175 ms,(j)0.4425 ms,(k)0.4675 ms,and(l)0.4925 ms after the start of simulation,respectively.Compressive stresses are negative.

The stress-strain response of the polyaxial models and the peak maximum principal stress versus the intermediate principal stress (σ1vs.σ2)are presented in Fig.23.The addition ofσ2initiallyincreased the peak stress(σ1),reaching the peak increase at σ2=5 MPa.In contrast,at higherσ2values(σ2?7.5),σ1decreased (Fig.23 right).Theσ1-σ2relationship showed a downward concave shape.This general trend is in good agreement with published experimental results(Mogi,1967;Haimson and Chang,2000; Haimson,2006;Colmenares and Zoback,2002;Paterson and Wong,2005;Chang and Haimson,2012)and rock failure criteria that incorporateσ2.(Wiebols and Cook,1968;Benz and Schwab, 2008).

Fig.13.3D fracture volume of the triaxial specimen withσ3=12.5 MPa during the simulation time at(a)0.3175 ms,(b)0.3675 ms,(c)0.4175 ms,and(d)0.4675 ms after the start of simulation,respectively.

Fig.14.Stress-strain curves of the 3D uniaxial and triaxial compression test simulations.

Fig.15.Stress-strain curves of the 2D uniaxial and biaxial compression test simulations.

Fig.16.Failure envelope of the 3D Brazilian,uniaxial,and triaxial compression simulations for 3D and 2D simulations(Lisjak et al.,2014b)and experiments(Popp et al., 2008;NAGRA,2008).Error bars correspond to standard deviation of the laboratory experiments.

Table 2Comparison of emergent properties in 3D simulations,2D simulations,and experiment(used as model input).The percent values in parentheses show the difference of the variable with respect to the experiment.

Fig.17.Major principal stress(σ1)and the fracture pattern of the polyaxial specimen withσ2=0 MPa along two orthogonal views during the simulation time at(a)0 ms,(b) 0.05 ms,(c)0.1 ms,(d)0.15 ms,(e)0.18 ms,(f)0.19 ms,(g)0.2 ms,(h)0.21 ms,(i)0.2185 ms,(j)0.227 ms,(k)0.2355 ms,and(l)0.244 ms after the start of simulation,respectively. Compressive stresses are negative.

Nonetheless,comparison with literature indicates that the extent of the change may be underestimated by the code(similar to the triaxial models).Using another FDEM code,Cai(2008)also underestimated the in f l uence of the intermediate principal stress. Mogi(1967,1971)among others showed that the increase of peak stress is more pronounced at higherσ3values and for brittle rather than ductile material.Therefore,to further study the effect ofσ2on peak strength,various values ofσ3should be used(versusσ3=0MPa used in this work).However,this was outside the scope of the current study.

Fig.18.3D fracture volume of the polyaxial specimen withσ2=0 MPa during the simulation time at(a)0.2 ms,(b)0.21 ms,(c)0.227 ms,and(d)0.244 ms after the start of simulation,respectively.

Since the Mohr-Coulomb failure criterion,which is used in our FDEM code,ignores the effect of the intermediate principal stress, the peak strength should be insensitive toσ2.However,the code successfully captured the dependence of peakσ1onσ2,albeit to a lower degree,as an emergent property,meaning that this behaviour was not a result of the constitutive relationships implemented in the model but an emerging property of the simulations.Regardless,to capture the in f l uence of the intermediate principal stress in full,a failure criterion that incorporatesσ2in its formulation should be used(e.g.Drucker-Prager,Wiebols and Cook,or modi f i ed Mohr-Coulomb formulations consideringσ2(Singh et al.,2011)).Pan et al.(2012)argued that,if using Mohr-Coulomb,σ2dependence is also related to rock heterogeneity. However,their heterogeneous simulations using Mohr-Coulomb still underestimatedσ2dependence compared with simulations using Drucker-Prager criterion.

Fig.19.Major principal stress(σ1)and the fracture pattern of the polyaxial specimen withσ2=5 MPa along two orthogonal views during the simulation time at(a)0.0625 ms,(b) 0.125 ms,(c)0.174 ms,(d)0.1865 ms,(e)0.199 ms,(f)0.2115 ms,(g)0.224 ms,(h)0.2365 ms,(i)0.249 ms,(j)0.2615 ms,(k)0.274 ms,and(l)0.2865 ms after the start of simulation, respectively.Compressive stresses are negative.

Fig.20.3D fracture volume of the polyaxial specimen withσ2=5 MPa during the simulation time at(a)0.1865 ms,(b)0.2115 ms,(c)0.2365 ms,and(d)0.2865 ms after the start of simulation,respectively.

Peculiarly,the specimens withσ2=10-12.5 MPa showed a higher stiffness(Fig.23).The reason for this artefact remains unclear and requires further research.

Fig.21.Major principal stress(σ1)and the fracture pattern of the polyaxial specimen withσ2=10 MPa along two orthogonal views during the simulation time at(a)0.075 ms,(b) 0.1 ms,(c)0.125 ms,(d)0.15 ms,(e)0.175 ms,(f)0.2 ms,(g)0.225 ms,and(h)0.25 ms after the start of simulation,respectively.Compressive stresses are negative.

6.Concluding remarks

The 3D FDEM was used to study the mechanical response and damage evolution of Opalinus Clay during classic laboratory tests.To better reproduce the response of the rock,a transversely isotropic formulation for elastic deformation was incorporated in the code.

The Brazilian disc simulation showed reasonable fracture trajectories,while the tensile strength was slightly larger than the range of experimental values.The uniaxial tests showed a realistic fracturing process with fractures coalescing into noticeable fracture planes.The effect of increasing con f i ning pressure was captured in both the brittle-ductile transition of the stress-strain response and fracture patterns.Similar to the 2D simulations(Lisjak et al., 2014b),the failure envelope of the 3D simulations tended to be in the lower range of the experimental values as the con f i ning pressure increased.However,unlike the 2D simulations,the 3D resultsunderestimated the friction angle of the rock.The choice of numerical parameters requires a comprehensive calibration study.

Fig.22.3D fracture volume of the polyaxial specimen withσ2=10 MPa during the simulation time at(a)0.175 ms,(b)0.2 ms,(c)0.225 ms,and(d)0.25 ms after the start of simulation,respectively.

The polyaxial(true triaxial)simulations successfully illustrated the in f l uence of the intermediate principal stress(σ2)on the fracturing process of the rock.It was shown that the major planes of failure will be parallel to theσ2direction.An increase inσ2values would align the fractures more parallel to theσ1direction,which can be related to spalling-type failure on tunnel sidewalls.The polyaxial results also exhibited the dependence of the peak strength on increasingσ2,following a downward concave shape where the peak increase was observed atσ2=5 MPa.Further investigation,in particular using variousσ3values,is required to study this effect.Overall,the 3D simulations presented promising results for the use of 3D FDEM in modelling the mechanical response and damage evolution of the rock.

Fig.23.Stress-strain curves of the 3D polyaxial compression test simulations.

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.

Barenblatt G.The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics 1962;7:55-129.

Benz T,Schwab R.A quantitative comparison of six rock failure criteria.International Journal of Rock Mechanics and Mining Sciences 2008;45(7):1176-86.

Bieniawski ZT.Mechanism of brittle fracture of rock:Parts I,II,and III.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1967;4(4):395-430.

Cai M.In f l uence of intermediate principal stress on rock fracturing and strength near excavation boundaries-Insight from numerical modeling.International Journal of Rock Mechanics and Mining Sciences 2008;45(5):763-72.

Chang C,Haimson B.A failure criterion for rocks based on true triaxial testing.Rock Mechanics and Rock Engineering 2012;45(6):1007-10.

Colmenares LB,Zoback MD.A statistical evaluation of intact rock failure criteria constrained by polyaxial test data for f i ve different rocks.International Journal of Rock Mechanics and Mining Sciences 2002;39(6):695-729.

Dugdale DS.Yielding of steel sheets containing slits.Journal of the Mechanics and Physics of Solids 1960;8(2):100-4.

Elmo D.Evaluation of a hybrid FEM/DEM approach for determination of rock mass strength using a combination of discontinuity mapping and fracture mechanics modelling,with particular emphasis on modelling of jointed pillars.PhD Thesis. Exeter,UK:Camborne School of Mines,University of Exeter;2006.

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.

Haimson B,Chang C.A new true triaxial cell for testing mechanical properties of rock,and its use to determine rock strength and deformability of westerly granite.International Journal of Rock Mechanics and Mining Sciences 2000;37(1-2):285-96.

Haimson B.True triaxial stresses and the brittle fracture of rock.Pure and Applied Geophysics 2006;163(5-6):1101-30.

Hamdi P,Stead D,Elmo D.Damage characterization during laboratory strength testing:a 3D-f i nite-discrete element approach.Computers and Geotechnics 2014;60:33-46.

Hillerborg A,Modéer M,Petersson PE.Analysis of crack formation and crack growth in concrete by means of fracture mechanics and f i nite elements.Cement and Concrete Research 1976;6(6):773-81.

Ida Y.Cohesive force across the tip of a longitudinal-shear crack and Grif f i th’s speci f i c surface energy.Journal of Geophysical Research 1972;77(20):3796-805.

Itasca.PFC3D(particle f l ow code in 3 dimensions).Minneapolis,USA:Itasca Consulting Group Inc.;2012.

Itasca.UDEC,universal distinct element code.Minneapolis,USA:Itasca Consulting Group Inc.;2013a.

Itasca.3DEC,three dimensional distinct element code.Minneapolis,USA:Itasca Consulting Group Inc.;2013b.

Jia P,Yang TH,Yu QL.Mechanism of parallel fractures around deep underground excavations.Theoretical and Applied Fracture Mechanics 2012;61:57-65.

Kazerani T.Effect of micromechanical parameters of microstructure on compressive and tensile failure process of rock.International Journal of Rock Mechanics and Mining Sciences 2013;64:44-55.

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

Labuz JF,Shah SP,Dowding CH.The fracture process zone in granite:evidence and effect.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1987;24(4):235-46.

Lisjak A.Investigating the in f l uence of mechanical anisotropy on the fracturing behaviour of brittle clay shales with application to deep geological repositories. PhD Thesis.Toronto,Canada:University of Toronto;2013.

Lisjak A,Grasselli G,Vietor T.Numerical investigation of failure mechanisms around unsupported circular excavations in clay shales.International Journal of Rock Mechanics and Mining Sciences 2014a(in press).

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.Investigating the in f l uence of micro-scale heterogeneity and microstructure on the failure and mechanical behaviour of geomaterials.PhD Thesis.Toronto,Canada:University of Toronto;2012.

Mahabadi OK,Grasselli G,Munjiza A.Y-GUI:a graphical user interface and preprocessor for the combined f i nite-discrete element code,Y2D,incorporating material heterogeneity.Computers and Geosciences 2010;36(2):241-52.

Mahabadi OK,Lisjak A,Munjiza A,Grasselli G.Y-Geo:a new combined f i nitediscrete element numerical code for geomechanical applications.International Journal of Geomechanics 2012;12: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.

Mogi K.Effect of the intermediate principal stress on rock failure.Journal of Geophysical Research 1967;72(20):5117-31.

Mogi K.Fracture and f l ow of rocks under high triaxial compression.Journal of Geophysical Research 1971;76(5):1255-69.

Mogi K.Experimental rock mechanics.Leiden,Netherlands:Taylor&Francis/A.A. Balkema;2007.

Munjiza A,Owen DRJ,Bicanic N.A combined f i nite-discrete element method in transientdynamicsoffracturingsolids.Engineering Computations1995;12(2):145-74.

Munjiza A.The combined f i nite-discrete element method.Chichester,UK:John Wiley&Sons;2004.

Munjiza A,Andrews KFR,White JK.Combined single and smeared crack model in combined f i nite-discrete element analysis.International Journal for Numerical Methods in Engineering 1999;44(1):41-57.

National Cooperative for the Disposal of Radioactive Waste(NAGRA).RA experiment. Updated review of the rock mechanics properties of the Opalinus Clay of the Mont Terri URL based on laboratory and f i eld testing.Technical Report TR 2008-04.National Cooperative for the Disposal of Radioactive Waste(NAGRA);2008.

Pan PZ,Feng XT,Hudson JA.Study of failure and scale effects in rocks under uniaxial compression using 3D cellular automata.International Journal of Rock Mechanics and Mining Sciences 2009;46(4):674-85.

Pan PZ,Feng XT,Hudson JA.The in f l uence of the intermediate principal stress on rock failure behaviour:a numerical study.Engineering Geology 2012;124(4):109-18.

Paterson MS,Wong TF.Experimental rock deformation-the brittle f i eld.2nd ed. New York,USA:Springer;2005.

Popp T,Salzer K,Minkley W.In f l uence of bedding planes to EDZ-evolution and the coupled HM properties of Opalinus Clay.Physics and Chemistry of the Earth, Parts A/B/C 2008;33(Supp.1):S374-87.

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

Rock f i eld Software Ltd.ELFEN.Swansea,UK:Rock f i eld Software Ltd.;2003.

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.

Singh A.Mechanics of solids.Horsham,UK:PHI Learning Ltd.;2007.

Singh M,Raj A,Singh B.Modi f i ed Mohr-Coulomb criterion for non-linear triaxial and polyaxial strength of intact rocks.International Journal of Rock Mechanics and Mining Sciences 2011;48(4):546-55.

Tatone BSA.Investigating the evolution of rock discontinuity asperity degradation and void space morphology under direct shear.PhD Thesis.Toronto,Canada: University of Toronto;2014.

Wiebols GA,Cook NGW.An energy criterion for the strength of rock in polyaxial compression.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1968;5(6):529-49.

Xu X,Needleman A.Numerical simulations of dynamic crack growth along an interface.International Journal of Fracture 1996;74(4):289-324.

Yong S,Kaiser PK,Loew S.In f l uence of tectonic shears on tunnel-induced fracturing. International Journal of Rock Mechanics and Mining Sciences 2010;47(6):894-907.

Zhang Y.Modelling hard rock pillars using a synthetic rock mass approach.PhD Thesis.Burnaby,Canada:Simon Fraser University;2014.

Omid Mahabadi holds a PhD degree in civil engineering from the University of Toronto specialising in computational geomechanics.For his PhD thesis,he investigated the in f l uences of micro-scale heterogeneity and microcracks on the failure behaviour of a crystaline rock using laboratory and numerical tools.He uses his knowledge of rock mechanics and numerical modelling by means of hybrid continuum-discontinuum numerical methods to analyse phenomena involving rock fracturing.

*Corresponding author.Tel.:+1 647 478 9767.

E-mail address:omid.mahabadi@geomechanica.com(O.Mahabadi).

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