Mjid Noorin Bidgoli, Zhihong Zho, Lnru Jing

aDepartmentofLandandWaterResourcesEngineering,EngineeringGeologyandGeophysicsResearchGroup,RoyalInstituteofTechnology(KTH),Stockholm,Sweden

bDepartmentofGeologicalSciences,StockholmUniversity,Stockholm,Sweden

Numerical evaluation of strength and deformability of fractured rocks

Majid Noorian Bidgolia,?, Zhihong Zhaob, Lanru Jinga

aDepartmentofLandandWaterResourcesEngineering,EngineeringGeologyandGeophysicsResearchGroup,RoyalInstituteofTechnology(KTH),Stockholm,Sweden

bDepartmentofGeologicalSciences,StockholmUniversity,Stockholm,Sweden

A R T I C L E I N F O

Articlehistory:

Received 15 April 2013

Received in revised form 26 August 2013

Accepted 2 September 2013

Strength

Deformability

Fractured rocks

Discrete element methods (DEM)

Failure criteria

Knowledge of the strength and deformability of fractured rocks is important for design, construction and stability evaluation of slopes, foundations and underground excavations in civil and mining engineering. However, laboratory tests of intact rock samples cannot provide information about the strength and deformation behaviors of fractured rock masses that include many fractures of varying sizes, orientations and locations. On the other hand, large-scale in situ tests of fractured rock masses are economically costly and often not practical in reality at present. Therefore, numerical modeling becomes necessary. Numerical predicting using discrete element methods (DEM) is a suitable approach for such modeling because of their advantages of explicit representations of both fractures system geometry and their constitutive behaviors of fractures, besides that of intact rock matrix. In this study, to generically determine the compressive strength of fractured rock masses, a series of numerical experiments were performed on two-dimensional discrete fracture network models based on the realistic geometrical and mechanical data of fracture systems from field mapping. We used the UDEC code and a numerical servo-controlled program for controlling the progressive compressive loading process to avoid sudden violent failure of the models. The two loading conditions applied are similar to the standard laboratory testing for intact rock samples in order to check possible differences caused by such loading conditions. Numerical results show that the strength of fractured rocks increases with the increasing confining pressure, and that deformation behavior of fractured rocks follows elasto-plastic model with a trend of strain hardening. The stresses and strains obtained from these numerical experiments were used to fit the well-known Mohr-Coulomb (MC) and Hoek-Brown (H-B) failure criteria, represented by equivalent material properties defining these two criteria. The results show that both criteria can provide fair estimates of the compressive strengths for all tested numerical models. Parameters of the elastic deformability of fractured models during elastic deformation stages were also evaluated, and represented as equivalent Young’s modulus and Poisson’s ratio as functions of lateral confining pressure. It is the first time that such systematic numerical predicting for strength of fractured rocks was performed considering different loading conditions, with important findings for different behaviors of fractured rock masses, compared with testing intact rock samples under similar loading conditions.

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

1. Introduction

As is well known, natural rock masses consist of intact rock blocks separated by discontinuities such as joints, bedding planes, sheared zones, and faults. The presences of various discontinuities, the inherent complexity of their geometrical parameters, and the difficulties for estimation of their geomechanical and geometrical properties, make it difficult to measure directly mechanical properties of fractured rocks under ordinary laboratory conditions. The main reasons for such difficulties are: (1) the need for testing fractured rocks of large volumes of hundreds of cubic meters, equal to or larger than their representative elementary volume (REV); and (2) the fracture system geometry remains unknown before testing. Clearly laboratory tests are not practical and different techniques are needed for reliable predictions.

Over the years, various attempts have been made to study the strength and deformability of fractured rock masses. Currently, methods available for estimating the strength and deformability of fractured rock masses fall into two broad categories, namely direct and indirect methods. Direct methods are the experimentalmethods in laboratory or in situ tests. Important laboratory tests on intact rock samples have been conducted and comprehensively reported in the literature. However, laboratory tests of intact rock samples cannot provide information of strength and deformation behaviors of fractured rocks, due to the existence of fractures of varying sizes, orientations and locations at larger scales. To obtain realistic results for strength and deformation behaviors of fractured rocks, large volumes of rock containing fractures should be tested at desired stress levels, which is almost impossible to be carried out in conventional laboratory facilities today, but is possible by using direct in situ field tests. However, in situ field tests are usually very difficult to control the initial and boundary (loading) conditions and are time-consuming, and economically costly.

Indirect methods commonly include empirical, analytical and numerical methods. Each of these methods has its own limitations and advantages.

One of the popular and simple indirect methods for estimating strength and deformability of fractured rocks is the empirical methods using the rock mass classification systems such as RQD, RMR, and Q-system, which are based on the engineering experiences obtained from the past projects. In this approach, rock mass properties are linked to a representative rock mass classification index that reflects the overall rock mass quality. The main shortcoming of this approach is that it lacks a proper mathematical platform to establish constitutive models and the associated properties of the fractured rocks, so that the second law of thermodynamics should not be violated, since complex properties of a rock mass cannot be satisfactorily represented in this method for establishing constitutive models and their associated properties concerned quantitatively with a proper mathematical logic.

Analytical methods are very useful in geomechanics because they provide results that can highlight impacts of the most important issues or variables that determine the solution of a problem. Analytical methods attempt to calculate strength and deformability of fractured rocks from the strength and deformation properties of the discontinuities and of the intact rock matrix, but are applicable only with simple and regular fracture system geometry. These limitations make this approach impossible for fractured rocks containing complex fracture systems.

Numerical methods can be used to calculate strength and deformability of fractured rocks with more flexibility, by representing different mechanical and geometric features of the fractures and the intact rock matrices. With almost daily improvements of efficiencies of numerical solution methods and increase of computing power, numerical modeling methods have been developed to estimate the strength and deformability of fractured rocks by using various discrete and continuum modeling methods. The FEM (finite element method) is the most widely applied numerical method for rock engineering problems (Pouya and Ghoreychi, 2001; Sitharam, 2009), besides the DEM (discrete element method). Since the FEM models are based on an overall continuum material assumption, effective and reliable considerations of effects of a large number of fractures of different sizes, orientations and behaviors are still difficult.

The DEM was introduced by Cundall (1971) and further developed by Cundall and co-workers (Lemos et al., 1985; Lorig et al., 1986; Cundall, 1988; Hart et al., 1988). A comprehensive presentation of the DEM can be found in Jing and Stephansson (2007). The method is a powerful technique to perform stress analyses for blocky rock masses formed by fractures, since its advantage of explicit representations of both the fracture system geometry and constitutive behaviors of fractures and intact rock matrix. Therefore, both deterministic and stochastic approaches can be applied for such evaluations. Since fracture systems in rock masses are geometrically complex and largely hidden in subsurface without being exposed, a large number of discrete fracture network realizations, based on the probabilistic distribution functions of geometrical parameters, are needed as the geometric models for statistical numerical modeling of fractured rocks (Priest, 1993). Some of recent publications using DEM were given by Park et al. (2006), Christianson et al. (2006), Kim et al. (2007), Noel and Archambault (2007), Zhang et al. (2007), Cundall et al. (2008), Singh and Singh (2008), Esmaieli et al. (2010), Wu and Kulatilake (2012), and Khani et al. (2013).

Min and Jing (2003) developed a numerical modeling approach of scientific originality for evaluating equivalent hydro-mechanical properties of fractured rocks, using DEM. However, for mechanical properties, the study presented was limited to equivalent elastic properties, not strength, of the fractured rock concerned. Baghbanan (2008) tried to attack this problem by adopting an M-C type of equivalent behavior, but no journal publications appeared in the literature on this subject at that time. The paper aims to solve this problem by extending the research considering effects of different loading conditions that may be considered in future laboratory or field testing of fractured rocks of large volumes with realistic fracture system geometry. Since most of the civil and mining engineering projects are built in or on rock masses, reliable estimation and good understanding of the strength and deformability of fractured rocks remain as a challenging issue for safe and economical design, construction and stability evaluation of slopes, foundations and underground excavations.

The first aim of this research is to extend the original numerical platform (e.g. Min and Jing, 2003; Baghbanan, 2008) for predicting strength and deformability of fractured crystalline rocks, in a general sense, due to its ability for more realistic and complex geometrical representation of fractured rocks. The second aim is to test fractured rock behaviors with different loading conditions, i.e. under controlled axial load and axial velocity, in a similar way as testing intact rock samples in laboratories. The reason is to understand any differences when testing large volumes of fractured rocks with fracture and intact rock matrix of different constitutive behaviors, since such tests have not been performed yet.

There are several empirical failure criteria developed for representing strength of rock masses. The M-C failure criterion and the H-B failure criterion (Hoek and Brown, 1980; Hoek, 1983; Hoek et al., 1992, 2002) are the two most commonly accepted failure criteria in the international rock mechanics community, with the H-B criterion widely applied to hard rocks such as granites. The third aim is to see whether these two popular failure criteria may or may not yield very different estimations of compressive strength, using the stress and strain results of numerical models of fractured rock concerned.

The universal distinct element code, UDEC (Itasca Consulting Group Inc., 2004), was used to perform numerical uniaxial and biaxial compressive tests on fractured rock models containing a large number of fractures of varying sizes, created using stochastic discrete fracture network (DFN) method and realistic fracture system information. A few deterministic fracture system models extracted from randomly generated fracture system realizations were adopted, and statistical models using Monte Carlo simulations with multiple fracture system realizations will be studied at the next step, and the results will be published in due course.

For simplicity, the term ‘fracture’ and ‘fractures’ are adopted as the general term for all types of discontinuities of rocks such as faults, joints, fracture zones, etc., unless specified separately.

Fig. 1. Flowchart for a numerical stress–strain analysis processes in a fractured rock.

2. Methodology

In this study, a systematic investigation was conducted using numerical experiments of typical laboratory compression tests to determine, generically, the compressive strength and deformation parameters of the fractured rock, as equivalent properties at its REV size. During loading on a numerical model of the fractured rock concerned, both rock matrix and fractures will deform or be displaced, governed by the equations of motions of the rock blocks and constitutive models, material parameters for rock matrix and fractures, and the initial and boundary conditions. Fig. 1 shows a flowchart used in this study for a numerical stress–strain analysis process in a fractured rock.

2.1.DFNmodelfornumericalexperiments

At the first stage of this research, a few DFN models were generated to represent the fractured rock masses and after that, the generated geometry was used to create DEM model for the numerical experiments using the UDEC code. It should be noted that it is generally accepted that the sizes of DFN models must not be less than its REV of the models concerned. REV is de fined as the minimum volume (or a range) beyond which the characteristics of the domain remain basically constant. Therefore, when the sizes of DFN models are not less than their REVs, the equivalent properties will become scale-independent. Min and Jing (2003) and Min et al. (2004) conducted numerical studies to establish elastic compliance tensor and permeability tensor for fractured rock masses, by investigating the scale-dependent equivalent permeability of fractured rock at the Sella field site, Cumbria, England. Their results showed that an acceptable REV scale is above 5 m × 5 m for the fracture systems with constant apertures, for both elastic compliance tensor and permeability tensor of the concerned fractured rock as an equivalent continuum.

In the investigation of strength and deformability of fractured rock, the same field data from Sella field site, as reported in Min and Jing (2003), were used. Three square DEM models of fracture systems were generated with side length of 2 m × 2 m, 5 m × 5 m, and 10 m × 10 m, respectively, as extracted from the center of an original parent model of fracture system, based on the same fracture system model data as was used in Min and Jing (2003). The geometric parameters for generating fracture network realizations were based on the field mapping results of a site characterization at the Sellafield area, undertaken by the United Kingdom Nirex Limited. The basic information of the identified four sets of fractures is shown in Table 1.

Fig. 2. A DEM model with size of 2 m × 2 m before (the upper) and after (the lower) the fracture system regularization.

Then the DFN models were used to generate DEM models with internal discretization of finite difference elements, for stressdeformation analyses. Before performing the analyses, the DFN models were regularized by deleting the isolated fractures and dead-ends, so that the resultant fractures were all connected and each fracture contributes to form two and just two opposing surfaces on two adjacent blocks. Fig. 2 shows a DEM model of 2 m × 2 m in size before and after the fracture system regulation as an example. The DEM model with size of 2 m × 2 m was used only for demonstrating the different results obtained when model sizes are less than 5 m × 5 m, and the DEM model with size of 10 m × 10 m was used to ensure the validity of the 5 m × 5 m REV (see Section 3). The REV of the fractured rock concerned was already established in Min and Jing (2003) as 5 m × 5 m so that there was no need for further study on scale effects.

Table 1 Fracture parameters used for DFN generation (Min and Jing, 2003).

2.2.DEMmodelfornumericalexperiments

A number of general and important assumptions related to the set-up of the DEM models about rock matrix and fractures, as adopted in Min and Jing (2003) and Baghbanan (2008), were adopted in this study:

(1) The numerical model was de fined in a two-dimensional (2D) space for a generic study.

(2) Simulations were performed under quasi-static plane strain conditions for deformation and stress analyses, without considering effects of gravity.

(3) Fractured rock was a hard rock mass, containing rock matrix and fracture, without considering strain-softening.

(4) Rock matrix was a linear, isotropic, homo-geneous, elastic, and impermeable material.

(5) The fractures follow an ideal elasto-plastic behavior of an M-C model in the shear direction and a hyperbolic behavior (Bandis’Law) in the normal direction.

(6) The initial aperture of fractures (without stress) was a constant.

(7) Coupled hydro-mechanical effects on the fractures were neglected in the current study.

The above assumptions are based on measured data from site investigations, and are necessary for a numerical prediction for strength and deformability of fractured rocks, without attention of application to site-speci fic case studies.

The basic information about the intact rock, the granite matrix, and mechanical properties of fractures that were used for modeling in UDEC is shown in Table 2. This information was based on the laboratory test results reported in Sella field site investigation, which was used in Min and Jing (2003).

2.3.Modelingprocedure

In this study, similar to the standard compression test of small intact rock samples in laboratory, a series of numerical experiments, namely uniaxial and biaxial compression tests, were performed, following the modeling procedure (see Fig. 1) on three DEM models of varying sizes to determine the strength and deformability of fractured rocks.

The applied axial compressive stress loading condition was similar to that of the standard confined compression tests on axisymmetric intact rock samples, by varying confining pressures and axial loads. Fig. 3 shows the typical physical set-up and boundary conditions of uniaxial (Fig. 3a) and biaxial compression (Fig. 3b) tests, respectively. For both uniaxial and biaxial compression tests, the bottom of the DEM models was fixed in they-direction and an axial load (σy) was applied on the top of the DEM model. Varying confining pressure (σx) was applied laterally on the two vertical boundary surfaces of the model, as in the biaxial compression tests. For the uniaxial compression tests, the two vertical sides of the DEM model were kept as free surfaces. The DEM models were loaded sequentially with a constant and very small axial load increment (Δσy), equal to 0.05 MPa, in every loading step of calculation in the vertical direction, the same as a conventional uniaxial or triaxial loading tests on intact rock samples.

The axial stress loading process was controlled by a velocity monitoring scheme during simulation. The velocities (in the bothx-andy-directions) at a number of carefully specified monitoring points were checked to ensure that they become zero or very close to zero at the end of every loading step so that a quasi-static state of equilibrium of the model was reached under the applied boundary conditions, since the simulated tests should be quasi-static tests for generating static behaviors of the models. Six parallel sampling lines within each model were placed in bothx- andy-directions, with the same distance in between. Therefore, thirty-six points were defined at intersections of the horizontal and vertical monitoring lines. These points plus one point at the center of the DEM model were the monitoring points in this study. Fig. 4 shows a DEM model with a size of 10 m × 10 m and positions of monitoring points into the DEM model for using velocity monitoring technique during loading compression tests.

Fig. 3. Typical set-ups and boundary conditions for numerical experiments. (a) Set-up for uniaxial compression tests, (b) set-up for biaxial compression tests.

Vertical and horizontal velocities (iny-velocity andxdirections), vertical and horizontal displacements (y- andxdisplacements), normal and shear stresses (σyy, σxxand τxy) weremonitored at all monitoring points at each loading step during the uniaxial and biaxial compression tests. The mean values of normal stresses and strains in thex- andy-directions were calculated at the end of each loading step. The average stresses and strains were computed by taking the average values obtained from the monitoring points by using the FISH algorithm, the programming language embedded within the UDEC code.

Table 2 Mechanical properties of intact rock and fractures (Min and Jing, 2003).

Fig. 5 shows curves of velocity versus time inx- andy-directions at 6 selected monitoring points located on two horizontal and vertical lines within a DEM model during a few loading compression tests. It can be observed that values of velocities at the defined monitoring points (in the bothx- andy-directions) became very close to zero at the end of every loading step. Using the same velocity monitoring grid, the deformation and stress of each DEM model were evaluated in order to calculate the average stress and strain values of the tested model, which were then used to evaluate the equivalent strength and deformability parameters when different strength criteria were adopted.

Fig. 4. DEM model with size of 10 m × 10 m and position of monitoring points into model. (a) Fracture system model after regularization, (b) locations and numbering of the monitoring points.

In order to keep a servo-controlled loading condition, a new FISH program was developed and inserted in the UDEC model to simulate a standard servo-controlled test similar to the standard servo-controlled tests of small intact rock samples in laboratory, to minimize the in fluence of inertial effects on the response of the model, by setting the upper and lower limits for unbalanced forces. Cyclic loading rate was kept in a range of maximum and minimum unbalanced forces in UDEC program to avoid sudden (violent) failure of the DEM models during cycles of uniaxial and biaxial compression tests.

It should be noted that the equivalent strength and deformability of the fractured rock, as an equivalent continuum, were the concern of research, not the complete constitutive model of the fractured rock as an equivalent continuum under any stress paths. Therefore, the loading needs to be stopped when the peak strength of the model was reached, without model collapse or appearance of very large shear displacements along the fractures or large block motion, which will make the equivalent continuum assumption of the fractured rock invalid, and the homogenization (averaging) for equivalent parameter evaluation could not be applied.

3. Stress–strain behaviors of the fractured rock under compression

3.1.Stress–strainbehavioroffracturedrockunderaxialstress loadingcondition

Fig. 6 shows the results of stress–strain behaviors of the fractured rock mass models under uniaxial compression tests without con fining pressure, as the curves of axial stresses versus axial strains for DEM models with varying sizes. Figs. 7–9 show the numerical test results with different con fining pressures of 0.5 MPa, 1 MPa, 1.5 MPa, 2 MPa, 2.5 MPa and 3 MPa, respectively, as the curves of axial stresses versus axial strains in the DEM models with sizes of 2 m × 2 m, 5 m × 5 m and 10 m × 10 m, respectively. These curves were used to evaluate strength behaviors of the fractured rocks after the models reached their peak strengths.

It can be seen from Figs. 7–9 that the DEM models deform linearly and elastically at axial stresses below the yield strength, depending on the con fining pressure. Further compression leads to inelastic deformation up to the peak strength. With increase of confining pressure, the strength of the DEM models increases and the stress–strain curves follow an elasto-plastic behavior with a strain hardening trend. Also, behaviors of the DEM models change with increase in the model size, but the change becomes insigni ficant between models of size of 5 m × 5 m and 10 m × 10 m. Therefore, the model of the established REV size of 5 m × 5 m is adequate for evaluating the equivalent strength and deformability of fractured rock concerned.

Fig. 5. Curves of velocity versus time in x- and y-directions at the six monitoring points. The numbers of the monitoring points and their locations in the 10 m × 10 m are shown in Fig. 4b.

Fig. 6. Axial stress versus axial strain curves for DEM model with varying sizes under uniaxial compression conditions, without confining pressure.

Fig. 7. Axial stress versus axial strain curves for DEM model with size of 2 m × 2 m under different confining pressures.

Fig. 8. Axial stress versus axial strain curves for DEM model with size of 5 m × 5 m under different confining pressures.

Fig. 9. Axial stress versus axial strain curves for DEM model with size of 10 m × 10 m under different confining pressures.

Fig. 10. Axial stress versus axial strain curves for DEM model with size of 5 m × 5 m under different confining pressures, with constant velocity condition at the top surface.

Fig. 11. Axial stress versus axial strain curves for DEM model with size length of 10 m × 10 m under different confining pressures, with constant velocity condition at the top surface.

Fig. 12. Comparison of axial stress versus axial strain curves for DEM model with size of 5 m × 5 m under different confining pressures, between using constant normal velocity loading condition (solid lines) and the constant normal stress loading condition (dash lines).

Fig. 13. Comparison of axial stress versus axial strain curves for DEM model with size length of 10 m × 10 m under different confining pressures, between using constant normal velocity loading condition (solid lines) and the constant normal stress loading condition (dash lines).

3.2.Effectofloadingconditionsonstress–strainbehaviorsof fracturedrockmass

Besides the axial compressive stress loading condition, another common and popular method for numerical experiments is the constant velocity boundary condition. In this context, for studying the effect of different axial loading conditions on stress–strain behaviors of fractured rock masses, numerical experiments were also performed by applying a constant velocity in they-direction at the top boundary of the models. The downward velocity loading conditions applied were similar to the previously mentioned axial stress boundary condition. The bottom of the DEM model was fixed in they-direction and the numerical experiments were performed by varying the confining pressure.

Figs. 10–11 show the results of average axial stress versus axial strain curves with different confining pressures of 0.5 MPa, 1 MPa,1.5 MPa, 2 MPa, 2.5 MPa and 3 MPa, with a constant velocity boundary condition, for two DEM models of 5 m × 5 m and 10 m × 10 m, respectively. Figs. 12–13 compare the obtained stress–strain curves using the constant velocity boundary condition (solid line) and the axial stress boundary condition (dash line) with different con fining pressures.

Fig. 14. Normal stress versus shear displacement curves of direct shear tests of rock fractures under different system stiffness (Skinas et al., 1990).

Numerical results show a higher average axial stress under constant velocity test condition than that under constant axial stress condition. The main reason is that the higher normal stress of fractures were induced during shear under constant velocity conditions, due to the 5?dilation angle (see Table 2), a phenomenon observed in many shear tests of rock joints under constant strain conditions that are equivalent to the constant velocity condition as applied in this paper, and reported widely in the literature (Fig. 14). The amount of difference is decided by the values of dilation angle. Minor difference may occur when the dilation angle is very small. There may be other reasons, such as block interlocking and stress concentration at fracture intersections, which may also contribute to such differences locally, but play a less important role compared with shear dilation effects on normal stress of fractures.

4. Strength of fractured rock

The M-C and H-B failure criteria were selected to be fitted for representing the equivalent strength of the fractured rock concerned, due to their wide acceptance in the international rock mechanics community.

4.1.Thefailurecriteria

The M-C failure criterion is an empirical linear failure criterion that has been adopted for different rocks and soils. It can be expressed in a functional relation as

where τmaxis the shear strength, σnis the normal stress,cis the cohesion, and ? is the internal friction angel. This criterion can be expressed in terms of principal stresses as

Fig. 15. Strength curves for DEM models in the normalized principal stress space. (a) With model size of 5 m × 5 m, (b) with model size of 10 m × 10 m.

where σ1is the major principal stress at failure or elastic strength, and σ3is the minor principal stress or confining pressure.

The M-C failure criterion can be applied for both intact rocks and rock masses, with the parametercand ? changes representing effects of fracture and intact rock properties on the overall equivalent strength of the fractured rock mass concerned.

The H-B failure criterion is an empirical nonlinear failure criterion that is proposed for failure of intact rocks and rock masses. It can be expressed in terms of principal stresses (Hoek and Brown, 1980) aswhere σciis the uniaxial compressive strength of the intact rock;mandsare the material constants, wheres= 1 for intact rocks.

Fig. 16. Estimated values of equivalent Young’s modulus for DEM model with varying size under different confining pressures (Young’s moduli at zero confining pressure are 43 and 72 MPa for the two models with sizes of 5 m × 5 m and 10 m × 10 m, respectively).

Fig. 17. Estimated equivalent values of Poisson’s ratio for DEM model with varying sizes under different confining pressures.

Table 3 Equivalent material parameters of Mohr-Coulomb and Hoek-Brown failure envelops.

The strength parameters of cohesion (c) and internal friction angle (?) of M-C criterion, and themandsfor H-B criterion were used for fitting the two respective strength criteria curves as de fined by Eqs. (2) and (3), with results obtained from models using constant axial stress conditions.

4.2.StrengthcurvefittingwithM-CandH-Bcriteria

The curve fitting with M-C and H-B failure criteria was illustrated in Fig. 15 for the DEM models of 5 m × 5 m and 10 m × 10 m in size, under different con fining pressures. Both the M-C and H-B strength envelops made acceptable fitting to the numerical data, with insigni ficant difference between them, despite the fact that the M-C criterion is a linear and H-B criterion is a nonlinear one. The data used for generating these curves are obtained from the normal stress conditions.

The equivalent material parameters for the two criteria, derived from the fitting to the strength criteria, are given in Table 3. The results clearly show that the differences between strength parameter values of DEM models are basically minor between model sizes of 5 m × 5 m and 10 m × 10 m. The difference between the correlation coef ficient root values (R) of two failure envelops is also insigni ficant.

4.3.Estimationofdeformabilityparametersofthefracturedrock concerned

In addition of the estimated strength envelops, the deformability parameters of fractured rock concerned, as the equivalent Young’s modulus and Poisson’s ratio, using the DEM models of different sizes and under different con fining pressures, were calculated, as illustrated in Figs. 16 and 17, respectively. The Young’s modulus was calculated as the averaged local slope of the stress–strain curves of the DEM models during the stage of elastic deformation, and the Poisson’s ratio was calculated as the ratio of the mean transverse strain to the mean axial strain of the DEM models. The data used for generating these curves are obtained from the normal stress conditions.

As we can see from Fig. 16, the Young’s modulus of the fractured rock increases gradually with increase of con fining pressure, but the effect of DEM model sizes is not very signi ficant. Fig. 17 shows that the Poisson’s ratio of the fractured rock decreases gradually with increase of con fining pressure and size of DEM models. Generally, the magnitude of Young’s modulus for fractured rocks is less than the Young’s modulus for intact rock (Fig. 16), but the magnitude of Poisson’s ratio for fractured rocks is much larger than that for intact rock (Fig. 17). The general trends of the curves are the same and converge when the con fining pressure reaches 3.0 MPa, but with a considerable difference between the two models of different sizes after the con fining pressure is larger than 1.5 MPa. The obtained results of larger Poisson’s ratio indicated that fractures affected the deformability of the rock mass much more signi ficantly that of strength, so that care should be taken when developing constitutive models of fractured rocks as equivalent continua.

5. Discussions and conclusions

A systematic 2D numerical procedure to predict strength behavior and deformability parameters of fractured rocks, by using the DEM, was developed in this study for the first time, as an extension to the research performed in Min and Jing (2003) and in Baghbanan (2008). 2D numerical experiments were performed on three geometric models of a fractured rock with varying sizes and with realistic fracture system geometry data from an in situ fracture mapping. Results obtained from these numerical experiments were used to fit the M-C and H-B failure criteria, and to calculate the deformability parameters, the Young’s modulus and Poisson’s ratio, respectively. It is noted that, due to the lack of measured data support from laboratory or in situ experiments of testing volumes not less than the REV sizes of the granite rocks at the Sella field site, this study has been performed in a generic form in nature and the results have only conceptual values. The main conclusions are summarized as follows:

(1) The DFN-DEM is a suitable and flexible numerical approach to predict the behaviors and properties of fractured rocks that cannot be obtained by conventional laboratory tests using small intact rock samples. This method provides an important extension to the comprehensive modeling procedure as developed in Min and Jing (2003) and Baghbanan (2008), compared with empirical models of rock mass classi fication or analytical solutions considering only regular system geometry, despite the fact that the procedure requires much more computing time compared with that used by the empirical and analytical methods.

(2) The axial stress loading condition with servo-controlled constant axial strain loading condition generated different stress–strain behaviors from that under constant normal stress loading conditions, with the same testing model geometry and size, due to the fact that the effect of shear dilation induced increase of normal stresses of fractured under shearing. Whether such difference may also be caused by the effects of stress concentration at fracture intersections or block rotations remains an important issue for further investigations.

(3) The results show that the model size (or scale) has a signi ficant in fluence on the strength and deformation behaviors of fractured rocks if the model size is less that the properly derived REV size of the rock concerned. In this paper, the strength of fractured rocks decreases slightly with the model size up to the established REV size, after that the changes become insigni ficant. Also, deformability parameters of fractured rocks, including the Young’s modulus and Poisson’s ratio, change signi ficantly with con fining pressures. These findings, rarely reported in the literature, have signi ficant impact on developing constitutive models of fractured rocks as equivalent continua, since the REV size, strength and deformability parameters are the three issues that must be readily understood for constitutive model development for fractured rocks whose properties have signi ficant dependence on understanding of stress and size effect.

(4) The results show that the strength and deformation behaviors of fractured rocks are nonlinear over the concerned range of stress and are also dependent upon con fining pressure. With increase of con fining pressure, the strength of fractured rocks increases and deformation behavior of fractured rocks follows an elasto-plastic model with a strain hardening trend. The results show that the testing volume, loading conditions (both axial stress and velocity loading conditions) and adequate quantitative knowledge on fracture system geometry and their mechanical behaviors play a signi ficant role for designing future physical tests for estimating the strength and deformability of fractured rocks, very different and much more challenging compared with testing intact rock samples. The mechanical behavior of the individual fractures plays a signi ficant role in understanding the strength and deformability of fractured rocks, besides the fracture system geometry, as demonstrated by the effect of dilation angle of the research.

(5) Our model did not show strain-softening, except in one case for a model less than REV size (2 m × 2 m) at an early stage of loading (Fig. 6), since we need to stop loading when the peak strength of the model was reached in order to maintain a physical basis of equivalent continuum assumption of the rock mass concerned. Strain softening may occur with continued loading.

(6) Both the M-C and H-B criteria give a fair estimate of the compressive strength of the rock concerned for almost all cases tested in this study. On the other hand, the H-B criterion is in essence a nonlinear failure envelope and is more flexible for modeling different fracture system geometries and stress conditions.

The scienti fic originality of this paper is its extension of the mathematical platform established in Min and Jing (2003) for logical representation of complex and realistic fracture system geometry, from in situ fracture mapping, for estimating strength and deformability of fractured rocks. Similar early research used more regular or much simpli fied fracture systems, and the effects of such simpli fications cannot be properly estimated. However, there are some outstanding issues that remain to be addressed in future.

(1) The study presented is based on an assumption that the initial aperture of fractures is a constant. So, more modeling is needed to perform numerical experiments when initial aperture of fractures is not constant but correlated to fracture size.

(2) In this research, the models have not shown strain-softening, but it may occur after continued loading or change of fracture system or shape and geometry of the testing volume. By continued loading, strain-softening may occur, but it is not required at this stage of the research, since deriving a meso-scale comprehensive constitutive model of the fractured rocks is a future work. However, the subject is new, first tested by numerical modeling in this research, and needs continued investigations.

(3) Partial cracking and complete crushing of rock blocks during loading processes were not considered. This issue may affect evaluating the equivalent strength and deformability of fractured rocks to a certain extent, but may not be an important factor since the dominating factor for strength and deformation of the fractured rocks is the displacements of fractures, which is the main mechanism of the energy dissipation according to the first principle of energy minimization. In addition, this simplifying assumption was needed by the current version of the UDEC code that does not have the ability to consider block cracking, and this is an issue for future work.

(4) Coupled hydro-mechanical effects on the fractures were neglected for evaluation of strength and deformation behaviors of fractured rocks at this stage of research. Further studies are needed to study water pressure effects on strength and deformation behaviors of fractured rocks.

(5) Stochastic analysis using multi-fracture system realizations needs to be performed for a more comprehensive understanding of the uncertainty of the predicting of fractured rock behaviors.

(6) 3D investigations are necessary to eliminate the limitations caused by the assumption of 2D space under plane strain loading conditions.

Baghbanan A. Scale and stress effects on hydro-mechanical properties of fractured rocks masses. Stockholm: Ph. D. Thesis. Royal Institute of Technology; 2008.

Christianson MC, Board MP, Rigby DB. UDEC simulation of triaxial testing of lithophysal tuff. In: Proceedings of the 41st U.S. Symposium on Rock Mechanics (USRMS) and 50th Anniversary of the U.S. Rock Mechanics Symposium. Golden Rocks 2006. Golden, Colorado: American Rock Mechanics Association, 2006: Paper No. 06-968; 2006.

Cundall PA. A computer model for simulating progressive large-scale movements in blocky rock system. In: Proceedings of International Symposium on Rock Mechanics. ISRM; 1971. p. 128–32.

Cundall PA. Formulation of a three-dimensional distinct element model. Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences 1988;25(3):107–16.

Cundall PA, Pierce ME, Mas Ivars D. Quantifying the size effect of rock mass strength. In: Proceedings of the 1st Southern Hemisphere International Rock Mechanics Symposium. Australian Centre for Geomechanics (ACG); 2008. p. 3–15.

Esmaieli K, Hadjigeorgiou J, Grenon M. Estimating geometrical and mechanical REV based on synthetic rock mass models at Brunswick Mine. International Journal of Rock Mechanics and Mining Sciences 2010;47(6):915–26.

Hart R, Cundall PA, Lemos J. Formulation of a three dimensional distinct element model. Part II. Mechanical calculation for motion and interaction of a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences 1988;25(3):117–25.

Hoek E. Strength of jointed rock masses. Geotechnique 1983;33(3):187–223.

Hoek E, Brown ET. Empirical strength criterion for rock masses. Journal of Geotechnical Engineering Division, ASCE 1980;106(9):1013–35.

Hoek E, Wood D, Shah S. A modified Hoek-Brown criterion for jointed rock masses. In: Proceedings of Rock Characterization Symposium; 1992. p. 209–14.

Hoek E, Carranza TCT, Corkum B. Hoek-Brown failure criterion-2002 edition. In: Proceedings of the 5th North American Rock Mechanics Society Meeting (NARMS-TAC); 2002. p. 267–73.

Itasca Consulting Group Inc. UDEC 4.0 user’s guide. Minneapolis, MN: Itasca Consulting Group Inc; 2004.

Jing L, Stephansson O. Fundamentals of discrete element methods for rock engineering: theory and applications (developments in geotechnical engineering). Amsterdam, Netherlands: Elsevier Science BV; 2007.

Khani A, Baghbanan A, Norouzi S, Hashemolhosseini H. Effects of fracture geometry and stress on the strength of a fractured rock mass. International Journal of Rock Mechanics and Mining Sciences 2013;60(1):345–52.

Kim BH, Cai M, Kaiser PK, Yang HS. Rock mass strength with non-persistent joints. In: Proceedings of the 1st Canada-US Rock Mechanics Symposium. Taylor and Francis Group; 2007. p. 241–8.

Lemos JV, Halt RD, Cundall PA. A generalized distinct element program for modeling jointed rock mass. In: Proceedings of the International Symposium on Fundamentals of Rock Joints. Centek; 1985. p. 335–43.

Lorig LJ, Brady BHG, Cundall PA. Hybrid distinct element-boundary element analysis of jointed rock. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1986;23(4):303–12.

Min KB, Jing L. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. International Journal of Rock Mechanics and Mining Sciences 2003;40(6):795–816.

Min KB, Jing L, Stephansson O. Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: method and application to the field data from Sellafield, UK. Hydrogeology Journal 2004;12(5): 497–510.

Noel JF, Archambault G. Boundary conditions effects on the formation of kink bands in fractured rock masses. In: Proceedings of the 1st Canada-US Rock Mechanics Symposium. Taylor and Francis Group; 2007. p. 27–31.

Park ES, Jeon YS, Ryu CH, Martin CD, Christiansson R. Numerical approach for a rock mechanics descriptive model. Tunnelling and Underground Space Technology 2006;21(3-4):446–7.

Pouya A, Ghoreychi M. Determination of rock mass strength properties by homogenization. International Journal for Numerical and Analytical Methods in Geomechanics 2001;25(13):1285–303.

Priest SD. Discontinuity analysis for rock engineering. London: Chapman and Hall; 1993.

Singh M, Singh B. Laboratory and numerical modeling of a jointed rock mass. In: Proceedings of the 12th International Conference of International Association for Computer Methods and Advances in Geomechanics (IACMAG). India Institute of Technology, Bombay; 2008. p. 1373–80.

Sitharam TG. Equivalent continuum analyses of jointed rock mass: some case studies. International Journal of the JCRM 2009;5(1):39–51.

Skinas CA, Bandis SC, Demiris CA. Experimental investigations and modeling of rock joint behavior under constant stiffness. In: Proceedings of the International Symposium on Rock Joints (ISRM). A.A. Balkema; 1990. p. 301–8.

Wu Q, Kulatilake PHSW. REV and its properties on fracture system and mechanical properties, and an orthotropic constitutive model for a jointed rock mass in a dam site in China. International Journal of Compute and Geotechnics 2012;43:124–42.

Zhang ZG, Qiao CS, Li X. Study on the strength of jointed rock mass. Journal of Key Engineering Materials 2007;353(1):381–4.

?Corresponding author. Tel.: +46 8 790 8661.

E-mail address: mnoorian@kth.se (M.N. Bidgoli).

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

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

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