Hybrid FiniteDiscrete Element Simulation of the E

Authors: Andrea Lisjak Bryan S A Tatone Omid K Mahabadi Giovanni Grasselli Paul Marschall George W Lanyon Rémi de la Vaissière Hua Shao Helen Leung Christophe Nussbaum

Publish Date: 2015/09/16

Volume: 49, Issue: 5, Pages: 1849-1873

PDF Link

Abstract

The analysis and prediction of the rock mass disturbance around underground excavations are critical components of the performance and safety assessment of deep geological repositories for nuclear waste In the short term an excavation damaged zone EDZ tends to develop due to the redistribution of stresses around the underground openings The EDZ is associated with an increase in hydraulic conductivity of several orders of magnitude In argillaceous rocks sealing mechanisms ultimately lead to a partial reduction in the effective hydraulic conductivity of the EDZ with time The goal of this study is to strengthen the understanding of the phenomena involved in the EDZ formation and sealing in Opalinus Clay an indurated claystone currently being assessed as a host rock for a geological repository in Switzerland To achieve this goal hybrid finitediscrete element method FDEM simulations are performed With its explicit consideration of fracturing processes FDEM modeling is applied to the HGA experiment an in situ test carried out at the Mont Terri underground rock laboratory to investigate the hydromechanical response of a backfilled and sealed microtunnel A quantitative simulation of the EDZ formation process around the microtunnel is first carried out and the numerical results are compared with field observations Then the recompression of the EDZ under the effect of a purely mechanical loading capturing the increase of swelling pressure from the backfill onto the rock is considered The simulation results highlight distinctive rock failure kinematics due to the bedded structure of the rock mass Also fracture termination is simulated at the intersection with a preexisting discontinuity representing a fault plane oblique to the bedding orientation Simulation of the EDZ recompression indicates an overall reduction of the total fracture area as a function of the applied pressure with locations of ineffective sealing associated with selfpropping of fractures These results are consistent with hydraulic testing data revealing a negative correlation between pressure values and an increase in the EDZ transmissivityThis work has been supported by the Swiss National Cooperative for the Disposal of Radioactive Waste NAGRA and the Mont Terri Consortium representing the partners of the HGA experiment ANDRA France BGR Germany NWMO Canada NAGRA and Swisstopo SwitzerlandIn a 2D FDEM simulation the modeling domain is discretized with a mesh consisting of threenoded triangular finite elements Fournoded interface elements referred to hereinafter as crack elements are embedded between the edges of all adjacent triangle pairs An explicit finitedifference integration scheme is applied to solve the equations of motion for the discretized system and to update the nodal coordinatesThe elastic deformation of the solid rock is modeled according to the continuum theory of linear elasticity using constantstrain triangular elements Munjiza 2004 Lisjak et al 2014b If transversely isotropic conditions are assumed for instance in the case of bedded or schistose rocks the elastic formulation is characterized by five independent elastic constants two Young’s moduli E P and E S and Poisson’s ratios nu P and nu S for the directions parallel and perpendicular to the plane of isotropy respectively and the shear modulus G SConstitutive behavior of the crack elements a Normal bonding stress versus fracture opening for pure compressive/tensile loading Mode I b shear bonding stress versus fracture slip for pure compressive/tensile loading Mode II c elliptical coupling relationship between fracture opening ang slip for mixedMode fracturing Arrows denote loading unloading and reloading paths The area under the curves in a and b is defined by G mathrm Ic and G mathrm IIc respectivelyThe introduction of a finite crack element stiffness represented by the fracture penalty coefficient p f is required by the explicit nature of the FDEM solver For practical purposes the crack element contribution to the overall model stiffness can be largely limited by adopting penalty values that are very high yet practical in terms of the required time step size Munjiza 2004 Mahabadi 2012 Tatone 2014FDEM modeling of strength anisotropy after Lisjak et al 2014a a Linear variation of cohesive strength parameters with the angle gamma between crack element and layering orientation b Example of mesh combining a Delaunay triangulation for the intralayer material with edges preferentially aligned along the isotropy directionFrequency distribution of crack element orientation a random Delaunay triangulation and b Delaunay triangulation preconditioned along bedding In the preconditioned mesh a spike in the orientation distribution is associated with the bedding plane orientation psi =50 circ Also two symmetric peaks can be observed at approximately gamma pm 60 circ due to the adoption throughout the model of a uniform crack element size with triangles tending to an equilateral shapeSensitivity of the simulated fracture pattern to the adopted mesh topology Comparison of fracture pattern evolution for a a random Delaunay triangulation and b a Delaunay triangulation preconditioned along bedding A uniform element size h = 0015 m and same material models and input parameters were used in both cases The results presented in Sects 4 5 and 6 are based on mesh type b c Conceptual drawing illustrating the artificial meshinduced surface roughness introduced along a potential fracture path coincident with the bedding direction when adopting left a completely random mesh versus right a mesh preconditioned along the bedding

Keywords: