Scholarly article on topic 'Fracture development around deep underground excavations: Insights from FDEM modelling'

Fracture development around deep underground excavations: Insights from FDEM modelling Academic research paper on "Earth and related environmental sciences"

CC BY-NC-ND
0
0
Share paper
Keywords
{Tunnelling / Caverns / "Rock fracturing" / "Excavation damaged zone (EDZ)" / "Hybrid finite-discrete element method (FDEM)" / "Numerical modelling"}

Abstract of research paper on Earth and related environmental sciences, author of scientific article — Andrea Lisjak, Daniel Figi, Giovanni Grasselli

Abstract Over the past twenty years, there has been a growing interest in the development of numerical models that can realistically capture the progressive failure of rock masses. In particular, the investigation of damage development around underground excavations represents a key issue in several rock engineering applications, including tunnelling, mining, drilling, hydroelectric power generation, and the deep geological disposal of nuclear waste. The goal of this paper is to show the effectiveness of a hybrid finite-discrete element method (FDEM) code to simulate the fracturing mechanisms associated with the excavation of underground openings in brittle rock formations. A brief review of the current state-of-the-art modelling approaches is initially provided, including the description of selecting continuum- and discontinuum-based techniques. Then, the influence of a number of factors, including mechanical and in situ stress anisotropy, as well as excavation geometry, on the simulated damage is analysed for three different geomechanical scenarios. Firstly, the fracture nucleation and growth process under isotropic rock mass conditions is simulated for a circular shaft. Secondly, the influence of mechanical anisotropy on the development of an excavation damaged zone (EDZ) around a tunnel excavated in a layered rock formation is considered. Finally, the interaction mechanisms between two large caverns of an underground hydroelectric power station are investigated, with particular emphasis on the rock mass response sensitivity to the pillar width and excavation sequence. Overall, the numerical results indicate that FDEM simulations can provide unique geomechanical insights in cases where an explicit consideration of fracture and fragmentation processes is of paramount importance.

Academic research paper on topic "Fracture development around deep underground excavations: Insights from FDEM modelling"

Accepted Manuscript

Fracture development around deep underground excavations: Insights from FDEM modelling

Andrea Lisjak, Daniel Figi, Giovanni Grasselli

PII: S1674-7755(14)00085-7

DOI: 10.1016/j.jrmge.2014.09.003

Reference: JRMGE 113

To appear in: Journal of Rock Mechanics and Geotechnical Engineering

Received Date: 5 August 2014

Revised Date: 29 September 2014 Accepted Date: 30 September 2014

Rock Mechanics and

Geotechnical

Engineering

—" —

Please cite this article as: Lisjak A, Figi D, Grasselli G, Fracture development around deep underground excavations: Insights from FDEM modelling, Journal of Rock Mechanics and Geotechnical Engineering (2014), doi: 10.1016/j.jrmge.2014.09.003.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Journal of Rock Mechanics and Geotechnical Engineering

Journal online: www.rockgcotcch.org

Fracture development around deep underground excavations: Insights from FDEM modelling

Andrea Lisjak^*, Daniel Figib , Giovanni Grassellib

a Geomechanica Inc., 90Adelaide Street West, Suite 300, M5H3V9, Toronto, ON, Canada b Department of Civil Engineering, University of Toronto, 35 St. George Street, M5S1A4, Toronto, ON, Canada Received 5 August 2014; received in revised form 29 September 2014; accepted 30 September 2014

Abstract: Over the past twenty years, there has been a growing interest in the development of numerical models that can realistically capture the progressive failure of rock masses. In particular, the investigation of damage development around underground excavations represents a key issue in several rock engineering applications, including tunnelling, mining, drilling, hydroelectric power generation, and the deep geological disposal of nuclear waste. The goal of this paper is to show the effectiveness of a hybrid finite-discrete element method (FDEM) code to simulate the fracturing mechanisms associated with the excavation of underground openings in brittle rock formations. A brief review of the current state-of-the-art modelling approaches is initially provided, including the description of selecting continuum- and discontinuum-based techniques. Then, the influence of a number of factors, including mechanical and in situ stress anisotropy, as well as excavation geometry, on the simulated damage is analysed for three different geomechanical scenarios. Firstly, the fracture nucleation and growth process under isotropic rock mass conditions is simulated for a circular shaft. Secondly, the influence of mechanical anisotropy on the development of an excavation damaged zone (EDZ) around a tunnel excavated in a layered rock formation is considered. Finally, the interaction mechanisms between two large caverns of an underground hydroelectric power station are investigated, with particular emphasis on the rock mass response sensitivity to the pillar width and excavation sequence. Overall, the numerical results indicate that FDEM simulations can provide unique geomechanical insights in cases where an explicit consideration of fracture and fragmentation processes is of paramount importance.

Keywords: tunnelling; caverns; rock fracturing; excavation damaged zone (EDZ); hybrid finite-discrete element method (FDEM); numerical modelling

1. Introduction

The stability of deep underground excavations is a common issue in a variety of rock engineering fields, including mining, tunnelling, hydroelectric power generation, and nuclear waste disposal. Furthermore, the deformation and failure of underground openings, such as boreholes, are of great importance in the drilling industry associated with hydrocarbon extraction and geothermal production. In tunnelling and mining operations, the stability of underground openings directly affects the choice of the excavation method and sequence, as well as the design of support and reinforcement measures. In the case of underground hydroelectric power stations, the rock mass behaviour is strongly affected by complex interaction mechanisms between multiple caverns. In the context of the deep geological disposal of nuclear waste, one main concern is that the disturbed zone around the excavations, namely the excavation damaged zone (EDZ), may negatively impact the hydro-mechanical behaviour of the rock mass, thus affecting its isolation properties and, as a consequence, the long-term safety of the repository.

Analytical solutions can be used to determine the stress and deformation fields around underground excavations (Brady and Brown, 2006). However, closed-form solutions are available only for simple excavation shapes (e.g. circular, elliptical) and under highly simplifying mechanical assumptions, such as perfect elasticity and homogeneity. Therefore, in engineering practice, numerical models are frequently

Corresponding author. Tel.: +1 416 799 9161. E-mail address: andrea.lisj ak@geomechanica. com

used to analyse and predict the rock mass behaviour. In computational geomechanics, the numerical approaches are commonly classified as (i) continuum methods and (ii) discontinuum (or discrete) methods (Jing and Hudson, 2002). Conventionally, numerical models based on continuum mechanics are employed to simulate rock mass response to excavation process (e.g. Mizukoshi and Mimaki, 1985; Eberhardt, 2001; Cai and Kaiser, 2014). However, their ability to consider the rock mass discontinuities remains somewhat limited. Although joint elements can be integrated into the continuum formulation (Hammah et al., 2008), typically only small deformations can be correctly captured due to the lack of contact detection and interaction algorithms (Cundall and Hart, 1992). On the other hand, discrete models may provide a more realistic representation of the physical behaviour observed in the field and, specifically, of the intrinsically discontinuous nature of the rock mass (Barton, 2011). Among the available discrete numerical approaches, the hybrid finite-discrete element method (FDEM) (Munjiza, 2004; Mahabadi, 2012) captures material failure by explicitly considering fracture nucleation and propagation, as well as the interaction of pre-existing and newly-created discrete rock blocks.

In this study, FDEM simulations are used to obtain unique insights into the failure process around deep underground excavations for three different geomechanical scenarios. Firstly, the stress-driven fracturing process of a circular shaft excavated in a homogeneous and isotropic rock is analysed. Secondly, the influence of mechanical anisotropy on the development of an EDZ around a tunnel in a layered rock formation is considered. Thirdly, the interaction mechanisms between two adjacent underground caverns are investigated, with particular emphasis on the rock mass response sensitivity to the in situ stress anisotropy, pillar

width and excavation sequence.

2. Review of available modelling approaches

Numerical modelling in rock engineering is a challenging task owing to several characteristics of the rock mass behaviour. Firstly, the stress-strain response of the rock material under uniaxial compression is highly non-linear (Martin, 1997; Jaeger et al., 2007). The initial strain hardening associated with the closure of voids and pre-existing microcracks is typically followed by a nearly linear stress-strain portion. Subsequently, the nucleation, propagation, and coalescence of microcracks lead to the loss of linearity, strain localization and the formation of macroscopic fractures. Upon reaching the peak strength, strain softening is associated with brittle rupture phenomena.

Secondly, the rock failure process is significantly influenced by confining pressure. Under unconfined compression (brittle) failure tends to occur in the form of axial splitting, while, under increasing confinement, the rock exhibits a more ductile behaviour accompanied by shear band formation.

Thirdly, the failure process observed at the laboratory-scale is further complicated at the rock mass level, where the behaviour is often influenced by the presence of discontinuities, such as joints, fractures, bedding planes, and tectonic structures. Discontinuities represent mechanical weaknesses of the rock mass and hence have a crucial effect on its deformability, strength, failure, and permeability properties (Hudson and Harrison, 1997). Moreover, the presence of discontinuities may add kinematic constraint on the deformation and failure mode of rock mass structures (Hoek et al., 1995). 2.1. Continuum approaches

The most commonly adopted numerical methods are the continuum-based approaches, such as the finite difference method (FDM), the finite element method (FEM), and the boundary element method (BEM). While FDM uses the differential form of the governing partial differential equations, FEM and BEM are based on their integral form and require solving a global equation system (Peiro and Sherwin, 2005). Continuum methods are suitable tools for simulating the stress and deformation fields around underground excavations. However, due to the lack of an internal length scale, standard strength-based, strain-softening constitutive relationships cannot reproduce the localisation of failure, as the underlying mathematical problem becomes ill-conditioned (de Borst et al., 1993).

To overcome the above limitations, different enrichment approaches, such as higher-order constitutive laws (e.g. Masin, 2005), Cosserat micro-polar models (e.g. Mühlhaus and Vardoulakis, 1987), non-local models (e.g. Bazant and Pijaudier-Cabot, 1988), and meshfree methods (e.g. Rabczuk and Belytschko, 2004, 2007; Zhuang et al., 2012, 2014) have been introduced. Recently, techniques such as the generalised finite element method (GFEM) (e.g. Strouboulis et al., 2000) and the extended finite element method (XFEM) (e.g. Möes and Belytschko, 2002), based on addition of non-polynomial shape functions to the classical FEM formulation, have been adopted for rock mechanics applications. Belytschko et al. (2001) used XFEM to investigate the stability of a tunnel in a jointed rock mass by modelling the fractures as interior displacement discontinuities. A similar approach was adopted by Deb and Das (2010) to numerically analyse a circular tunnel intersected by a joint plane. XFEM has also been successfully employed to simulate the propagation of cohesive cracks within continuum finite element models (Möes and Belytschko, 2002; Zhang and Feng, 2011). Unlike conventional fracture-mechanics-based studies (e.g. Steer et al.,

2011), in XFEM the discontinuities are completely independent of the finite element mesh and, therefore, remeshing is not required. However, the technique is, in general, not well suited to capture the interaction of multiple, arbitrarily located discontinuities, as well as large-scale material flow and motion (Karekal et al., 2011).

Another class of continuum-based approaches is represented by damage mechanics models, which capture the heterogeneous nature of rocks by statistically distributing defects into numerical domain. Several variations of this technique have been implemented in FEM (Tang and Kaiser, 1998), FDM (Fang and Harrison, 2002), smooth-particle hydrodynamics (SPH) (Ma et al., 2011), cellular automaton (CA) (Feng et al., 2006), and lattice (Blair and Cook, 1998) models. Among these implementations, the realistic failure process analysis (RFPA) code of Tang and Kaiser (1998) can provide an effective description of microscopic damage mechanisms by assuming a Weibull distribution of the mechanical parameters, including Young's modulus and strength properties (Zhu et al., 2005). Application of RFPA to simulate the evolution of the EDZ around a circular opening was illustrated by Zhu and Bruhns (2008) and Wang et al. (2009), in the presence of material anisotropy and under hydro-mechanically coupled conditions, respectively.

Rock mass discontinuities can be explicitly incorporated into continuum models by means of discrete joint (or interface) elements. This technique, originally proposed by Goodman et al. (1968) and known as the combined continuum-interface method (Riahi et al., 2010), is however suited only for a relatively low number of discontinuities. Alternatively, if the number of discontinuities is large and the discontinuities are not preferably oriented, homogenization techniques can be employed. That is, the rock mass is modelled as continuum with reduced deformation and strength properties accounting for the degrading effect of local geological conditions (Hoek et al., 2002; Hammah et al., 2008). Numerical homogenization of a continuum constitutive model can also be obtained from the results of discrete element simulations explicitly accounting for the presence of synthetic fracture networks (e.g. Beck et al., 2009). 2.2. Discontinuum approaches

In discrete (or discontinuous) modelling techniques, commonly known as the Discrete Element Method (DEM), the material is treated as an assembly of independent, rigid or deformable blocks or particles. Unique features of the DEM are the abilities to: (i) capture finite displacements and rotations of discrete bodies, including complete detachment, and (ii) automatically recognise new contacts as the simulation progresses (Cundall and Hart, 1992). Unlike continuum methods, which are based on constitutive laws, DEM relies on interaction laws. Based on the different solution strategies, DEMs can be divided into two main groups (Jing and Stephansson, 2007). The first group, usually referred to as the distinct element method, uses an explicit time-domain integration scheme with finite difference discretization to solve the equations of motion for rigid or deformable discrete bodies with deformable contacts (Cundall and Strack, 1979). The most widely used codes of this type are the universal distinct element code (UDEC) (Itasca, 2013) for blocky systems and the particle flow code (PFC) (Itasca, 2012) for granular systems. The second category uses an implicit (and thus unconditionally stable) time integration scheme and it is represented mainly by the discontinuous deformation analysis (DDA) method (Shi and Goodman, 1988). Further classification of DEMs is based on criteria such as the type of contact between bodies, the representation of deformability of solid bodies, and the methodology for detection and revision of contacts (Jing and

Stephansson, 2007).

While original applications of DEMs were mainly in the field of granular materials and jointed structures, further developments made DEMs also capable of explicitly simulating failure through intact rock material. Particularly, the concept of particle (or block) bonding, together with the introduction of cohesive contact models in DEMs, allowed the formation of new fractures to be captured. In this context, the FDEM (Munjiza, 2004; Mahabadi, 2012) adopted is a special type of discontinuum approach, whereby the simulation effectively starts with a continuous representation of the solid domain and, as the simulation progresses with time, new discontinuities are allowed to form upon satisfying some fracture criterion, thus leading to the formation of new discrete bodies. For a detailed review of discrete methods, and their application to underground structures, the reader is referred to Lisjak and Grasselli (2014).

3. Fundamental principles of FDEM

The modelling platform adopted for the numerical simulations was the open source FDEM software known as Y-Geo (Mahabadi et al., 2012). In Y-Geo, the modelling domain is discretised with a mesh consisting of three-node triangular elements with four-node interface (or crack) elements embedded between the edges of all adjacent triangle pairs (Fig. 1a). The progressive failure of rocks is simulated using a cohesive-zone approach, a technique originally introduced in the context of the elasto-plastic fracturing of ductile metals (Dugdale, 1960) and then extended to quasi-brittle materials, such as concrete and rocks (Hillerborg et al., 1976). During elastic loading, stresses and strains are assumed to be distributed over the bulk material (i.e. the continuum portion of the model), which is therefore treated as linear-elastic using the triangular elements. Impenetrability between these elements is enforced by a penalty-based contact interaction algorithm (Munjiza and Andrews, 2000). Upon exceeding the peak strength of the material (in tension, shear, or a mixed-mode), the strains are assumed to localise within a narrow zone, known as the Fracture Process Zone (FPZ). The mechanical response of the FPZ is captured by a non-linear interdependence between stress and crack displacement implemented at

three-noded, constant-strain elastic element

the crack element level.

The constitutive response of a crack element is defined in terms of a variation of the bonding stresses, s and t, between the edges of the triangular element pair as a function of the crack relative displacements, o and s, in the normal and tangential directions, respectively (Fig. 1b). In tension (i.e. Mode I), the response of each crack element depends on the cohesive tensile strength, ft, and the Mode I fracture energy, GIc. In shear (i.e. Mode II), the behaviour is governed by the peak shear strength, fs, and the Mode II fracture energy, GIIc. The peak shear strength is defined as

fs = c + sn tan ji (1)

where c is the cohesion, j is the internal friction angle, and sn is the normal stress acting across the crack element. The Mode I and Mode II fracture energies, Gic and Giic, represent the amount of energy, per unit crack length along the crack edge, consumed during the creation of a tensile and shear fracture, respectively. Upon breaking a crack element, a purely frictional resistance, fr, is assumed to act along the newly-created discontinuity:

fr = ontan jf (2)

where jf is the fracture friction angle. For mixed Mode I-II fracturing, an elliptical coupling relationship is adopted between crack opening, o, and slip, s (Fig. 1c). Although no deformation should, in theory, occur in the crack elements before the cohesive strength is exceeded, a finite cohesive stiffness is required by the formulation of FDEM. Such an artificial stiffness 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 cohesive contribution to the overall model compliance can be largely limited by adopting very high (i.e. dummy) penalty values (Munjiza, 2004; Mahabadi, 2012). Since fractures can nucleate only along the boundaries of the triangular elements, arbitrary fracture trajectories can be reproduced within the constraints imposed by the mesh topology. As the simulation progresses, through explicit time stepping, finite displacements and rotations of newly-created discrete bodies are allowed and new contacts are automatically recognised (Munjiza and Andrews, 1998).

Mode I

fs = c + crn tan cf>/

crack element yielding

crack element yielding

Mode I

four-noded o

crack element

pn, pf pf. penalty parameters fy c : cohesive strengths o, s : crack opening and sliding G|c, G||c : fracture energies h : element size <|>/, : friction angles

Fig. 1. Simulation of rock deformation and fracturing with FDEM. (a) Representation of a continuum using cohesive crack elements interspersed throughout a mesh of triangular elastic elements. Triangles are shrunk for illustration purposes. (b) Constitutive behaviour of the crack elements defined in terms of normal and tangential bonding stresses, cr and t, versus crack relative displacements, o and s (i.e. opening and slip). (c) Elliptical coupling relationship between o and s, for mixed-mode fracturing.

The FDEM formulation described above was originally introduced to model isotropic materials. However, additional capabilities have been recently introduced into the FDEM solver to capture the mechanical response of anisotropic media (Lisjak et al., 2014a, b). In particular, the modulus anisotropy is captured by a transversely isotropic, linear elastic constitutive law implemented at the triangular element level. In this case, the elastic deformation is fully characterised by five independent elastic parameters: two Young's moduli, Ep and Es, and Poisson's ratios, VP and vs, for the directions parallel and perpendicular to the plane of isotropy, and the shear

modulus, GS. The anisotropy of strength is instead introduced at the crack element level by specifying the cohesive strength of each crack element as a function of the relative orientation, g between the crack element itself and the bedding orientation (Fig. 2a). The cohesive strength parameters and the fracture energies are assumed to vary linearly between a minimum value for g=0° (i.e. ft,min, Cmin, Gfc,min, Gnc,min) to a maximum value for g=90° (i.e. ft,max, cmax, GIc,max, GIIc,max). Furthermore, the mesh topology combines a random triangulation for the intra-layer material (i.e. matrix) together with crack elements preferably aligned along the plane of isotropy (Fig. 2b).

-t Min -

Angle between Crack Element and Bedding Plane, y(°]

Fig. 2. FDEM modelling of strength anisotropy. (a) Linear variation of cohesive strength parameters with the angle, g between crack element and layering orientation. (b) Example of mesh combining a Delaunay triangulation for the intra-layer material with edges preferentially aligned along the isotropy direction (after Lisjak et al. (2014a)).

4. Case studies

4.1. Model description

4.1.1. Geometries

Three different geomechanical scenarios were considered in the FDEM simulations: (i) a 6-m-diameter vertical shaft sunk in a horizontally bedded rock formation (shaft model, Fig. 3a), (ii) a 3-m-diameter tunnel excavated in a bedded formation (tunnel model, Fig. 3b), and (iii) two adjacent horseshoe-shaped caverns excavated in an isotropic rock (cavern model, Fig. 3c). The openings were placed at the centre of a square domain with dimensions equal to 100 mx100 m, 50 mx50 m, and 500 mx500 m, for the shaft, tunnel, and cavern model, respectively. To maximise the model resolution in the EDZ, while keeping the run times within practical limits, a mesh refinement zone was adopted around the excavation boundaries, with an average element size of 0.06 m, 0.03 m, and 0.75 m, for the shaft, tunnel, and cavern model, respectively. The sensitivity of the model to variations in element size and topology was not investigated. In the tunnel model, the cross-section was assumed to be perpendicular to the strike of bedding planes inclined at y=33° from the horizontal. In the cavern model, five and four sub-domains were adopted for the powerhouse and transformer caverns, respectively, to analyse the effect of excavation staging on the fracturing process.

4.1.2. In situ stresses and boundary conditions

To correctly simulate the prior-to-excavation stress state, each model required two separate runs. In the first run, the vertical and horizontal in situ stress conditions, as reported in Table 1, were applied to the model without the insertion of crack elements. Gravity-induced stress gradients were neglected. As suggested by Hudson and Harrison (1997), in the cavern model, three different in situ stress fields were simulated (stress ratio K0=0.5, K0=1.0 and K0=1.33) in order to investigate the cases of pillar over-stressing (for K0<1) as well as stress shadowing (for K0>1). The first run was continued until the total kinetic energy of the system decayed to a negligible value (i.e. resulting stress waves were attenuated). The revised nodal coordinates corresponding to the system at rest (i.e. static equilibrium) were then obtained. Subsequently, these revised nodal coordinates were used as the current nodal coordinates (i.e. deformed mesh) of the second run in which the actual material strengths were assigned. By changing the far-field boundaries to be fixed in the horizontal and vertical directions, the first order in situ conditions were maintained while allowing the excavation to be initiated. Model relaxation induced by the artificial compliance of the crack elements was minimised by the choice of sufficiently large contact and fracture penalty values. It is noteworthy that only in-plane stresses were effectively used in the analysis, as the crack element formulation cannot account for the influence of an out-of-plane stress. 4.1.3. Excavation and support modelling

With the correct in situ stress conditions achieved, the openings were created using a core replacement technique. With this approach, the

three-dimensional supporting effect of the excavation face, which causes a gradual reduction of radial resistance around the excavation boundary, is captured by a fictitious, softening elastic material placed in the excavation core. In general, with this method, the deformation modulus of the excavated material is progressively reduced from the original rock mass value, corresponding to an undeformed section far ahead of the face, to a value that results in the wall displacements at the time of support installation. In this work, no attempt was made to match any real deformation measurements and, therefore, the modulus reduction ratio at the time of support activation was arbitrarily chosen. Since the deformation modulus of the excavation core was reduced over time in a stepwise fashion, the total kinetic energy of the model was again monitored to ensure that steady-state conditions were reached at every excavation stage. The final stage of the excavation sequence also involved the actual material removal and, for the shaft and tunnel models only, the activation of the support layer. To simplify the analysis, the effect of rock support was not considered in the cavern model. The application of shotcrete on the tunnel walls was modelled using constant-strain linear-elastic triangular elements. The support installation consisted of specifying the liner thickness and the installation time from a given core softening ratio. Since the delayed installation of shotcrete was accomplished by varying the elastic properties of the liner (from those of the rock mass to those of the shotcrete), the deformation in the liner had to be zeroed to avoid an artificial build-up of stress in response to an instantaneous increase of material stiffness in a pre-stressed medium. 4.1.4. Input parameters

Since the shaft was mined perpendicular to the layering strike, an

isotropic mechanical model was assumed with input parameters based on laboratory values for an indurated claystone from Northern Switzerland (unpublished report) (Table 2). For the tunnel model, the rock mass was modelled using an anisotropic strength and stiffness model with a layering thickness of 0.1 m (Fig. 2). The input elastic properties as well as the cohesive strength parameters were those obtained from the back-analysis of a test tunnel excavated in an anisotropic shale formation (Opalinus Clay) at the Mont Terri underground research laboratory (URL) (Switzerland) (see Lisjak (2013) and Lisjak et al. (2014c) for further details). The rock mass parameters for the cavern model were based on unpublished laboratory values of a gneissic rock.

4.2. Fracturing process around a circular shaft

The simulation results of the shaft model highlight the stress-driven nature of the rock mass failure process under homogeneous and isotropic conditions. Upon reducing the elastic modulus of the tunnel core, the N-S-oriented in situ maximum principal stress flows around the shaft boundary, resulting in the development of a compressive stress concentration in the sidewalls. The intensity of this stress concentration is such that shear (i.e. Mode II) fractures start to nucleate (Fig. 4a). Due to the stress-free surface created by the excavation process, a state of unconfined (or moderately confined) compression arises in proximity to the shaft walls. Consequently, the failure mode closely resembles that observed for rock specimens subjected to uniaxial compressive stress. In agreement with the Mohr-Coulomb failure criterion, conjugate shear cracks tend to develop at angle of 45°± j /2 to the vertical compressive stress. As the shaft face advances

Fig. 3. Geometry and boundary conditions of the FDEM models: (a) shaft, (b) tunnel, and (c) cavern.

Table 1. Summary of in situ stress conditions applied to the excavation models. Vertical and horizontal stresses are oriented along the principal directions. The stress ratio, Ko, is reported in brackets.

Model Vertical stress, Cv (MPa) Maximum horizontal stress, cth (MPa) Minimum horizontal stress, ah (MPa)

Shaft 19.6 15.7 (0.8)

Tunnel 6.5 4.5 (0.7)

8 4 (0.5)

Cavern 6 6 (1.0)

6 8 (1.3)

(i.e. the core modulus is further reduced), the shear fractures tend to propagate away from the excavation and, at the same time, tend to curve and realign themselves in the direction of the far-field maximum principal stress (Fig. 4b, c). As a result, a characteristic fracture pattern consisting of multiple families of cracks resembling logarithmic-spiral rupture surfaces is created. The mutual intersection of these slip lines tends to break up the rock mass by forming distinct blocks and fragments. In close vicinity to the excavation walls, the occurrence of

rock crushing and fine fragmentation is due to the higher stress concentration. At a distance from the excavation boundary, shearing of these fractures causes a local stress redistribution which tends to protect the intact rock. The activation of the shotcrete layer stabilises the rock fracturing process, thus allowing static equilibrium conditions to be reached (Fig. 4d). The final EDZ assumes an elliptical shape with major and minor axes oriented parallel to the minimum and maximum in situ principal stress directions, respectively (Fig. 5). The fractured zone extends for roughly 10 m (i.e. 1.7 times of shaft diameter) and 5 m (i.e. 0.8 times of shaft diameter) in the E-W and N-S directions, respectively. Arching of the maximum compressive stress, S, can be observed around the damaged area, while the confining stress, S3, decreases to zero inside the EDZ. Overall, the simulated fracturing process is in good qualitative agreement with the mechanisms described by Barton (1993) and with experimental observations of borehole stability in massive isotropic rocks (e.g. Addis et al., 1990) and of excavation-induced fracture networks in a claystone formation (Armand et al., 2014).

Table 2. FDEM input parameters of the three excavation models.

Input parameters

Continuum triangular elements Formulation type Bulk density, p (kg/m3) Young's modulus, E (GPa) Young's modulus parallel to bedding, EP (GPa) Young's modulus perpendicular to bedding, ES (GPa) Poisson's ratio, v

Poisson's ratio parallel to bedding, vP Poisson's ratio perpendicular to bedding, Vs Shear modulus, GS (GPa) Viscous damping coefficient, m (kg/(m s)) Crack elements Formulation type Tensile strength, f (MPa) Tensile strength parallel to bedding, ft,max (MPa) Tensile strength perpendicular to bedding, ft,min (MPa) Cohesion, c (MPa)

Cohesion parallel to bedding, cmin (MPa)

Cohesion perpendicular to bedding, cmax (MPa)

Mode I fracture energy, Gic (J/m2)

Mode I fracture energy parallel to bedding, GIc,max (J/m2)

Mode I fracture energy perpendicular to bedding, GIc,min (J/m2)

Mode II fracture energy, GIIc (J/m2)

Mode II fracture energy parallel to bedding, GIIc,min (J/m2)

Mode II fracture energy perpendicular to bedding, GIIc,max (J/m2)

Friction angle of intact material, j (°)

Friction angle of fractures, j (°)

Normal contact penalty, pn (GPa m)

Tangential contact penalty, pt (GPa/m)

Fracture penalty, pf (GPa)

Isotropic 2430 1 11.4

Isotropic 1.5

24 24 114 11.4

Anisotropic 2430

3.8 1.3

0.35 0.25 3.6

1.83X105

Anisotropic

1.8 0.44

19.5 1

Isotropic

Isotropic 8

1.25x10

a = 0.5, time step = 40k

as= 0.1, time step = 70k

as= 0.01, time step = 150k

equilibrium, time step = 490k

19.6 MPa |

15.7 MPa

Mode of Fracture j 3 m j — tensile-dominated - shear-dominated

Fig. 4. Shaft model: simulated evolution of fracture growth around the opening at increasing simulation times corresponding to different stages of the core modulus reduction sequence. The core modulus reduction ratio, a, is equal to the ratio of the core modulus to the rock mass modulus.

Maximum Principal Stress, a.

19.6 MPa I 15.7 MPa

Minimum Principal Stress, a3

Fig. 5. Shaft model: contours of maximum and minimum principal stresses associated with the excavation at equilibrium.

4.3. Influence of mechanical anisotropy on a circular tunnel

The numerical results of the tunnel model indicate that failure around an opening in a layered formation is triggered by the excavation-induced stress redistribution in combination with the lower strength of bedding planes favourably oriented for slip. As depicted in Fig. 6, fractures start to develop around the excavation boundary at approximately 0°<0<15o, 120°<6£195°, and 300°<6£360° in the form of shear-dominated (i.e. Mode II) fractures along the bedding direction. The polar orientation of these slip zones corresponds to critical values of relative orientation between compressive stress around the excavation boundary and bedding favourably oriented for slippage. As the simulation progresses, the slippage of bedding planes causes a local perturbation in the stress field which results in the nucleation of strain-driven, Mode I fractures in the direction perpendicular to the layering (Fig. 6b). Also, bedding-parallel sliding is simulated at about 70° and 250°.

Further rock mass deconfinement triggers further delamination of bedding planes (Fig. 6c) and the formation of wing-shaped fractured zones that tend to extend out in the direction parallel to the bedding to a distance of about 3 m from the sidewalls. After installing the support, the propagation of damage away from the opening is suppressed in favour of fragmentation in close proximity to the excavation boundary, until new equilibrium conditions are reached (Fig. 6d).

Stress, a (MPa)

0 10 20 30

a = 0.5, time step = 30k

as= 0.1, time step = 70k

as= 0.01, time step = 110k

equilibrium, time step = 250k al

6.5 MPa | 4.5 MPa

Mode of Fracture

I 1 m j - tensile-dominated

— shear-dominated

Fig. 6. Tunnel model: simulated evolution of fracture growth around the tunnel at increasing simulation times corresponding to different stages of the core modulus reduction sequence.

The total displacement field associated with the tunnel configuration at equilibrium (Fig. 7) indicates that at a distance from the excavation, the rock mass behaves elastically and, therefore, small strains, induced by the stress redistribution around the damaged zone, are simulated. Due to the highly anisotropic rock mass response, this distance varies from a minimum of 0.5 m to a maximum of 3 m in the direction parallel to bedding and in the sidewalls, respectively. Furthermore, elastic deformations of higher intensity are captured in the direction sub-perpendicular to the bedding orientation due to the high rock compressibility in the said direction. In the near-field excavation, an inner and an outer shell can be identified. The shape of the inner zone is roughly a 4.5 mx4.5 m square with edges oriented in the direction parallel and perpendicular to the bedding and centre coincident with the tunnel axis. In this zone, the rock mass deformation is governed by a combination of Mode I and Mode II fracturing and bulking, thus resulting in large displacements (i.e. d>3 cm). In the outer shell, while Mode II fractures can still nucleate, the relative sliding along the fracture surfaces is limited by higher values of confining stress, 03. Consequently, the growth of extensional fractures is effectively inhibited.

The redistribution of compressive stress in response to the tunnel excavation (Fig. 8) is influenced by the in situ stress anisotropy as well as the characteristic fracture pattern (bedding-parallel discontinuities and a heavily fractured zone around the tunnel). The lateral extension of the EDZ due to bedding delamination is suppressed by the re-orientation of 01 in the direction perpendicular to bedding. In proximity to the tunnel boundary, bedding plane slippage promotes a drastic reduction of confining stress, S3, with low to moderate negative values responsible for the observed extensional fracturing.

Fig. 7. Tunnel model: colour contour of total displacement associated with the tunnel configuration at equilibrium.

Maximum Principal Stress, g1

Minimum Principal Stress, cj3

Stress, a (MPa)

-2 0 5 10 15

Fig. 8. Tunnel model: colours contour of principal stresses associated with the tunnel configuration at equilibrium. Principal stress directions are indicated by short straight lines.

A quantitative comparison of the simulated damage pattern with specific in situ observations is beyond the scope of this work. Nevertheless, the failure mechanisms simulated here are in general agreement with a number of field and laboratory observations of excavations in laminated rock formations. In particular, the

characteristic shear failure of bedding planes was observed in Opalinus Clay during hollow cylinder experiments (Labiouse and Vietor, 2014), and around boreholes, microtunnels and drifts at the Mont Terri URL (Marschall et al., 2006; Blumling et al., 2007). Also, the importance of weakness planes in controlling the rock mass behaviour and the stability of underground openings is confirmed by observations from the construction of a hydroelectric tunnel in laminated sedimentary formations (Perras and Diederichs, 2009). Characteristic square-shaped fractured zones have also been reported in the hydrocarbon exploration industry when drilling horizontal boreholes in laminated shales (e.g. 0kland and Cook, 1998; Willson et al., 1999). 4.4. Interaction mechanisms between two adjacent caverns 4.4.1. Effect of in situ stress

The numerical results indicate a critical influence of the in situ stresses on the fracture development around the two caverns (Fig. 9). In general, for the adopted rock mass properties and cavern configuration, an isotropic stress field induces the lowest deviatoric stresses in the surrounding rock mass and, therefore, minimises damage development in the pillar. On the other hand, a vertically oriented maximum in situ principal stress leads to pillar over-stressing and failure. Also, higher stress concentrations tend to occur, in agreement with the analysis of Brady and Brown (2006), around high-curvature boundaries, which therefore become preferential loci of fracture initiation.

For the case of K0=0.5 (Fig. 9a), the stress channelling within the rock pillar causes the development of a through-going macroscopic fracture plane, resembling that often observed in rock specimens subjected to uniaxial compression. The EDZ starts to form with a fracture growing from the lower right corner of the powerhouse cavern to the upper left corner of the transformer chamber. The stress redistribution causes further fracturing to originate from the centre of the pillar and propagate towards the upper right sidewall of the powerhouse and the lower left sidewall of the transformer cavern. Although a low to slightly negative confining stress, 03, develops in the rock pillar, the fracturing process is dominated by Mode II failure, due to the relatively high tensile strength of the rock compared to the cohesion value (Table 1).

Under isotropic stress conditions (K0=1.0, Fig. 9b), damage develops around the lower corners of the two excavations and above the arched roof of the powerhouse cavern. Unlike the previous case (K0=0.5), the rock pillar remains substantially intact. For the case of K=1.33 (Fig. 9c), the simulated failure process leads to an EDZ network that is overall similar to the isotropic case. One notable difference is represented by the larger fractured areas simulated in the back and roof of the powerhouse cavern due to the less favourable orientation of the in situ stress field.

Stress, a (MPa) Mode of Fracture

1 1 — tensile-dominated

"2 0 5 10 15 ,, ■ + ^

- shear-dominated

Fig. 9. Effect of in situ stress anisotropy in the cavern model. Final stress distribution and fracture pattern for the cases of (a) Ko=0.5, (b) Ko=1.0, and (c) Ko=1.33. Colour contours on the left and right hand side represent the maximum and minimum principal stresses, si and s, respectively. Local principal stress directions are indicated by short straight lines.

4.4.2. Effect of pillar width

As described by Hoek (2006), the distance between the two caverns should be as small as possible to minimise the length of the busbars that connect the generators in the powerhouse cavern to the transformers in the adjacent cavern. On the other hand, this distance has to be large enough to preserve the structural integrity of the pillar. Therefore, the optimisation of the pillar width represents a crucial aspect of the design

of this type of underground structures. In this study, the pillar damage was simulated for cavern spacing values, s, of 25.5 m, 35.5 m and 45.5 m, corresponding to ratios of the cavern width to the pillar width equal to 1.2, 1.7 and 2.2, respectively. Only the case of K0=0.5 was analysed. As expected, the extent of the damaged area decreases as the spacing increases (Fig. 10). For s=25.5 m (Fig. 10a), the EDZ spans the entire width of the pillar and is characterised by heavy fragmentation

developing between the two caverns (see also Section 4.4.1). For s=35.5 m, the damaged zones of the two caverns are still inter-connected (Fig. 10b), however fracturing is sensibly less intense than the previous case and a narrow load-bearing zone is preserved at the centre of the pillar. Also, unlike the first case, the stress concentration along the external sidewalls is too low to induce significant fracturing. Lastly, for s=45.5 m (Fig. 10c), although a noticeable disturbance of the pre-excavation stress field is still presented, the stress redistribution does not cause a connected failure pattern to develop. Instead, an asymmetric EDZ fracture network is created: fracturing is concentrated along the inner sidewall of the transformer cavern, while, interestingly, the rock mass around the powerhouse cavern remains nearly intact.

Fig. 10. Effect of pillar width in the cavern model. Final stress distribution and fracture pattern are presented for a cavern spacing of (a) 25.5 m, (b) 35.5 m, and (c) 45.5 m. Colour contours represent the maximum principal stresses, s1. Local major principal stress directions are indicated by short straight lines.

4.4.3. Effect of excavation staging

The excavation of the entire cavern cross-section at once arguably represents an extreme case, which leads, as described above for K)=0.5, to a large over-stressed and fractured area within the rock pillar. To investigate the adoption of a more realistic excavation procedure, further simulations were carried out, for the case of K0=0.5 and s=35.5 m, whereby both caverns were excavated using a multi-stage sequence. The excavation process started with the top heading of the powerhouse cavern (stage I), followed by a sequential excavation of the remaining four sub-domains of both caverns (stage II-V). The simulation results indicate that the stress history plays an important role in controlling the evolution of failure and the final fracture pattern (Fig. 11). For the one-stage excavation (Fig. 10b), fractures initiate from the bottom right corner of the powerhouse cavern and propagate towards the back of the transformer cavern. In contrast, in the case of a staged excavation, the failure process starts from the upper right sidewall of the powerhouse cavern (Fig. 11a). The material removal is accompanied by the downward growth of shear fractures sub-parallel to the caverns' inner sidewalls (Fig. 11b, c). Compared to the single-stage simulation, the final fracture pattern of the multi-stage model (Fig. 11d) shows (i) a larger residual intact area in the centre of the pillar and (ii) more fragmentation in close proximity to the sidewalls, due to a repeated occurrence of high compressive stresses around the corners of the excavation benches.

Stress, q (MPa)_ Mode of Fracture

- tensile-dominated

O 1U 1 o

- shear-dominated

20 m I-1

8 MPaj

, 'ús ^ЛЬ-^Л i flT'1 '< b ¡ ь il'/ к%|('V1

il II I I ti l ' . ^

8 MPa 4 MPa

Stress, a (MPa)

Mode of Fracture

-2 0 5 10 15 — tensile-dominated

shear-dominated

Fig. 11. Effect of excavation staging in the cavern model. Failure sequence and maximum principal stress contour within the rock pillar for the case of K0=0.5 and spacing of 35.5 m after (a) stage II, (b) stage III, (c) stage IV, and (d) stage V. Local major principal stress directions are indicated by short straight lines.

5. Concluding remarks

Hybrid continuum-discontinuum simulations, based on the FDEM, were used to investigate excavation-induced fracturing processes around different types of underground structures. Three main geomechanical scenarios were considered.

Firstly, for a circular shaft excavated in a homogenous and isotropic medium, shear failure started in the regions of the highest excavation-induced compressive stress concentration. Subsequently, fractures tended to follow characteristic trajectories, resembling the logarithmic-spiral slip zones captured by conventional Mohr-Coulomb elasto-plastic models, and in agreement with experimental observations from borehole breakout experiments.

Secondly, in the case of a tunnel excavated in a laminated shale, the bedding induced mechanical anisotropy was shown to sensibly influence the locus of fracture initiation as well as the direction of fracture growth. Due to the lower shear strength along the bedding, the rock failure was strongly dependent upon the relative orientation between bedding planes and in situ principal stress directions. Furthermore, the initial shearing along bedding planes induced a tensile stress state in the perpendicular direction with consequent formation of secondary strain-driven extensional fractures.

The third scenario focused on the behaviour of two adjacent

horseshoe-shaped caverns. For the case of a vertically oriented in situ maximum principal stress, pillar over-stressing with formation of a through-going shear fracture plane was simulated, whereas isotropic in situ stress conditions resulted in the least amount of rock damage. A sensitivity analysis to the pillar width revealed that, to avoid the formation of an interconnected EDZ between the two caverns and, therefore, preserve the pillar load-bearing capacity; the cavern spacing should be greater than about two times the cavern width. Finally, the adoption of a multi-stage excavation sequence was shown to affect the fracture growth as well as the final damage pattern.

In conclusion, the modelling results indicate that FDEM simulations can provide unique geomechanical insights in all those cases where an explicit consideration of fracture and fragmentation processes is of paramount importance.

Conflict of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgements

This work has been supported by the Natural Science and

Engineering Research Council (NSERC) of Canada in the form of discovery grant No. 341275 and the Swiss National Cooperative for the Disposal of Radioactive Waste (NAGRA).

References

Addis MA, Barton NR, Bandis SC, Henry JP. Laboratory studies on the stability of vertical and deviated boreholes. In: Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, September 23-26, 1990.

Armand G, Leveau F, Nussbaum C, de La Vaissiere R, Noiret A, Jaeggi D, Landrein P, Righini C. Geometry and properties of the excavation-induced fractures at the Meuse/Haute-Marne URL drifts. Rock Mechanics and Rock Engineering 2014; 47(1): 21-41.

Barton NR. From empiricism, through theory, to problem solving in rock engineering. In: Harmonising Rock Engineering and the Environment. London, UK: CRC Press; 2011. pp. 3-16.

Barton NR. Physical and discrete element models of excavation and failure in jointed rock. In: Keynote Lecture Presented at ISRM International Symposium on Assessment and Prevention of Failure Phenomena in Rock Engineering, Istanbul, Turkey. 1993.

Bazant ZP, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. Journal of Applied Mechanics 1988; 55(2): 287-93.

Beck DA, Pfitzner MJ, Arndt SM, Fillery B. Estimating rock mass properties and seismic response using higher order, discontinuous, finite element models. In: Proceedings of the 3rd Canada-US Rock Mechanics Symposium. Toronto, Canada: University of Toronto Press; 2009. Paper 4189.

Belytschko T, Moes N, Usui S, Parimik C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering 2001; 50(4): 993-1013.

Blair S, Cook NGW. Analysis of compressive fracture in rock using statistical techniques: Part I. A non-linear rule-based model. International Journal of Rock Mechanics and Mining Sciences 1998; 35(7): 837-48.

Blümling P, Bernier F, Lebon P, Martin CD. The excavation damaged zone in clay formations time-dependent behaviour and influence on performance assessment. Physics and Chemistry of the Earth, Parts A/B/C 2007; 32(8-14): 588-99.

Brady BHG, Brown ET. Rock mechanics for underground mining. Dordrecht, Netherlands: Springer, 2006.

Cai M, Kaiser PK. In situ rock spalling strength near excavation boundaries. Rock Mechanics and Rock Engineering 2014; 47(2): 659-75.

Cundall PA, Hart RD. Numerical modelling of discontinua. Engineering Computations 1992; 9(2): 101-13.

Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979; 29(1): 47-65.

de Borst R, Sluis LJ, Mühlhaus HB, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Engineering Computations 1993; 10(2): 99-122.

Deb D, Das KC. Extended finite element method for the analysis of discontinuities in rock masses. Geotechnical and Geological Engineering 2010; 28(5): 643-59.

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

Eberhardt E. Numerical modelling of three-dimension stress rotation ahead of an advancing tunnel face. International Journal of Rock Mechanics and Mining Sciences 2001; 38(4): 499-518.

Fang Z, Harrison JP. Development of a local degradation approach to the modelling of brittle fracture in heterogeneous rocks. International Journal of Rock Mechanics and Mining Sciences 2002; 39(4): 443-57.

Feng XT, Pan PZ, Zhou H. Simulation of the rock microfracturing process under uniaxial compression using an elasto-plastic cellular automaton. International Journal of Rock Mechanics and Mining Sciences 2006; 43(7): 1091-108.

Goodman RE, Taylor RL, Brekke TA. A model for the mechanics of jointed rock. Journal of the Soil Mechanics and Foundations Division 1968; 94(3): 637-660.

Hammah RE, Yacoub T, Corkum B, Curran JH. The practical modelling of discontinuous rock masses with finite element analysis. In: Proceedings of the 42nd US Rock Mechanics Symposium and 2nd US-Canada Rock Mechanics Symposium, San Francisco, USA. Alexandria, Virginia, USA: American Rock Mechanics Association; 2008. ARMA 08-180.

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

Hoek E, Carranza-Torres CT, Corkum B. Hoek-Brown failure criterion - 2002 edition. In: Proceedings of the 5th North American Rock Mechanics Symposium. Toronto, Canada: University of Toronto Press; 2002. pp. 267-73.

Hoek E, Kaiser PK, Bawden WF. Support of underground excavations in hard rock. Leiden, Netherlands: Taylor & Francis/Balkema, 1995.

Hoek E. Practical rock engineering. 2006. http://www.rocscience.com/ hoek/pdf/Practical_Rock_Engineering.pdf.

Hudson JA, Harrison JP. Engineering rock mechanics, volume I. Pergamon, 1997.

Itasca. PFC2D (particle flow code in 2 dimensions). Minneapolis, MN, USA: Itasca Consulting Group Inc., 2012.

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

Jaeger JC, Cook NGW, Zimmerman RW. Fundamentals of rock mechanics, 4th ed. Malden, MA, USA: Blackwell Publishing, 2007.

Jing L, Hudson JA. Numerical methods in rock mechanics. International Journal of Rock Mechanics and Mining Sciences 2002; 39(4): 409-27.

Jing L, Stephansson O. Fundamentals of discrete element methods for rock engineering: Theory and applications, volume 85 (developments in geotechnical engineering). Amsterdam, Netherlands: Elsevier, 2007.

Karekal S, Das R, Mosse L, Cleary PW. Application of a mesh free continuum method for simulation of rock caving processes. International Journal of Rock Mechanics and Mining Sciences 2011; 48(5): 703-11.

Labiouse V, Vietor T. Laboratory and in situ simulation tests of the excavation damaged zone around galleries in Opalinus Clay. Rock Mechanics and Rock Engineering 2014; 47(1): 57-70.

Lisjak A, Grasselli G. A review of discrete modelling techniques for fracturing processes in discontinuous rock masses. Journal of Rock Mechanics and Geotechnical Engineering 2014; 6(4): 301-14.

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 2014a; 65: 96-115.

Lisjak A, Tatone BSA, 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 2014b; 47(1): 187-206.

Lisjak A, Garitte B, Grasselli G, Müller H, Vietor T. The excavation of a circular tunnel in a bedded argillaceous rock (Opalinus Clay): Short-term rock mass response and numerical analysis using FDEM. Tunnelling and Underground Space Technology 2014c; doi:10.1016/j.tust.2014.09.014 (in press).

Lisjak A. Investigating the influence 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. http://hdl.handle.net/1807/43649.

Ma GW, Wang XJ, Ren F. Numerical simulation of compressive failure of heterogeneous rock-like materials using SPH method. International Journal of Rock Mechanics and Mining Sciences 2011; 48(3): 353-63.

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

Mahabadi OK. Investigating the influence of micro-scale heterogeneity and

microstructure on the failure and mechanical behaviour of geomaterials. Ph.D. thesis. Toronto, Canada: University of Toronto, 2012. Http://hdl. handle.net/1807/32789.

Marschall P, Distinguin M, Shao H, Bossart P, Enachescu C, Trick T. Creation and evolution of damage zones around a microtunnel in a claystone formation of the Swiss Jura Mountains. In: Proceedings of the International Symposium and Exhibition on Formation Damage Control. Lafayette, Louisiana, USA: Society of Petroleum Engineers; 2006.

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

Masin D. A hypoplastic constitutive model for clays. International Journal for Analytical and Numerical Methods in Geomechanics 2005; 29(4): 311-36.

Mizukoshi T, Mimaki Y. Deformation behaviour of a large underground cavern. Rock Mechanics and Rock Engineering 1985; 18(4): 227-51.

Moes N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 2002; 69(7): 773-81.

Muhlhaus HB, Vardoulakis I. The thickness of shear bands in granular materials. Geotechnique 1987; 37(3): 845-57.

Munjiza A, Andrews KRF. NBS contact detection algorithm for bodies of similar size. International Journal for Numerical Methods in Engineering 1998; 43(1): 131-49.

Munjiza A, Andrews KRF. Discretised penalty function method in combined finite-discrete element analysis. International Journal for Numerical Methods in Engineering 2000; 49(12): 1495-520.

Munjiza A. The combined finite-discrete element method. Chichester, West Sussex, England: John Wiley & Sons Ltd., 2004.

0kland D, Cook JM. Bedding-related borehole instability in high angle wells. In: Proceedings of the SPE/ISRM Rock Mechanics in Petroleum Engineering, Society of Petroleum Engineers, Trondheim, Norway. 1998.

Peiro J, Sherwin S. Finite difference, finite element and finite volume methods for partial differential equations, chapter 8.2. Dordrecht, Netherlands: Springer, 2005.

Perras MA, Diederichs MS. Tunnelling in horizontally laminated ground. In: Diederichs MS, Grasselli G. (Eds.), Proceedings of 3rd Canada-US (CANUS) Rock Mechanics Symposium (RockEng09), Toronto, Canada. 2009.

Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering 2007; 196(29-30): 2777-99.

Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering 2004;

61(13): 2316-43.

Riahi A, Hammah ER, Curran JH. Limits of applicability of the finite element explicit joint model in the analysis of jointed rock problems. In: Proceedings of the 44th American Rock Mechanics Association (ARMA) Symposium. 2010.

Shi G, Goodman RE. Discontinuous deformation analysis - A new method for computing stress, strain and sliding of block systems. In: Proceedings of the 29th US Symposium on Rock Mechanics. Minneapolis, USA: A.A. Balkema; 1988. pp. 381-93.

Steer P, Cattina R, Lave J, Godardd V. Surface Lagrangian remeshing: A new tool for studying long term evolution of continental lithosphere from 2D numerical modelling. Computers and Geosciences 2011; 37(8): 1067-74.

Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized finite element method. Computer Methods in Applied Mechanics and Engineering 2000; 181(1-3): 43-69.

Tang CA, Kaiser PK. Numerical simulation of cumulative damage and seismic energy release during brittle rock failure - Part I: fundamentals. International Journal of Rock Mechanics and Mining Sciences 1998; 35(2): 113-21.

Wang SH, Lee CI, Ranjith PG, Tang CA. Modelling the effects of heterogeneity and anisotropy on the excavation damaged/disturbed zone (EDZ). Rock Mechanics and Rock Engineering 2009; 42(2): 229-58.

Willson SM, Last NC, Zoback MD, Moos D. Drilling in South America: A wellbore stability approach for complex geologic conditions. In: Proceedings of the Latin American and Caribbean Petroleum Engineering Conference, Society of Petroleum Engineers, Caracas, Venezuela. 1999.

Zhang YL, Feng XT. Extended finite element simulation of crack propagation in fractured rock masses. Materials Research Innovations 2011; 15: 594-6.

Zhu WC, Bruhns OT. Simulating excavation damaged zone around a circular opening under hydromechanical conditions. International Journal of Rock Mechanics and Mining Sciences 2008; 45(5): 815-30.

Zhu WC, Liu J, Tang CA, Zhao XD, Brady BH. Simulation of progressive fracturing processes around underground excavations under biaxial compression. Tunnelling and Underground Space Technology 2005; 20(3): 231-47.

Zhuang X, Augarde CE, Mathisen KM. Fracture modelling using meshless methods and level sets in 3D: Framework and modelling. International Journal for Numerical Methods in Engineering 2012; 92(11): 969-98.

Zhuang X, Zhu H, Augarde C. An improved meshless Shephard and least square method possessing the delta property and requiring no singular weight function. Computational Mechanics 2014; 53(2): 343-57.

Andrea Lisjak Bradley obtained a BSc (2006) and MSc (2008) in civil engineering from the University of Trieste, Italy, and a PhD (2013) in rock mechanics from the University of Toronto, Canada. Currently, he is a geomechanics specialist and NSERC R&D Fellow with Geomechanica Inc. in Toronto, Canada. Andrea's areas of expertise lie in the development and use of advanced hybrid continuum-discontinuum numerical methods to investigate failure processes in rocks. His research interests focus on the numerical simulation of fracture development around underground excavations, rock support, and hydro-mechanical behaviour of discontinuous rock masses.