A specific upscaling theory of rock mass parameters exhibiting spatial variability: Analytical relations and computational scheme

ARTICLE IN PRESS

International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125 www.elsevier.com/locate/ijrmms

A specific upscaling theory of rock mass parameters exhibiting spatial variability: Analytical relations and computational scheme G. Exadaktylosa, M. Stavropouloub, a

Mining Engineering Design Laboratory, Department of Mineral Resources Engineering, Technical University of Crete, GR-73100 Chania, Greece b Department of Dynamic, Tectonic and Applied Geology, Faculty of Geology and Geoenvironment, University of Athens, Greece Received 15 May 2007; received in revised form 29 October 2007; accepted 30 November 2007 Available online 31 January 2008

Abstract In this paper the representation of geological conditions in a numerical simulation model is considered. By the expression ‘‘geological conditions’’ we mean the 3D volume geometry of the geological formations, the spatial variability exhibited by the rock parameters inside each of these geological volumes and the necessary upscaling of the rock deformability and strength parameters that are determined in the laboratory from cores collected in the field. A specific theory is outlined of how to go from laboratory tests, geological information and field measurements and observations to the full-scale numerical or ‘‘ground model’’ that includes, apart from initial and boundary conditions and ground water, the rock constitutive laws and associated material parameters for use in simulation models. The term ‘‘specific’’ used in the title of this paper stems from the fact that other possible approaches for the same problem may be applied; i.e. empirical rock mass classification systems, explicit modeling of joints in the rock by the distinct element or finite element methods, homogenization techniques, etc. The method used to take into account the spatial variability of rock mass properties by virtue of Geostatistical Theory and 3D modeling tools is also outlined, with an example case study. r 2007 Elsevier Ltd. All rights reserved. Keywords: Rock mass; Upscaling; Geostatistics; Numerical modelling; Damage Mechanics; Size effect

1. Introduction Numerical simulation models based on the finite element method (FEM), finite differences method (FDM), boundary element method (BEM), distinct element method (DEM), etc., are essential tools in the pre-feasibility and feasibility stages of a project, for aiding the selection of the best design and construction method of a tunnel, cavern, mining system, etc. A numerical model allows the engineer to test various design scenarios before construction commences and to assess with respect to economy and stability different alternative construction sequences or support designs. Thus, the underground excavation system may be constructed in ‘‘virtual space’’ (i.e. the computer) prior to excavating it in reality.

Corresponding author. Tel.: +30 210 7274778; fax: +30 210 7274096.

E-mail address: [email protected] (M. Stavropoulou). 1365-1609/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2007.11.008

For numerical models based on Continuum Mechanics such as the FEM or the FDM the input data required include the elastic stiffness tensor that in the case of isotropic conditions is expressed by the Young’s modulus and Poisson’s ratio of the rock mass, and the loading function based on a yield surface, with parameters describing it. The yield surface may be either isotropic (i.e. Hoek and Brown, Mohr–Coulomb, Drucker–Prager, etc.) or anisotropic (e.g. a multi-laminate model or others). The flow rule is also needed, describing the stress–strain behavior at yield and the hardening/softening function for non-perfectly plastic rock masses. Since rock masses are discontinuous there are two basic ways to represent jointed media in numerical simulation within the finite element or finite difference analyses, namely, joints are modeled as contacts using contact elements [1], or the jointed mass is replaced with equivalent solid medium and is simulated with a multi-laminate model [2]. In addition, knowledge of the virgin stress field (i.e. the stress tensor that exists in the

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

ground before excavation) and ground water conditions is required for the final assembly of the conceptual ‘‘ground model’’. On the top of the numerical simulation tools lies the geological model. This is because the cornerstone of a Rock Engineering project is the geological database upon which the definition of rock types, joints and material properties is based. The consideration of geological conditions in a numerical simulation program and appropriate rock mass parameters taking into account the presence of joints are important and challenging tasks, since geological and laboratory data have to be transformed into ‘‘digital information’’ that is required by the model. Even the most sophisticated simulation analysis can become a meaningless exercise if the geological information upon which it is based is inadequate or inaccurate. Here we describe a procedure for arriving at the input data required by the ground model. First, we analyze the information that is usually available at the outset of an underground engineering project. Then we look at the development of possible geomaterial constitutive models calibrated on standard tests in the lab that could be employed in a hierarchical fashion (i.e. from simple linear elastic and perfectly plastic laws of the Mohr–Coulomb type with five (5) parameters up to nonlinear yield laws with thirty (30) parameters and from isotropic to anisotropic cases with seven (7) parameters for each of the maximum three joint sets, and finally at the parameter determination at the scale of the underground excavation. This is performed with an appropriate upscaling theory that is based on the size effect exhibited by brittle or quasi-brittle rocks, fundamental Damage Mechanics principles [3–5] and rock mass classification systems. 2. Creation of the ‘‘digital’’ volume geological model 2.1. Handling of information obtained at design phase At the design phase of the project the following data are available: (a) The constitutive behavior of the rock at lab scale (i.e. of the order of 10 cm) in the frame of isotropic or anisotropic continuum elastoplasticity and elastoviscoplasticity theories, calibrated on several standard testing configurations (uniaxial and triaxial compression, indirect tension or direct tension, shear, etc.) on representative intact rock samples [6]. Some rock parameters may also be evaluated based on correlations with index tests that are more easily applied in the field on rock cores extracted from boreholes, such as Schmidt rebound hammer, point load or others. (b) The mechanical properties of the joints at the lab scale, e.g. normal and shear stiffness and strength properties, respectively, that are derived from the laboratory shear tests on representative specimens collected in the field or from cores.

1103

(c) From mapping on exposed surfaces, from (oriented) boreholes and from other exploration methods (i.e. geophysical surveys) the geometrical properties of the joints are collected and adequately described in a statistical context (probability density functions). These properties are the number of joint sets and their 3D orientation, persistence, frequency of joints and RQD, aperture, roughness, joint wall strength, water inflow, etc. The guidelines for the collection of these joint data are provided by ISRM [6]. These set of data may be subsequently used in rock mass classification systems, such as the rock mass rating (RMR), the tunnelling quality index (Q) or the geological strength index (GSI), that assign properly weight factors to each joint set property measured in the field in order to derive some rating or index of the rock mass quality. The flowchart of Fig. 1 depicts a possible design approach for underground excavations in rocks. It is worth noticing here that this flowchart complies with a certain path of the updated Rock Engineering Flowchart presented in Fig. 9 of the paper by Hudson and Feng [7]. One first starts with geological mapping and drill coring data or even seismic wave data from the field. Based on these data the conceptual geological model is created, and a procedure is employed to assign to every discretization element of the numerical simulation tool the lithology according to the geological model. This is necessary for assigning a constitutive law and parameters to each geological material, separately, in a subsequent stage. A geological longitudinal section for example (Fig. 2a) may be transformed into a ‘‘digital’’ volume model comprised from the individual volumes of each geological bed (Fig. 2b). Finally, a geological-discretized volume model is created as illustrated in Fig. 2d in which a given geological material is assigned to each grid cell of a finite element mesh (Fig. 2c). In a subsequent step this material is then assigned a material law and appropriate upscaled parameters. The next step is to assign to each element, apart from the lithology, the parameters needed in the constitutive model Field mapping, drill coring etc. Assembly of 3D geological model

Attribute lithology to the grid of the 3D numerical model

Evaluation of mechanical properties (upscaling) on the grid of the ground model

Observations during construction (feedback) Initial & Boundary conditions

Numerical computations with the ground model

Kriging code

Assign RMR (or Q or GSI) to the 3D grid

Fig. 1. Structure of the combined geological–geostatistical–numerical concept for tunnel design.

ARTICLE IN PRESS 1104

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

Fig. 2. Example of a ‘‘digital’’ geological model: (a) longitudinal geological section created from geological mapping and boreholes (courtesy of GISA, Spain); (b) geological block model in ACIS format; (c) import of geological volumes in GID pre-processor (courtesy of CIMNE, Spain) and discretised geological volumes and (d) discretized geological model with geological material assigned to grid points.

of each distinct lithological unit. The simplest way to do this is to use RMR (hereafter we consider the modified RMR proposed in [8]), Q [9] or GSI [10] rock mass quality indices to derive a single constant scalar damage factor for isotropic rock mass or a vector for anisotropic ones, which in turn is used with the values of the elasticity and strength parameters derived from small-scale lab experiments on intact rocks to perform the upscaling. The interpolation of the values of RMR, Q or GSI indices from the measurement points along the boreholes to each grid point (or Gauss integration point or element’s centroid) of the numerical model may be performed by the

Kriging technique and particularly with the aid of a special numerical code called KRIGSTAT, as was described in [11]. In the example under consideration the experimental semivariograms of the RMR data are found independently for combined GR1+Pf formations and GR2 formations, representing fresh granodiorite with interventions of porfyrite and weathered granodiorite, respectively. It is recalled here that the semivariogram of a regionalized variable, that is, functions whose values depend on spatial position, is commonly used to describe the spatial property of the domain (assuming of course that such a spatial dependence exists). It is a measure of the dissimilarity

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

between properties (the covariance between two values) that are located a distance h apart. The estimate of variance is repeated for many values of h to represent the spatial correlation of the property being analyzed. This similarity measure is denoted by the symbol g(h) and it is plotted on an xy plot with the x-axis being the distance h and g(h) on the y-axis. The experimental semivariogram is defined as half of the average squared difference between two attribute values approximately separated by vector hi, which is also called the lag. Subsequently, an exponential semivariogram model was selected from the library of models included in KRIGSTAT code and was fitted to each of the computed experimental semivariograms of each geological formation, that is    h 2 gðhÞ ¼ C 0 þ s 1  exp  (1) L in which C0 denotes the nugget, L is the range of influence and s represents the sill. The experimental semivariogram

1105

of the RMR of GR1+Pf formation and the best fit of the exponential model are shown in Fig. 3a. The validation (evaluation) of the Kriging model is performed in this case by virtue of the Q2 statistic and by the cross-validation technique [12]. The best validation results that are supported from geological considerations are obtained by using an anisotropic semivariogram model with an approximately four times larger correlation radius in the horizontal Oxy-plane compared to that in the vertical plane. As may be seen from Fig. 3b, the Kriging model passes this specific validation test successfully for this layer. Additionally, the results from a cross-validation analysis were also positive. In that case the geostatistical model is validated by re-estimating the known sampled values under implementation conditions as close as possible to those of the forthcoming production runs, i.e. semivariogram model, type of Kriging and the search strategy [13]. The semivariogram and Q1 test for the layer GR2 overlying GR1 formation are shown graphically in Figs. 4a and b. Also, in this case, Q1 and Q2 statistical tests referring to the

Singuerlin: GR1 layer, RMR 600 Experimental semivariogram Model semivariogram

Semivariogram, (h)

500 400 300 200 100 0 0

20

40

60 80 Lag, h [m]

100

120

140

Q2 validation: GR1 layer, RMR 4 3.5 3

Q2 pdf Q2 value Lower 95% confidence limit Upper 95% confidence limit

pdf(Q2)

2.5 2 1.5 1 0.5 0 0.5

0.6

0.7

0.8

0.9

1 Q2

1.1

1.2

1.3

1.4

1.5

Fig. 3. GR1 layer: (a) experimental semivariogram of RMR with a fitted exponential model (range ¼ 30 m, sill ¼ 360, nugget ¼ 25); (b) probability density function (pdf) of Q2 statistic and estimated Q2 value for n ¼ 151 and (c) cross-validation of the kriging model (n ¼ number of measurements).

ARTICLE IN PRESS 1106

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

Singuerlin: GR2 layer RMR 140 Experimental semivariogram Model semivariogram

Semivariogram, (h)

120 100 80 60 40 20 0 0

20

40

60

80

100

Lag, h [m] Singuerlin: GR2 RMR 6 Q1 pdf Lower 95% Upper 95% Q1 value

5

pdf(Q1)

4 3 2 1 0 -0.4

-0.3

-0.2

-0.1

0 Q1

0.1

0.2

0.3

0.4

Fig. 4. GR2 Layer: (a) semivariogram of RMR with a fitted exponential model (range ¼ 18 m, sill ¼ 75, nugget ¼ 20); (b) probability density function (pdf) of Q1 statistic and estimated Q1 value for n ¼ 151; (b) probability density function of Q2 statistic and estimated Q2 value for n ¼ 147.

semivariogram are successful. A similar anisotropic semivariogram model has been employed that was found to give better results compared with an isotropic model, as was done for GR1 formation. The validated semivariograms shown in Figs. 3a and 4a are subsequently exploited for the interpolation of the RMR rock quality index inside rock layers GR1, GR2 and of porphyry at the grid nodes of the imported GID mesh (e.g. Figs. 2c and d). The visualization of Kriging results referring to the RMR of granodiorite and porphyry in longitudinal sections and in the entire 3D model are displayed in Figs. 5a–c. The mean square estimation error (Kriging variance) spatial distribution is also displayed in Figs. 5a and b. In addition, the spatial distribution of RMR predicted by the model inside the tunnel and the estimated Kriging variance are illustrated in Figs. 6a and b. Subsequently, the initial and boundary conditions and field stresses are prescribed, ground water is added and the numerical or ‘‘ground model’’ based on finite elements, finite differences or coupled finite-boundary elements is

ready to perform the necessary computations for the various tunnel alignment, rock excavation and support scenarios (e.g. Fig. 1). After the commencement of tunnel excavation, the numerical model results are compared with available monitoring and face mapping data collected during the construction of the underground excavation. Monitoring data usually include convergence, strain and stress increment measurements, whereas face mapping data are used for the rock mass characterization. This set of ‘‘hard’’ data may be used for the upgrading of the ground model by applying appropriate back-analysis procedures (based on Simplex or other methods) (Fig. 1). If the upscaling problem is solved, then the definition of material parameters in the above manner in an FEM or an FDM model is fairly straightforward as each element (or indeed each Gauss Point inside an element) can be assigned different material properties. Since the BEM has a boundary discretization only, and no elements exist in the rock mass the inclusion of geological conditions requires further elaboration. As outlined in more detail in [14], the

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1107

Fig. 5. Spatial distribution of RMR and Kriging variance distribution in a GID mesh in (a) GR1+Pf formations, (b) GR2 formation and (c) spatial distribution of RMR in the entire model.

idea is to use internal cells for this purpose. These cells actually resemble finite elements but are significantly different because no additional degrees of freedom are introduced as the cells are only used to integrate the volume terms that arise. The internal cell mesh has the advantage that it can be unstructured (i.e. cells do not have to connect at nodal points as they have to do in the FEM) so that cell meshes are very easy to construct. The proposed procedure is to provide the ‘‘upscaling software’’ with coordinates of points and to receive the material law and parameters.

3. A specific upscaling theory 3.1. Hierarchical rock constitutive models Constitutive models for rock can be divided into elastic and elasto(visco-) plastic models. Elastic models can be isotropic or anisotropic. Isotropic models only need two parameters (E, n), whereas in the case of a stratified or a transversely isotropic material, five parameters (E, n, E0 , n0 , G0 ) are needed; in orthotropic rocks nine (9) parameters are needed, etc. [15].

ARTICLE IN PRESS 1108

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

Fig. 6. 3D perspective views of (a) RMR distribution of layers GR1 and GR2 and (b) mean square estimation error (Kriging variance) inside the tunnel.

However, elastic models are not very realistic for modeling the behavior of the rock mass since failure and inelastic deformations play a significant role. For example, the properties of joints can be significantly altered as asperities are sheared off after failure and during inelastic deformation. These phenomena can be modeled by appropriate plasticity evolution laws (for example strainhardening/softening behavior) or using the theory of Damage Mechanics or both. In most cases the behavior of the rock mass is affected by the presence of joints and is therefore essentially anisotropic. The term ‘‘joint’’ is used here to cover all discontinuities; these might technically be joints, faults, bedding planes, slickensides, compaction bands or other surfaces of weakness [16]. In many cases, especially when there are a large number of joint sets (i.e. more than three), the geomaterial behavior could be approximated with an isotropic elastoplastic model. The most well-known failure conditions that are used in rock mechanics are Mohr–Coulomb, Hoek and Brown and the multi-laminate models [17]. Of these three models only the last one considers the anisotropy of the rock mass. Gens et al. [18] proposed a simple model with a single yield surface, having linear elastic and perfectly plastic behavior with a linear yield criterion of the Mohr–Coulomb type, which requires only five parameters (two elasticity moduli, cohesion, angle of internal friction and dilation angle). The next most sophisticated model proposed by the same investigators is the model with up to three discontinuity families and with hardening/softening behavior and residual parameters and up to three families of discontinuities, requiring up to 30 parameters for the intact rock and seven additional parameters for each joint set plus a viscoplastic parameter for the joints. The hyperbolic M-C criterion (HM-C) tends asymptotically to the linear M-C model for high enough compressive hydrostatic stress and gives good predictions of the rock strength at low or even extensional hydrostatic stresses. This multi-axial criterion is defined by a hyperbolic yield surface with a rounded tip and is given by the

following expression [18]: f ¼ J 2  F 1 F p ¼ 0.

(2)

In this equation, F1 is related pffiffiffiffiffito the shape and position of the failure surface in the J 2  I 1 plane (I1 is the first invariant of the stress tensor sij, that is I1 ¼ skk; J2 is the second invariant of the deviatoric stress tensor sij, i.e. J 2 ¼ sij sji =2) while Fp relates to the same surface seen in the octahedral (or p) plane. The two functions are defined as follows: F 1 ¼ ðc  p tan fÞ2  ðc  pT tan fÞ2 ,

(3a)

F p ¼ ð1  Y cos 3yÞZ ;

(3b)

YZX0;

1oY o1.

In the expression for F1, c denotes the cohesion, f the internal friction angle, p ¼ I1/3 is the mean stress, I1 is the first invariant of the Cauchy stress tensor and pT is the socalled tension limit that is the strength of the rock in hydrostatic extension. In Fp, Y and Z are parameters that influence the shape of the trace of the yield or failure surface on the deviatoric plane and y denotes the Lode angle that is defined as [19] pffiffiffi 3 3 J3 cos 3y ¼ , (4) 2 J 3=2 2

wherein J3 is the third invariant of the deviatoric stress tensor and is given by 1 J 3 ¼ sij sjk ski ¼ s1 s2 s3 . (5) 3 In Fig. 7a, the HM-C failure criterion was fit to triaxial compression and indirect tension test data on Serena pffiffiffiffiffi sandstone, in which T ¼ J 2 denotes the deviatoric shearing stress intensity. In the same figure it may be observed that while the M-C failure criterion fits the uniaxial and triaxial compression test data, it fails to predict the tensile strength of the sandstone. The Fp on the deviatoric plane is shown in Fig. 7b. One may note that, for zero values of parameters Y or Z, the deviatoric crosssection becomes a circle and the hyperbolic Drucker– Prager surface would be recovered, while by setting pT ¼ c cot f the linear M-C yield criterion is obtained.

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1109

pffiffiffiffiffi Fig. 7. (a) M-C and HM-C fitted failure criteria of Serena sandstone in the J 2  I 1 plane (with negative compression); (b) representation of the HM-C failure criterion in the p-plane and si (i ¼ 1, y, 3) denote the principal deviatoric stresses (the trisectrice is normal to the paper and stress unit in MPa). The encircled numbers 1, y, 6 designate distinct regions with prescribed relative magnitudes of the principal stresses.

The failure criterion defined previously requires five parameters, namely c, f, pT, Y and Z. The expressions that give the parameters c and pT in terms of the uniaxial tensile strength (UTS), the uniaxial compressive strength (UCS), internal friction angle f, Y and Z may be derived from the following conditions: pffiffiffiffiffi I1 UCS (6) J 2 ¼ pffiffiffi ¼  pffiffiffi 3 3

we find

under uniaxial compression with y ¼ 601 and F p ðp=3Þ ¼ ð1 þ Y ÞZ , and pffiffiffiffiffi I1 UTS J 2 ¼ pffiffiffi ¼ pffiffiffi (7) 3 3

When the material is flowing inelastically the inelastic part of the deformation is defined by the flow rule that rests on the following constitutive assumptions: (i) principal axes of plastic strain rates are coaxial to the principal axes of the Cauchy stress and (ii) plastic strain rates are normal to the plastic potential surface q in stress space, that does not necessarily coincide with the yield surface. Hence, the

under uniaxial tension with y ¼ 601 and F p ð0Þ ¼ ð1  Y ÞZ . Solving the above two equations,

½3ð1  Y ÞZ  tan f2 UTS 2 þ ½tan f2  3ð1 þ Y ÞZ UCS2 , 6ðUCS þ UTS Þ tan f qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3c þ 9c2 þ 6cUCS tan f þ UCS 2 ðtan f2  3=ð1 þ Y ÞZ Þ . pT ¼ 3 tan f c¼ 

ð8Þ

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1110

plastic strain can be written in incremental form as follows: DðplÞ ij ¼ Dc

qq ; qsij

DcX0,

(9)

wherein the Greek capital letter D is used to denote increments, the superscript (pl) denotes plastic strains and c is the plastic multiplier. The rate form of the flow rule is essential to incremental plasticity theory, because it allows the history dependence of the response to be modeled. The plastic potential and the yield function may be identical, i.e. qf, only if the measured dilatation and strength responses of rock are identical. Such models are called ‘associated flow’ plasticity models. For most purposes a non-associated flow rule of the HM-C type is judged to be adequate, i.e. qq ¼ q0 dij þ q1 sij þ q2 Sij ; qsij   2 3=2 Sij ¼ J 2 sik skj  J 2s dij 3

ð10Þ

with 2 q0 ¼  tan cm ðc  p tan cm ÞF p , 3 q1 ¼ 1, pffiffiffi 3 3 F 1 F p ZY ð1  Y cos 3yÞ1 , q2 ¼  2

ð11Þ

wherein dij denotes Kronecker’s delta and cm denotes the mobilized dilatancy angle of the rock. The evolution of the yield surface is controlled by the evolution of each one of the parameters pT, c, f and cm in terms of the history variable c. Two types of history variables may be considered, i.e. accumulated deviatoric plastic strains and accumulated plastic work. The first history variable may be expressed adequately with the shear strain invariant, which is expressed in terms of the deviatoric plastic strain eðplÞ as follows: ij qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðplÞ ðplÞ ðplÞ gðplÞ ¼ 2eðplÞ (12) eðplÞ ij eji ; ji ¼ ji  kk dij =3.

Finally, from Eqs. (10), (11) and (13) the following relations are found: qv ¼  2 tan cm ðc  p tan cm ÞF p ,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

qg ¼ 2 F 1 F p f1 þ 3YZð1  Y cos 3yÞ1 ½1 þ 3YZð1  Y cos 3yÞ1 g. ð15Þ As has been proposed by Gens et al. [18] a hardening/ softening function is used which consists of three different parts, namely, initial hardening, softening and final residual. The hardening section is represented by a parabola, while the softening part by a polynomial of third degree. For the illustration of the theory we employ here a simplified version of the hardening/softening rule that is described simply by straight lines (e.g. Fig. 8) and is represented by the following parameters: F0 stands for cohesion, internal friction angle, tension limit or dilatancy angle at the initial yield state; Fp denotes peak cohesion or internal friction angle, tension limit or dilatancy angle, respectively; Fr stands for residual cohesion or internal friction angle, tension limit or dilatancy angle of the rock; gðplÞ is the upper bound of the shearing strain intensity p corresponding to the peak cohesion, or internal friction angle, tension limit or dilatancy angle; and gðplÞ represents r the lower bound above which the strength of the rock is characterized by the residual cohesion, internal friction angle, tension limit or dilatancy angle. 3.2. Determination of rock mass parameters After a model is selected it is necessary to determine the parameters for the model at the scale of the project (i.e. 1, 10, 100 m and so on). Considering the information available at design stage outlined in Section 2.1, the rock mass properties are determined by ‘‘upscaling’’ from laboratory specimens considering the effect of the geological conditions.

The plastic shearing strain intensity is an (isotropic) macroscopic measure of plastic slip that occurs at intergranular (or interblock) boundaries and across microcracks at the microscale or joints at the macroscale. According to the flow rule (9) we have DvðplÞ ¼ qv Dc;

DgðplÞ ¼ qg Dc;

qq , qskk sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi qq 1 qq qq 1 qq qg ¼ 2  dij  dij . qsij 3 qskk qsij 3 qskk qv ¼

ð13Þ From the above relations (13) the mobilized dilatancy parameter d can be found as follows: q d ¼ v. (14) qg

Fig. 8. Evolution law of strength parameter F representing cohesion, internal friction angle or tension limit, with respect to the history variable g(pl).

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

3.2.1. Size effect exhibited by the intact rock Before considering the joints that are added to the rock mass as we move from the lab scale (or mesoscale) to the scale of the discretization element (macroscale), one must take into account the ‘‘size effect’’ exhibited by the brittle or quasi-brittle rocks. It is recalled here that the class of semi-brittle solids is characterized by limited plastic flow, which is both due to nucleation and a sustaining mechanism for the growth of microcracks. Size effect is a subject of increasing interest due to the fact that current applications in rock mechanics involve a variety of length scales ranging from a few hundreds of microns (grain size) to centimeters (lab test specimens) up to several meters (scale of model discretization element). This range of scales and the corresponding necessity for modeling have revealed that there is a strong connection between the various length scales involved and that the response may depend on size for otherwise dimensionally similar structures. The size effect of intact rock strength is well known in the mechanics of brittle and quasi-brittle materials such as rocks and may be attributed to one of the following causes: (a) the failure of rocks in uniaxial compression in the form of spalling cracks with subsequent buckling of the individual rock slender columns that subsequently fail due to buckling [20]; (b) the transformation of stored volume deformation energy into fracture surface energy for the creation of new cracks [21]; (c) the fractal character of crack surfaces produced at rock failure [21]; or (d) the stress or strain gradient effect, e.g. [22–27]. In the simplest case of a scale- or size-dependent response, one may envision the interaction between the geometric length of the specimen/externally applied load and an internal length/internal force associated with the underlying microstructure (e.g. grain, microcrack, pore, etc.). The interaction between macroscopic and micro-

scopic length scales in the constitutive response and the corresponding interpretation of the associated size effects may be modeled through the introduction of higher order stress or strain gradients in the respective constitutive equations. In turn, it is quite interesting that one may find a statistical interpretation of the stress or strain gradients based on the stochastic character of stresses and strains due to heterogeneity of the representative elementary volume (REV) [28]. The experimentally observed size effect exhibited by the UCS of intact marble and sandstone specimens with the same height:diameter ratio equal to 2 is illustrated in Fig. 9. This size effect may be described by any of the models discussed above. For illustration purposes the mathematical model proposed by Bazant et al. [20]—although slightly corrected—may be employed, i.e. UCS ¼ UCS d þ C 1 H 2=5 ; pffiffiffi 3=5 4=5 ! 5ð2 3 Þ ðp2 E 0 K 4IC Þ1=5 , C1 ¼ 6

ð16Þ

wherein E 0 ¼ E=ð1  n2 Þ, KIC is the fracture toughness of the rock necessary for the formation of new mode-I cracks, H is the height of the specimen and UCSd is the ultimate UCS of the rock for sufficiently large height d of the rock (REV) with the geometrical similarity preserved (i.e. shape and height:diameter ratio). The above size effect law has been found to be valid for a marble and sandstone as it is displayed in Fig. 9. Also there is experimental evidence that this inverse 2/5 relationship has been found to hold true in borehole breakouts [29]. By inserting into Eq. (16) the representative elastic and fracture toughness properties of a given rock the constant C1 may be calculated and then the mean UCS d can be found from the respective representative value determined

Marble (-2/5) power law (marble) Sandstone (-2/5) power law (sandstone) Lac du Bonnet granite (-2/5) power law (granite)

250

200

UCS [MPa]

1111

150

100

50

0 0

10

20

30 H [cm]

40

50

60

Fig. 9. Size effect exhibited by the UCS of Dionysos marble and Serena sandstone [11], as well as best-fitted (2/5) size effect law on the experimental data for Lac du Bonnet granite [30].

ARTICLE IN PRESS 1112

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

at the lab scale (i.e. most frequently on specimens of 10 cm height and with a height to diameter ratio of 2:1), e.g. UCS d ¼ UCS 10  C 1 102=5 ,

(17)

whereas UCS10 in MPa is the value of the UCS determined from the laboratory tests on specimens of 10 cm height (with height:diameter ¼ 2:1) and C1 has units MPa cm2/5. For example, for a marble with E ¼ 67 GPa, n ¼ 0.3 and KIC ¼ 0.6 MPa m1/2, it turns out from the experimental results that C1 ¼ 84.8 MPa cm2/5, and for sandstone with E ¼ 30 GPa, n ¼ 0.1 and KIC ¼ 1.5 MPa m1/2, C1 ¼ 147.8 MPa cm2/5. From the above data it may be observed that the ratio of UCSd =UCS10 varies between 0.3 and 0.5. Hence, in the absence of lab tests at various scales, the following lower value of rock mass UCS may be employed to be representative for the intact rock UCS for design purposes: UCS d =UCS 10 ffi 0:3.

(18)

According to (a) the following parameters of Lac du Bonnet granite listed in Table 1 below (the first two taken [30] and the fracture toughness from [31]), and (b) the test data presented in [30] and shown in Fig. 9, the calibrated (2/5) size effect law gives C 1 ¼ 246:11 MPa cm2=5 . This value of C1 gives an asymptotic large-scale value of UCSd ¼ 100.7 MPa with UCSd/UCS10 ¼ 0.5, which coincides with the maximum value of the ratio as was mentioned above. The UCS350 i.e. at the size of the mineby tunnel (diameter ¼ 3.5 m) is predicted by the same formula to be approximately 125 MPa, which gives a ratio

Value

Young’s modulus, E (Gpa) Poisson’s ratio, v Fracture toughness, KIC (Mpa m1/2)

50 0.26 2.46

c1 UTS 10 ¼ UTS d þ pffiffiffiffi , L

(19)

where c1 is a constant with units of stress multiplied by the square root of length. For both marbles it turned out that UTS d =UTS 10 ffi 0:38.

Table 1 Avergage mechanical properties of Lac du Bonnet granite [30,31] Mechanical property

UCS350/UCS10 ¼ 0.63. Read et al. [32] on p. 263 mention that ‘‘macroscopic failure associated with notch formation initiated between 0.2 and 0.5 m behind the test mine-by tunnel where the peal tangential stress was only about 120 MPa, i.e. 0.6UCS10 from undamaged lab samples of Lac du Bonnet granite’’. One may note the agreement of the above estimation with that of the size-effect law proposed here. Cai et al. [33] suggested that the field rock strength at URL could be higher than 169 MPa and that the maximum wall stress for the mine-by tunnel is 169 MPa, under elastic and smooth tunnel wall condition. However, this strength should not be interpreted to be the UCS of the rock due to the ‘‘structural effect’’ caused both by the hole curvature and the stress-gradient effect of the hoop stress close to the boundary of a hole. It is expected that the UTS also exhibits a similar size effect. To test this hypothesis, two series of Four Point Bending tests on Lorano Carrara marble and Dionysos marble are considered, taken from [34]. The size effect is examined by utilizing four types of geometrically similar beams. The depth was kept constant while length L and height were changed maintaining the same ratio. As illustrated in Fig. 10, the following simple scaling law with the typical inverse square root scale dependence suggested by linear elastic fracture mechanics (LEFM) was employed:

(20)

Recapitulating, in this work we consider the values for both ratios UCSd =UCS10 and UTS d =UT S10 to be equal to 0.3, even if in the literature lower or higher (i.e. 0.6/0.7) values than this have been reported. If we further make the reasonable assumption that the dimensionless internal friction angle of the intact rock does

12

Flexural strength [MPa]

10 8 6 4

Test data, Lorano Carrara marble -(1/2) power law Lorano Carrara marble

2

Test data, Dionysos marble -(1/2) power law Dionysos marble

0 0

20

40

60 80 Beam length, L [cm]

100

120

Fig. 10. Size effect exhibited by the flexural strength of Lorano (Carrara) marble and Dionysos marble in four-point bending.

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

not exhibit a size effect, then according to Eq. (8) the largescale cohesion and tension limit of the HM-C failure model for the intact rock are given by the following relationships: cd ¼ 

pTd ¼

½3ð1  Y ÞZ=2  tan f2 UT S2d þ ½tan f2  3ð1 þ Y ÞZ=2 UCS2d , 6ðUCSd þ UT Sd Þ tan f

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   3cd þ 9c2d þ 6cd UCS d tan f þ UCS 2d tan f2  ð1þY3 ÞZ=2 3 tan f

considered through appropriate tensorial analysis). If the area dA with outward unit normal nj of the REV with position vector xi of Fig. 11a is loaded by a force dFi the usual apparent traction vector si ¼ sij nj is si ¼ lim

dA!S

.

ð21Þ 3.2.2. Effect of joints Next, we turn to the problem of the effect of joints in the rock mass. The effects of joints on rock mechanics properties include increased deformability, decreased strength, increased permeability, anisotropy, scale effects and additional kinematic possibilities. The upscaling problem common to many aspects of rock mechanics is translating knowledge of microcracks to knowledge of excavation project-scale fractures. Currently, it is reasonably feasible to study microfracture orientation, frequency and permeability characteristics in the laboratory, but this is not yet possible for the field scale, the scale of most interest. Attempts to scale between the laboratory and the field using power laws and fractal dimensions have sometimes proven successful, but not always. Fractal and power laws have proven successful when the geologic mechanisms causing the studied characteristics were similar on both scales. When the large-scale discontinuities were formed by different mechanisms, e.g., thermal stresses locally and tectonic globally, these methods have found limited success. Other possible approaches that one could employ to derive a model at the scale of the discetization element of the numerical model, called herein the ‘‘macroscale’’, are fracture mechanics (FM) and continuum damage mechanics (CDM). Although it may be shown that the two theories are equivalent, it is most appropriate to start with the latter for practical purposes since the former is based on not easily comprehensible parameters such as the stress intensity factor, each crack should be considered explicitly, nonlinearity is induced by the contact of joint lips and so on. Hence, a possible upscaling methodology may be founded on two basic hypotheses, namely: Hypothesis A. In a first approximation upscaling may be based on the scalar damage parameter D or vector damage parameter D ¼ Dn for the anisotropic case of joint-induced anisotropy of the rock mass (n is the unit normal vector of the plane of interest). It is noted that joints may also alter the constitutive law of the rock during the transition from the lab scale to the scale of the discretization element of the numerical model, i.e. they may induce anisotropy or they may induce nonlinearity and/or creep. Let us assume for the sake of simplicity that there are more than three joint systems in the rock mass so that it behaves like an isotropic material (the case of anisotropic geomaterial may be easily

1113

dF i ; dA

i ¼ 1; 2; 3,

(22)

where S is the representative area of the REV. The value of the dimensionless scalar damage parameter D(ni,xi) may be defined as follows: dAD , (23) dA where dAD is the total area occupied by the joints. In a similar fashion as it is displayed in Figs. 11b and c, the vector damage parameter may be referred to the directed area S n. At this point we may introduce an ‘‘effective traction vector’’ s~ i that is related to the surface that effectively resists the load, namely D¼

s~ i ¼ lim

dA!S dA

dF i ;  dAD

i ¼ 1; 2; 3.

(24)

From relations (22)–(24) it follows that si ; i ¼ 1; 2; 3. (25) s~ i ¼ 1D According to the above definitions the elastic deformation of the intact rock can be described by the relation s~ ij ! ðelÞ ij , which is obtained from elasticity, and the relation si ! s~ i , which is obtained by employing the concept of damage [3–5], wherein the superscript (el) denotes elastic strains. The damage parameter D should be multiplied with the joint closure factor Dc ð0pDc p1Þ in order to take into account partial contact of joint lips in compressive loading normal to the plane of the joint. In Appendix A, it is shown how this closure factor is related to the initial normal stiffness of the joint, the relative joint closure with respect to the maximum joint closure and the Young’s modulus of the intact rock. It is recalled here that Bandis et al. [35] have derived empirical formulae for maximum joint closure and initial stiffness in terms of index properties for the joint roughness coefficient (JRC), joint contact strength (JCS) and average aperture thickness. In Appendix A, it is shown that Dc is close to unity for low compressive stress regimes or for low initial normal stiffness of the joints; therefore from hereafter we consider Dc ffi 1. In thermodynamics, the axiom of the ‘‘local state’’ assumes that the thermomechanical state at a point is completely defined by the time values of a set of continuous state variables depending on the point considered. This postulate applied at the mesoscale imposes that the constitutive equations for the strain of a microvolume element are not modified by a neighboring microvolume element containing a microcrack or inclusion. Extrapolating to the macroscale, this means that the constitutive equations for the strain written for the surface dA  dAD

ARTICLE IN PRESS 1114

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

Fig. 11. (a) Representative elementary volume (REV) of damaged rock due to discontinuities, (b) specimen weakened by cracks with their surfaces normal to the Ox3 direction and (c) corresponding area reduction in Cauchy’s stress tetrahedron due to damage.

are not modified by the damage or that the true stress loading on the rock is the effective stress s~ i and no longer the stress si. The above considerations lead to the ‘‘Strain Equivalence Principle’’ [5], namely: ‘‘Any strain constitutive equation for a damaged geomaterial may be derived in the same way for an intact geomaterial except that the usual stress is replaced by the effective stress’’. The above principle may be applied directly on the isotropic elastoplastic geomaterial obeying the HM-C yield criterion. In the presence of damage, assuming that it is constant, the coupling between the damage and the plastic strain is written in accordance with the above principle, that is to say the yield function is written in the same way as for the non-damaged material, except that the stress is replaced by the effective stress according to Eq. (25), that is to say

1 f~ ¼ s~ij s~ji  ðcd  p~ tan fÞ2  ðcd  pTd tan fÞ2 F p 2 h 2 1 1 p tan f 3f~ ¼ sij sji  c  d 2 1D ð1  DÞ2

2 ðcd  pTd tan fÞ F p

3f~ ¼

1 fJ 2  ½ðcd ð1  DÞ  p tan fÞ2 ð1  DÞ2

 ðcd ð1  DÞ  pTd ð1  DÞ tan fÞ2 F p g.

ð26Þ

Hence, the relevant parameters of the failure function will be given by F~ 1 ¼ ðcm  p tan fÞ2  ðcm  pTm tan fÞ2

(27)

with cm ¼ cd ð1  DÞ;

pTm ¼ pTd ð1  DÞ.

(28)

Since for the hardening–softening model the same yield function is valid for the initial yield, ultimate failure and residual phase, the same relations (27) and (28) also hold true for the values of these rock mass strength parameters, respectively, i.e. cm0 ¼ cd0 ð1  DÞ; cmp ¼ cdp ð1  DÞ;

pTm0 ¼ pTd0 ð1  DÞ, pTmp ¼ pTdp ð1  DÞ,

cmr ¼ cdr ð1  DÞ;

pTmr ¼ pTdr ð1  DÞ.

ð29Þ

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

In addition, from relations (21) and (28) it is inferred that UCS m ¼ UCSd ð1  DÞ;

UTS m ¼ UTSd ð1  DÞ.

(30)

It may be noted that from the above arguments it turns out that the internal friction angle of the rock is not affected by the joints. This result is in agreement with the experimental data [e.g. 36], which indicated that the internal friction angle of granulated specimens of Wombeyan marble has the same value as the corresponding intact specimens but lower values of cohesion. From this fact the same investigators suggested that non-frictional processes play a greater role in the yielding of intact materials than in discontinua. Fig. 12 shows how the failure line in the p–T plane is affected by the damage of the rock mass. It may be observed that the above arguments from Damage Mechanics lead to a self-similar yield surface for the rock mass. For the linear elastic isotropic rock mass the elastic strain tensor is related to the effective stress tensor through the elastic stiffness of the uncracked rock as follows: ðelÞ ij ¼

1þn n s~ ij  s~ kk dij . E E

(31)

decrease, whereas closed cracks will cause it to increase [37,38], we will assume as a first approximation that the damage does not influence Poisson’s ratio. Mohammad et al. [39] performed a thorough review of a large number of case studies in which the lab parameters have been upscaled to the size of the numerical model. With regard to the Poisson’s ratio they inferred that the above approximation is valid, i.e. nm  n.

(35)

Next, the general equations for plasticity for the rock mass are derived. It turns out from Eq. (9) that the increment of the plastic strain tensor for the damaged rock is given as follows:  qq~ Dc 2 ðplÞ  tan cm Dij ¼ Dc ¼ qs~ ij 1  D 3 pffiffiffi 3 3 ðcm  p tan cm ÞF p dij þ sij  2 ½ðcm  p tan cm Þ2  ðcm  pTm tan cm Þ2    2 1=2 3=2 Z1 J 2 sik skj  J 2 dij ZY ð1  Y cos 3yÞ . 3 ð36Þ

Due to Eq. (25) the above relationship takes the form ðelÞ ij

1115

1þn n ¼ sij  skk dij Em Em

(32)

with the effective elasticity modulus given by the relation E m ¼ Eð1  DÞ.

(33)

Hence, the elastic stiffness tensor takes the form   Em 2n ðelÞ dij dkl þ dik djl þ dil djk . C~ ijkl ¼ 2ð1 þ nÞ ð1  2nÞ

Further, it may be easily shown that Prager’s consistency condition that is expressed as follows: qf qf Dc ¼ 0 Ds~ ij þ qs~ ij qc

results in the following relation for the plastic multiplier: Dc ¼ h1i

(34)

Although effective medium theories predict that damage in the form of open cracks will cause the Poisson ratio to

(37)

qf~ ðelÞ C Dkl =H p , qs~ ij ijkl

(38)

where Hp is the plastic modulus defined as H p ¼ H t þ H 0 40

(39)

14 12 10 8 Damage increases

D=0 D=0.5 D=0.8

Deviatoric shearing stress intensity, T

16

6 4 2

tension limit pT, intact rock

0 Mean stress, p Fig. 12. Schematic representation of the hyperbolic Mohr–Coulomb criterion applied to damaged geomaterial in the T–p plane (y ¼ p/3) (compressive stresses are taken as negative quantities here).

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1116

with

(e.g. Eq. (26)):

qf~ ðelÞ qq~ C , qs~ kl klmn qs~ mn q~ g qf~ qf~ ¼ Ht ¼  qc ð1  DÞ qgðplÞ

f~ ¼

H0 ¼

and h1i ¼

8 > < 1 if f~ ¼ 0 ^

qf~ qs~ ij

ð40Þ

C ðelÞ ijkl Dkl 40; qf~ qs~ ij

> : 0 if f~o0 or f~ ¼ 0 ^

(41)

C ðelÞ ijkl Dkl p0;

wherein

1 fJ 2  ½ðcm ðgðplÞ Þ  paðgðplÞ ÞÞ2 ð1  DÞ2  ðcm ðgðplÞ Þ  pTm ðgðplÞ ÞaðgðplÞ ÞÞ2 F p g

in which cm ðgðplÞ Þ; pTm ðgðplÞ Þ; aðgðplÞ Þ indicate the hardening (softening) functions of the rock mass strength parameters with respect to the plastic shearing strain intensity g(pl) and a ¼ tan f denotes the stress inclination in the p–T plane. The above plasticity equations are necessary because based on them it may be shown that the effect of damage is

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    3YZ 3YZ 2 2 1þ . q~ g ¼ 2 ½ðcm  p tan cm Þ  ðcm  pTm tan cm Þ F p 1 þ ð1  Y cos 3yÞ ð1  Y cos 3yÞ The second relation of (40) was derived by recourse to the second of relationships (13) and (36). In case of hardening H t 40, while for H t o0 softening is said to take place. The gradient to the failure surface of the jointed rock mass is given by the equation

ð44Þ

(42)

to reduce the history parameter g(pl) as follows: g~ ðplÞ ¼ ð1  DÞgðplÞ .

(45)

ð43Þ

This means that the hardening/softening function of the intact geomaterial contracts in a self-similar fashion as damage increases, since both the strength and the history parameters deteriorate by the same amount. This selfsimilarity property directly derived from the theory is illustrated in Fig. 13. Therefore, the most important result of the above analysis is the demonstration that the constitutive equation of the rock mass may be derived from the respective model inferred from laboratory experiments, by considering the size effect and the effect of damage as it was described above.

Furthermore, the derivative qf~=qgðplÞ appearing in the second of Eq. (40) can be easily found by differentiating the following expression of the yield function

Hypothesis B. A universal relationship between the damage parameter D and RMR (or equivalently with Q or GSI rock quality indices) exists for all rock masses.

 qf~ 1 2  tan fðcm  p tan fÞF p dij þ sij ¼ qs~ ij 1D 3 pffiffiffi 3 3 ½ðcm  p tan fÞ2  ðcm  pTm tan fÞ2   2   2 1=2 3=2 Z1 ZY ð1  Y cos 3yÞ J 2 sik skj  J 2 dij . 3

2.5

2

1.5 F/F0

D=0 (intact) D=0.5 D=0.9

1

0.5

0 0

0.5

1

1.5

2 2.5 g(pl)/g(pl)-peak

3

3.5

4

4.5

Fig. 13. Self-similar contraction of the hardening/softening function of the strength parameter (i.e. cm, pTm, f, cm) with increasing damage levels of the rock.

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

This hypothesis is based on the fact that RMR, Q or GSI indices do take into account explicitly the most important features of the rock mass joints that are responsible for the deterioration of its parameters compared to those of the intact state; hence, they may be linked to the damage parameter D with an appropriate relationship. Such a function must have a sigmoidal shape resembling a cumulative probability density function giving D in the range of 0–1 for RMR or GSI varying between 100 and 0 or for Q varying from 1000 to 0.001, respectively. On the grounds of the above considerations, the Lorentzian cumulative density function (obtained directly from the Cauchy probability density function) is proposed, namely (    ) ^ b^ RMR  c p D ¼ 1  a^ þ tan1 þ , (46) p 2 d^ where c^ denotes the location parameter, specifying the location of the peak of the distribution, d^ is the scale parameter which specifies the half-width at half-maximum, respectively, b^ is a proportionality constant and a^ can be found in terms of the other three parameters by setting D ¼ 1 for RMR ¼ 100. It is noticed that (46) is not used here as a statistical function but merely as a deterministic function. In this formula the RMR does not include the correction term due to unfavorable tunnel orientation with respect to joints and the ground water (hence RMR considers only the joints, i.e. RMR89). As it is illustrated in Fig. 14 the ^ c^; d^ were three unknown parameters of this function b; calibrated on the empirical relationship proposed by Hoek and Brown [40] for the dependence of UCS m =UCS d ¼ 1  D on RMR, with UCS m being the rock mass UCS, namely ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi UCS m pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (47) ¼ eððRMR5Þ100Þ=9 ¼ eðGSI100Þ=9 . UCS 10

Based on this best-fit procedure, performed for RMRo80, the following values for the unknown parameters were found, namely c^ ¼ 72:9968; d^ ¼ 12:9828; and a^ ¼ 0:06111.

b^ ¼ 1:2377 ð48Þ

These values of parameters mean that D ffi 0:5 for RMR ffi 70.

(49)

In the same figure the equation proposed in [41] has also been plotted for comparison purposes. It may be observed from this figure that in contrast to the other two empirical relations, for RMR ¼ 100 i.e. for intact rock, the ratio UCSm =UCS10 is not equal to 1 but to 0.3, which is the assumed ratio of the large-scale UCSd over the lab-scale UCS10 (i.e. UCSd =UCS10 ). It is noted here that the calibration of the Lorentzian curve on Hoek’s and Brown’s test data does depend on the value of the ratio UTS d =UTS 10 ; however, the prediction of strength of the rock mass does not depend on the ratio at hand for RMR values below 80. In case some investigator assumes a ratio UTS d =UTS 10 ¼ 0:6 for example, then another calibration of the Lorentzian function parameters that gives UTS m =UTS 10 ¼ 0:6 for RMR tending to 100 should be performed. If Hypothesis A holds true, then for a rock mass exhibiting isotropy and linear elasticity, the damage parameter D influences to the same extent both the strength and the Young’s modulus as it is explicitly demonstrated by Eqs. (29), (30) and (33), respectively. In Fig. 15, the ‘‘integrity’’ of the rock mass (equal to (1D)) with respect to RMR as is predicted by the proposed empirical Lorentzian law (Hypothesis B) is plotted. It is recalled here that in contrast to strength, elastic moduli do not display a size effect that is Em/E tends to 1 for RMR tending to 100. The most important result that gives additional validity to the Hypotheses A and B comes from

0.5 Lorentzian law Ramamurthy & Arora (1994)

0.4

UCSm/UCS10

Hoek & Brown (1997) 0.3

0.2

0.1

0 0

20

1117

40

60

80

100

RMR Fig. 14. Calibration of the Lorentzian damage law on Hoek’s and Brown’s empirical equation.

ARTICLE IN PRESS 1118

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125 1 Lorentzian law Hoek & Diederichs (2006)

0.9 0.8

Integrity, 1-D

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

10

20

30

40

50 RMR

60

70

80

90

100

Fig. 15. Plot of rock mass integrity (1D) with respect to RMR based on Lorentzian law and on empirical equation pertaining to the deterioration of Young’s modulus proposed in [42].

recent results published in [42]. In this paper, the authors exploit a large database of in situ measurements of deformation moduli of rock masses and GSI from China and Taiwan have proposed a new relationship between Em and GSI, namely E m ¼ Ef0:02 þ ½1 þ eð60GSIÞ=11 1 g.

(50)

The above empirical relationship of the moduli ratio with GSI substituted by RMR89 is plotted in Fig. 15 for comparison with the relationship proposed here. The close agreement of the two relationships is illustrated in this figure. Therefore, it may be remarked here that the empirical relationships proposed by the above-referenced investigators regarding the deterioration of the strength and deformability of rock mass treated distinctively are equivalent and hold true both for the strength and deformability in a unique fashion as it is predicted by the Damage Mechanics theory. Furthermore, D as given by relationships (46) and (48) can be expressed as a function of the Q index by virtue of the following empirical relationship linking the latter index with RMR: Q  10ðRMR50Þ=15 .

(51)

Another attractive feature of the Damage Mechanics approach followed here is that the second hypothesis stated above may be discarded if the damage parameter D of rock mass could be directly estimated from seismic P or S-wave velocities with seismic measurements in boreholes, employing either the cross-hole technique or the suspension method developed in [43]. In the first method the in situ velocity of seismic waves is L V~ S ¼ , t

(52)

where L is the distance traveled by the S-wave between the receiver and the transmitter in the neighboring boreholes and t is the first arrival time of the S-wave at the receiver. In the second method, the soil or rock mass volume involved in the measurement is relatively small (L ffi 1 m), a rather sophisticated wave generation process is applied and a more direct measurement of local shear wave velocity is performed. Assuming again that Poisson’s ratio remains unaffected by the damage, the damage parameter can be estimated with the following simple Damage Mechanics formula [5]: D¼

2 2 V~ V~ Gm E m ¼ ¼ 1  S2 ¼ 1  P2 , G E VS VP

(53)

where Gm denotes the shear modulus of the rock mass, G the shear modulus of the intact rock and VS and VP denote the shear and compressional wave velocities, respectively, measured on intact rock cores in the laboratory. In order to evaluate the applicability of formulae (46) and (48) that are the bases of Hypothesis B, we exploit the following empirical formula proposed in [9]: V~ P  3:5 þ log10 Q.

(54)

Taking VP ¼ 6500 m/s, as has been suggested in [44], then the damage of the rock mass can be estimated through formula (53) in terms of Q; further, the damage with respect to Q can be independently estimated from the Lorentzian law (46) by virtue of the inverse relation (51) linking RMR with Q. The comparison of the damage predicted by the above two methods is illustrated in Fig. 16. It may be seen from this figure that both are in good agreement. This observation is encouraging for the future estimation of D through seismic wave measurements.

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1119

1 0.9

Damage Mechanics formula

0.8

Barton's (2002) empirical eqn

Damage parameter, D

0.7 0.6 0.5 0.4 0.3 0.2 0.1

0.0001

0.0010

0.0100

0 0.1000 1.0000 10.0000 Tunnelling Quality Index, Q

100.0000

1000.0000

10000.0000

Fig. 16. Dependence of damage (D) on Q rock quality index as is independently predicted by Barton’s empirical formula [9] and Damage Mechanics. Table 2 Characterization of the rock masses at the Kannagawa and Kazunogawa test sites using the GSI system and the damage theory Geological formations CG1 UCS (Mpa) mi c (Mpa) f (deg) E (Gpa)a GSI D cm (Mpa) fm (deg) Em (GPa)

111 22 19.6 50.8 50 74.0 0.25 4.4 50.8 37.4

Test data

CG2

5.2b 57.0b 45.3 (6.2)c

162 19 29.1 49.5 72.9 65.0 0.51 4.25 49.5 35.5

Test data

FSI

3.4 57.0 33 (3.6)

126 19 22.6 49.5 56.7 65.0 0.51 3.3 49.5 27.6

Test data

MI

3.4 57.0 24.4 (2.5)

48 9 9.9 41.8 21.6 54.0 0.73 0.81 41.8 5.9

Test data

CH

1.9 40.0 11.7 (1.7)

108 19 19.4 49.5 48.6 60.0 0.64 2.12 49.5 17.7

Test data

CM

Test data

1.5 58.0 12.9 (2.84)

108 19 19.4 49.5 48.6 46.0 0.84 0.91 49.5 7.6

0.8 55.0 7.9 (1.2)

a

Estimated from formula (60). Estimated from in situ block shear tests. c Estimated from in situ plate bearing tests. Number in parenthesis indicated the standard deviation. b

4. Evaluation of the specific damage theory predictions on in situ experiments In order to compare the above upscaling theory comprised of the two fundamental hypotheses, we first consider the data presented in [45]. In that paper the GSI system was applied to characterize the jointed rock masses at two underground powerhouse cavern sites in Japan. GSI values were obtained from the block volume and joint condition factors, denoted by the symbols Vb and Jc, respectively, which in turn are determined from site construction documents and field mapping data. Based on GSI values and intact rock strength properties, equivalent Mohr–Coulomb strength parameters and elastic modulus of the jointed rock mass were calculated and

compared to in situ test results from two case studies. The first one is the Kannagawa pumped hydropower project located in Gumma Prefecture in Japan. The rock mass at the site consists of conglomerate, sandstone and mudstone. The percentages of conglomerate in CG1 and CG2 rock mass zones are about 93% and 62%, respectively, and FSI and MI stand for sandstone and mudstone formations. The second one is the Kazunogawa power station in which the rock mass consists of sandstone and composite rock of sandstone and mudstone, described as two groups (CH and CM) of rock mass types. Table 2 presents the data for these two projects. UCS has been estimated from laboratory tests, mi is the parameter used in the Hoek–Brown failure criterion [40] and has been estimated from laboratory triaxial compression tests. The equivalent Mohr–Coulomb

ARTICLE IN PRESS 1120

G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

along with the respective measured UCSm are shown in Table 3 for comparison purposes.

cohesion and internal friction angle of each rock type have been estimated and presented in Table 2 below. The authors did not present any data for the Young’s modulus of the rocks derived from laboratory tests. For this purpose we employ the following empirical relationship between Young’s modulus and UCS10 deduced from regression analysis of experimental data from a variety of rocks [46]: E ffi 400C500  UCS10 .

5. Final numerical model ready for simulations Returning back to the creation of the ‘‘digital ground model’’, based on the above specific upscaling theory the ‘‘mapping’’ of material parameters in the above manner at each element of the numerical model shown in Fig. 2d is fairly straightforward. In order to illustrate the effectiveness of the proposed methodology, the GID discretization mesh with lithology and RMR (e.g. Figs. 2c, d and 5a–c) was first successfully transferred to FLAC3D [48] as it is displayed in Figs. 18a and b, respectively. In a subsequent step the parameters of the GR1, GR2 and Pf formations determined in the lab by fitting to the experimental results obtained from a series of Brazilian, uniaxial compression and triaxial compression tests on intact rocks a linear elastic and perfectly plastic linear Mohr–Coulomb model with tension cut-off (i.e. only two elasticity moduli and three strength parameters assuming zero dilatancy angle) were reduced based on the spatial distribution of RMR and relationship (46) linking D with RMR. The representative mechanical properties of the intact rocks at lab scale are shown in Table 4. As an example of the results of the above pre-processing analysis, the spatial distribution of rock mass cohesion and bulk modulus in the volume of the model are displayed in Figs. 18c and d, respectively. Since there are no distinct faults or major discontinuities in these formations that could be modeled explicitly, the ‘‘ground model’’ is now ready for simulations of the various excavation scenarios for the investigation of the rock mass behavior provided that ground water is added plus initial stress state and boundary conditions.

(55)

Also, from the empirical relationship RMR  GSI þ 5

(56)

the constant damage parameter D can be estimated for each geological formation from the GSI value and formula (46). The estimations of cm, fm, Em based on the above Damage Mechanics theory are illustrated in Table 2 in order to be compared with in situ measurements performed on large-scale specimens from these geological formations. It is recalled here that according to the previously presented damage theory the internal friction angle of the rock mass remains unaffected by the joints. In a subsequent step, the eight sets of in situ UCS data presented in the dissertation of Palmstrom [47] are considered. These data included large samples tested in a laboratory, in situ tests on pillars in a mine (one ‘sample’) and tests combined with back analysis from a slide in a mine. Also, for these cases the values found for the block volume (Vb) and the joint condition factor (Jc) for the samples have been found and presented in Table 3. By virtue of the GSI quantitative diagram based on the two indices Vb and Jc, the GSI regions of each case have been delimited and presented in Fig. 17. For the estimation of the Damage parameter D for each test case we employ the minimum GSI value as indicated in Table 3. Then, the UCSm can be predicted by taking into account the UCS of the intact rock (presented in Table 3 for each case), the size effect and the damage imparted to the large-scale specimens by the joints based on Eq. (46). The predicted UCSm

6. Conclusions In this paper, we have presented a specific upscaling theory on how information that is available about the rock

Table 3 Characterization of the rock masses using the GSI system and the damage theory Case no

Case studies

UCS (MPa)

Vb (cm3)

Jc

GSI

Integrity (1D)

Predicted UCSm (MPa)

1 2 3

Stripa granite Panguna andesite Sandstone pillar, Laisvall mine 100,000–300,000

200 265 210 0.75–2.25

5000–15,000 2–6

1.5–2.5 4–6

45 35

0.17 0.08

10.46 6.93

7.55 3.7

160

50 8000–20,000

0.22 0.2–0.3

30

0.06

2.80

1.85

55

5000–10,000

1.5–2

45

0.17

2.88

3

100

5000–10,000

2–2.5

47

0.19

5.91

8

8 65

1000–5000 5000–10,000

5–10 3.5–4.5

55 55

0.29 0.29

0.73 5.92

1.33 6.8

4 5 6 7 8

Langsele mine, grey schist and greenstone Caledonian clay-schist parallel to schistocity Caledonian clay-schist perpendicular to schistocity Mesozoic sandstone Paleozoic siltstone

14.74

Measured UCSm (MPa)

20

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1121

Fig. 17. Quantitative GSI chart [43] and the seven considered case studies.

mass from lab tests and in situ surveys can be converted into data that are suitable for a numerical model. Linking the geological conceptual model assembled from geological mapping, core drilling investigation and other measurements (e.g. geophysical data, borehole logging measurements, etc.), with the numerical model in the 3D space is essential for more realistic modeling of excavations in rock masses. In addition, a methodology for the representation

of geological information in the numerical model has been given. The above approximate theory for the estimation of model parameters may be very useful in rock engineering, since: (i) The state of joints distributed in the rock mass can be represented by the constant damage parameter which

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1122

Block Group Live mech zones shown gr1 qcs pf gr2

Live mech zones shown 2.5000e+001 to 3.0000e+001 3.0000e+001 to 3.5000e+001 3.5000e+001 to 4.0000e+001 4.0000e+001 to 4.5000e+001 4.5000e+001 to 5.0000e+001 5.0000e+001 to 5.5000e+001 5.5000e+001 to 6.0000e+001 6.0000e+001 to 6.5000e+001 6.5000e+001 to 7.0000e+001 7.0000e+001 to 7.5000e+001 7.5000e+001 to 7.8000e+001 Interval = 5.0e+000

Block contour of RMR

1.8500e+004 to 5.0000e+005 5.0000e+005 to 1.0000e+006 1.0000e+006 to 1.5000e+006 1.5000e+006 to 2.0000e+006 2.0000e+006 to 2.5000e+006 2.5000e+006 to 3.0000e+006 3.0000e+006 to 3.5000e+006 3.5000e+006 to 4.0000e+006 4.0000e+006 to 4.5000e+006 4.5000e+006 to 4.7027e+006 Interval = 5.0e+005

7.0915e+006 to 2.5000e+009 2.5000e+009 to 5.0000e+009 5.0000e+009 to 7.5000e+009 7.5000e+009 to 1.0000e+010 1.0000e+010 to 1.2500e+010 1.2500e+010 to 1.5000e+010 1.5000e+010 to 1.7500e+010 1.7500e+010 to 2.0000e+010 2.0000e+010 to 2.2500e+010 2.2500e+010 to 2.3472e+010 Interval = 2.5e+009

Fig. 18. (a) Imported rock formations in FLAC3D; (b) spatial distribution of RMR inside GR1+Pf and GR2 formations on the vertices of tetrathedra of an FLAC3D discretized model; (c) distribution of cohesion (in Pa) in the GR1, Pf and GR2 formations and (d) distribution of bulk modulus (in Pa) in the GR, Pf and GR2 granite formations.

may be scalar or vector, and its mechanical effects are described by introducing the net stress which acts on the ‘‘net’’ surface excluding the joint surfaces. (ii) The constitutive equation is required only for the intact rock, which is obtained from a calibration procedure of the model on laboratory test results. Then the constitutive equation of the rock mass is derived by considering the size effect and the effect of damage as it was described above.

(iii) It is obvious that the ‘‘engineering judgment’’ plays a role in the process of assessing the damage of the rock mass. (iv) Rock mass characterization systems are placed in a ‘Mechanics of Materials’ context since they are linked with the damage parameter, which is a central concept in Damage Mechanics theory. (v) The damage parameter to be incorporated in to the constitutive rock models in the numerical simulation

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

(vi)

(vii)

(viii)

(ix)

(x)

(xi)

codes is linked with empirical approaches, which are faster although much less accurate, based on rock mass classification systems for rock support measures. RMR, GSI or Q rock mass quality indices that quantify the rock mass structure and sampled from the boreholes or from rock mass exposed surfaces may be interpolated on a regular or an irregular mesh grid of each geological formation separately by virtue of measurements from boreholes or rock exposures and the Kriging technique [12,13]. Hence, damage parameter spatial distribution is also evaluated based on relationship (46) that links damage with RMR. The disadvantage of the rock mass classification systems approach, namely that complex properties of a rock mass cannot be satisfactorily described by a single number, is surpassed with this approach in which the complexity of the constitutive model is considered. Also the hindrance that the same rating can be achieved by various combinations of classification parameters, even though the rock mass behavior could be different, is also eliminated because we may choose a different constitutive law based on experimental evidence. The other drawback that the user is led directly from the geological characterization of the rock mass to a recommended ground support without the consideration of possible failure modes is also eliminated. Using the proposed upscaling scheme the initial values of the model parameters for a back-analysis algorithm aimed at identifying these parameters are given, hence facilitating the quick convergence of the algorithm. Finally, as it was demonstrated here, the in situ seismic wave measurements may be used for the evaluation of the damage state of geological formations in order to bypass the second hypothesis of the relation between D and RMR.

Benchmark tunnel sites have already been considered for the application of this type of analysis. The work performed in the project will lead to realistic results and knowledge of the applicability of the proposed methodology for certain problems. These benchmark cases will also indicate the applicability of the proposed upscaling

1123

relations for the deformability and strength of rock masses based on already existing rock mass indices (RMR, Q, GSI or other) and the Damage Mechanics theory. Acknowledgments The authors acknowldge the financial support from the EC 6th Framework Project TUNCONSTRUCT (Technology Innovation in Underground Construction) with Contract no. NMP2-CT-2005-011817. They also express their appreciation to Henning Schwarz of GISA, who kindly provided data from the L9 case study for the demonstration of the methodology. The authors also acknowledge financial support by PYTHAGORAS through the contract MIS 97510 ‘‘PYTHAGORAS II: Funding of research groups in Technical University of Crete-M2.2’’, program ‘‘Operational Program for Education and Initial Vocational Training’’ (O.P. ‘‘Education’’), Third Community Support Framework co-funded by the European Social Fund and National Resources— (EPEAEK-II) PYTHAGORAS. Appendix A In Section 3.2.2, the effective stress was defined as follows: si s~ i ¼ ; i ¼ 1; 2; 3. (A.1) 1D

Fig. A1. Specimen with a mated joint subjected to uniaxial compression.

Table 4 Representative mechanical properties of the intact rocks at the laboratory scale Test results

Uniaxial compressive strength Tensile strength Young’s modulus Poisson’s ratio Internal friction angle

Description

Average Average Average Average Average

(MPa) (MPa) (GPa) (deg)

Rock type Fresh granodiorite (GR1)

Porphyry (Pf)

Weathered granodiorite (GR2)

102 11 49 0.25 42

70 11 45 0.25 40

15 1.8 7.5 0.3 35

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

1124

Fig. A2. Dependence of the joint closure parameter on relative joint deformation.

This expression holds true both for tension and for compression if the joints remain open. But above some compressive load the joints may close and so the effective area in compression S~ is such that ~ S  S D oSoS,

(A.2)

where S~ ¼ S  Dc S D ,

(A.3)

wherein Dc is the joint closure parameter. From relation (A.3) it may be seen that the actual damage parameter in cases of joint closure will not be D but equal to DcD. Still damage D can be found by the empirical relationship (48). The stress–displacement constitutive relation s ¼ s ðdÞ, wherein d denotes displacement, of the jointed sample of Fig. A1 may be written as follows: s ¼

E ð1  Dc Þd. ‘

(A.4)

Bandis et al. [35] showed that for mated (correlated lips) joints (Fig. A1) the following hyperbolic empirical law fitted their experimental data well: s ¼

d , a  bd

(A.5)

where a=b ¼ d m with dm denoting the maximum joint closure. By setting 1 ¼ K n0 a

(A.6)

we may prove that Kn0 is the initial normal joint stiffness for zero compressive load. It may also be shown that the normal stiffness of the joint is governed by the following

relationship: Kn ¼

qs K n0 ¼ . qd 1  ðd=d m Þ2

(A.7)

Then from Eqs. (A.3)–(A.7) it can be proved that the following expression of the damage parameter Dc in terms of the joint properties, characteristic rock volume size and Young’s modulus of the intact rock holds true: Dc ¼ 1 

lK n0 . E½1  ðd=d m Þ

(A.8)

This expression for the joint closure parameter is plotted in Fig. A2 for three values of the ratio lKn0/E. References [1] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. Proc ASCE 1968;94(EM3):637–59. [2] Zienkiewicz OC, Pande GN. Time dependent multi-laminate model of rocks—a numerical study of deformation and failure of rock masses. Int J Numer Anal Methods Geomech 1977;1:219–47. [3] Kachanov LM. Introduction to continuum damage mechanics. Dordrecht: Kluwer; 1986. [4] Krajcinovic D. Damage mechanics. Mech Mater 1989;8:117–97. [5] Lemaitre J. A course on damage mechanics. Berlin: Springer; 1992. [6] Brown ET, editor. Rock characterization, testing and monitoring: ISRM suggested methods. Oxford: Pergamon; 1981. [7] Hudson JA, Feng XT. Updated flowcharts for rock mechanics modeling and rock engineering design. Int J Rock Mech Min Sci 2007;44:174–95. [8] Bieniawski ZT. Engineering rock mass classifications. New York: Wiley; 1989. [9] Barton N. Some new Q-value correlations to assist in site characterization and tunnel design. Int J Rock Mech Min Sci 2002;39: 185–216. [10] Hoek E. Strength of rock and rock masses. ISRM News J 1994; 2(2):4–16. [11] Stavropoulou M, Exadaktylos G, Saratsis G. A combined threedimensional geological–geostatistical–numerical model of under-

ARTICLE IN PRESS G. Exadaktylos, M. Stavropoulou / International Journal of Rock Mechanics & Mining Sciences 45 (2008) 1102–1125

[12] [13] [14]

[15] [16] [17] [18]

[19] [20]

[21] [22] [23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

ground excavations in rock. Rock Mech Rock Eng 2007;40(3): 213–43. Kitanidis PK. Introduction to geostatistics. Cambridge, MA: Cambridge University Press; 1997. Deutsch CV, Journel AG. GSLIB: geostatistical software library and user’s guide. New York: Oxford University Press; 1992. Beer G, Duenser C. Advanced numerical simulation of the tunnel excavation/construction process with the boundary element method. In: Proceedings of the ISRM congress, Lisbon, 2007. Lekhnitskii SG. Anisotropic plates. New York: Gordon & Breach; 1968. Pollard DD, Aydin A. Progress in understanding jointing over the past century. Geol Soc Am Bull 1988;100:1181. Beer G. Programming the boundary element method. New York: Wiley; 2001. Gens A, Carol I, Caballero A, Gonza´lez N, Segura JM. Specifications for web-based library of constitutive models for rock and soft soils. Report on model classes. Technology Innovation in Underground Construction (TUNCONSTRUCT) Project Deliverable D1.1.1.3, 2005. Chen WF, Han DJ. Plasticity for structural engineers. New York: Springer; 1988. Bazant ZP, Lin F-B, Lippmann H. Fracture energy release and size effect in borehole breakout. Int J Numer Anal Methods Geomech 1993;17:1–14. Bazant ZP, Chen E-P. Scaling of structural failure. Appl Mech Rev 1997;50(10):593–627. Aifantis EC. Strain gradient interpretation of size effects. Int J Fract 1999;95:299–314. Exadaktylos GE, Vardoulakis I. Microstructure in linear elasticity and scale effects: a reconsideration of basic rock mechanics and rock fracture mechanics. Tectonophysics 2001;335(1–2):81–110. Aifantis EC. Higher order gradients and size effects. In: Carpinteri A, editor. Size-scale effects in the failure mechanisms of materials and structures. New York: Chapman & Hall; 1995. p. 231–42. Efremidis GT, Aifantis EC. The coefficient of geostatic stress: gradient elasticity vs classical elasticity. J Mech Behav Mater 2007;18: 43–54. Tsagrakis I, Konstantinidis A, Aifantis EC. Strain gradient and wavelet interpretation of size effects in yield and strength. Mech Mater 2003;35:733–45. Frantziskonis G, Konstantinidis A, Aifantis EC. Scale-dependent constitutive relations and the role of scale on nominal properties. Eur J Mech A/Solids 2001;20:925–36. Frantziskonis G, Aifantis EC. On the stochastic interpretation of gradient-dependent constitutive equations. Eur J Mech A/Solids 2002;21:589–96. Papanastasiou P. Interpretation of the scale effect in perforation failure. In: Exadaktylos G, Vardoulakis I, editors. Bifurcations, instabilities, degradation in geomechanics. Berlin: Springer; 2007. p. 53–70. Martin CD, Chandler NA. The progressive fracture of Lac du Bonnet granite. Int J Rock Mech Min Sci 1994;31(6):643–59.

1125

[31] Lajtai EZ, Schmidtke RH, Bielus LP. The effect of water on the deformation and fracture time-dependency of a granite. Int J Rock Mech Min Sci Geomech Abstr 1987;24(4):247–55. [32] Read RS, Chandler NA, Dzik EJ. In situ strength criteria for tunnel design in highly-stressed rock masses. Int J Rock Mech Min Sci 1998;35(3):261–78. [33] Cai M, Kaiser PK, Tasaka Y, Maejima T, Morioka H, Minami M. Generalized crack initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int J Rock Mech Min Sci 2004;41(5):833–47. [34] Biolzi L, Meda A, Cardani G. Standaridization of the four-point bending test on natural building stone specimens. In: Characterization of mechanical properties and damage of natural building stones in historical monuments, EC Project, Contract no. SMT4-CT962130, Department of Structural Engineering, Polytechnic Milano, 1999. [35] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci 1983;20(6):249–68. [36] Gerogiannopoulos NG, Brown ET. The critical state concept applied to rock. Int J Rock Mech Min Sci Geomech Abstr 1978;15:1–10. [37] Zimmerman RW. The effect of microcracks on the elastic moduli of brittle materials. J Mater Sci Lett 1985;4:1457–60. [38] Jaeger JC, Cook NGW, Zimmerman RW. Fundamentals of rock mechanics. 4th ed. Oxford: Blackwell; 2007. [39] Mohammad N, Reddish DJ, Stace LR. The relation between in situ and laboratory rock properties used in numerical modeling. Int J Rock Mech Min Sci 1997;34(2):289–97. [40] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci 1997;34(8):1165–86. [41] Ramamurthy T, Arora VK. Strength predictions for jointed rocks in confined and unconfined states. Int J Rock Mech Min Sci 1994; 31(1):9–22. [42] Hoek E, Diederichs MS. Empirical estimation of rock mass modulus. Int J Rock Mech Min Sci 2006;43:203–15. [43] Tatsuoka F. Experimental research for the last 35 years by a geotechnical engineering researcher. In: Proceedings of the seminar ‘‘Deformation and Strength Characteristics of Granular Materials,’’ March 19–April 6, Institute Henri Poincare, Paris, France, 2005. [44] Bhasin R, Barton N, Loset F. Engineering geological investigations and the application of rock mass classification approach in the construction of Norway’s underground Olympic stadium. Eng Geol 1993;35:93–101. [45] Cai M, Kaiser PK, Uno H, Tasaka Y, Minami M. Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the GSI system. Int J Rock Mech Min Sci 2004;41:3–19. [46] Tatsuoka F, Shibuya S. Deformation characteristics of soils and rocks from field and laboratory tests. Report of the Institute of Industrial Science, University of Tokyo, vol. 37, no. 1, serial no. 235, March 1992. [47] Palmstrom A. RMi—a rock mass characterization system for rock engineering purposes. PhD thesis, University of Oslo, Norway, 1995. [48] Itasca Consulting Group Inc. FLAC3D, Fast Lagrangian analysis of continua in 3 dimensions: user’s guide, Minneapolis, 2002