Numerical modeling of rock mechanical behavior and fracture propagation by a new bond contact model

International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

Numerical modeling of rock mechanical behavior and fracture propagation by a new bond contact model Mingjing Jiang a,b,n, He Chen a,b, Giovanni B. Crosta c a

Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai, 200092, China Key laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai, 200092, China c Department of Earth and Environmental Sciences, Università degli Studi di Milano Bicocca, Milano, 20126, Italy b

art ic l e i nf o

a b s t r a c t

Article history: Received 15 July 2013 Received in revised form 23 March 2015 Accepted 25 March 2015

This paper presents a numerical investigation into mechanical behavior and rock fracture using the distinct element method (DEM). Based on a series of laboratory tests on the bonded granules idealized by two glued aluminum rods, a normal force dependent bond contact model was proposed and implemented into a two-dimensional (2D) DEM code. This 2D DEM code was used to carry out a series of numerical simulations, including uniaxial and biaxial compression tests, direct tension and Brazilian tension tests, whose results were compared with experimental data to calibrate the proposed model and identify the microscopic deformability and strength parameters. In addition, the validated model was then used to simulate crack propagation and rock fracture in single-flawed and double-flawed samples with different flaw inclinations, and the simulations were then compared with experimental observations. The numerical results demonstrate that the proposed bond model incorporated into the distinct element method is able to capture the main mechanical behaviors of crystalline rocks (Lac du Bonnet granite and Hwangdeung granite). The limitations associated to a low strength envelope angle and high ratio between tensile strength and compressive strength frequently formed in DEM simulations of rock behaviors are solved. In addition, five distinct stages of the stress–strain curve and the marked stresses can be clearly identified by lateral strain, volumetric response and the cumulative number of bond failures. The mechanical behavior and rock fracture patterns obtained from DEM simulations are in good agreement with experimental observations and the model captures both the geometrical control and the general strength value for crystalline rocks. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Distinct element method Bond contact model Crack propagation Crack initiation stress Crystalline rock Cemented particles

1. Introduction Rock is one of most important and complex materials in geotechnical engineering. A correct design of rock slopes, mines and nuclear waste repositories requires a series of laboratory tests and numerical models of rock behavior under complex loading conditions [1,2]. Numerical analysis methods are useful tools that can support the design and implementation of geomechanical engineering. However, as a result of a variety of geological processes, rock contains a large number of joints, fissures, weak surfaces and faults. These pre-existing discontinuities or cracks make of the rock a discontinuous, inhomogeneous, anisotropic and non-elastic compound structure [3] bringing some major difficulties in numerical study. Hence, a comprehensive study of numerical simulation on

n Correspondence to: College of Civil Engineering, Tongji University, 200092 Shanghai, China. Tel.: þ86 21 65980238; fax: þ 86 21 65985210. E-mail address: [email protected] (M. Jiang).

http://dx.doi.org/10.1016/j.ijrmms.2015.03.031 1365-1609/& 2015 Elsevier Ltd. All rights reserved.

rock mechanical behavior is essential due to the complexity of related problems. Some fundamental researches have been carried out on the numerical simulation of rocks in order to investigate crack propagation and rock fracture. Reyes and Einstein [4] developed an analytical model for simulating secondary crack coalescence formulated by combining a smeared crack/damage mechanics approach with a strain based failure criterion. Bobet and Einstein [5] studied crack propagation and coalescence in rock specimens containing two inclined flaws that were either open or closed by employing a new stress-based crack initiation criterion. Scavia [6] and Vasarhelyi and Bobet [7] studied the mechanical behavior of rocks and the crack growth between two bridged flaws by the displacement discontinuity method (DDM). Tang et al. [8,9] and Wong et al. [10] carried out a series of numerical simulations on models containing pre-existing flaws under uniaxial compressive loading in order to investigate the failure mechanism and fracture mode using a rock failure process analysis code (RFPA2D). In addition, the maximum tangential stress theory [11], the maximum energy release rate theory [12], the maximum energy

176

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

were simulated. DEM macroscopic responses were calibrated by the published experimental data [42,43] in order to identify a set of microscopic material parameters. By means of calibrated parameters, crack propagation and rock fracture in single-flawed and doubleflawed DEM samples under uniaxial compressive loading were modeled and compared against experimental observations [44].

2. A bond contact model for rock 2.1. Experimental set-up Fig. 1 presents the analogous cementation sample used in the laboratory tests consisting of a pair of aluminum rods glued together by the rock-like material in the inter-particle voids, i.e. cement, with a bond thickness of 0.6 mm and a width of 3 mm. The aluminum rods have a 12 mm diameter and are 50 mm long. This assembled sample was used to carry out a series of tests under various loading conditions, as illustrated in Fig. 2. Five different types of test were designed to characterize the mechanical behavior of the cemented samples. The following loading paths are referred to as tension, compression, shear under different normal forces, rolling under different forces and shear-rolling tests under different normal forces. 2.2. Mechanical response The mechanical response of the bond contact model was established firstly through reasonable theoretical derivation [45–47], and then it was verified by a series of laboratory tests on the cemented samples [40]. The mechanical response is characterized by three stiffness coefficients kn, ks, km, three bond strength coefficients Rt, Rs, Rr, and the related residual strength Rsr for friction and Rrr for rolling, as illustrated in Fig. 3. Fig. 3(a) presents the normal mechanical response which is composed of tension and compression directions. In tension, the normal contact force Fn between two particles is a product of normal stiffness coefficient kn and tension displacement un before its value reaches the normal bond strength Rt. When the force arrives at the normal bond strength Rt, the bond is broken and the normal force abruptly reduces to zero. In compression, the contact force is always a product of kn and compression displacement, un. Fig. 3(b) presents the tangential force-displacement relationship. The shear force increases linearly with the increasing of tangential displacement us up to the shear strength, Rs. After the peak value is reached, the bond breaks and tangential force drops to the residual frictional shear strength, Rsr. In Fig. 3(c), the mechanical response on rolling contact direction is similar to that on tangential direction, the moment, M, is a product of rolling stiffness, km, and relative rolling rotation θ, and then it drops to the residual strength, Rrr, that mimics particle shape effect on rolling. Mineral grains

Cement 50

density theory [13] and G-criterion [14] were also proposed for studying crack propagation and rock fracture mode. Distinct element method (DEM) is another suitable method for studying the failure behavior of rocks being able to describe the crack propagation and rock fracture both from the macroscopic and microscopic point of views. In general, any model intending to reproduce the failure process of rocks should allow for the development of a continuous sliding surface [15]. In fact, DEM is able to deal with material failure from micro crack to macro failure [16] and does not require the formulation of complex constitutive models [17,18]. This approach has been adopted to analyze penetration mechanism of cone-penetration tests [19], slope instability problems [20–26], constitutive models for granular material [27,28], the effective stress and shear strength functions of unsaturated granulates [29], strain localization [30,31]. This method has also been applied to model the mechanical behavior of rock considered as a cemented granular material [32,33]. The bonded particle model (BPM) proposed by Potyondy and Cundall [33] has drawn attention for its ability to reproduce the main mechanical behavior of Lac du Bonnet granite using DEM uniaxial and triaxial compression tests. The approach and benefits of BPM are compelling, but the calibration of microscopic parameters to match macroscopic response is inadequate [34]. In order to apply DEM to simulate rock behavior successfully, Wang and Tonon [35,36] developed a distinct element code and applied it to model the mechanical behavior of a crystalline rock such as Lac du Bonnet granite. Their simulations reveal that three stages of stress–strain curve can be well identified and tensile failure appears first, followed by mobilization of residual friction. Christian et al. [37] proposed a progressive failure model and reproduced many features observed during rock failure experiment. However, the bond contact models used by most researchers in their DEM simulations are necessarily hypothesized without any direct experimental verification. Brittle rock can be represented as an assembly of cemented mineral grains that can be described by its internal microstructures. The intact crystalline rock consists of a variety of mineral grains of different sizes, and of close contacts between mineral grains and microscopic defects in the form of cracks and holes [38]. The mechanical behavior of such geo-material is controlled by irregular pre-existing defects, as well by mineral grains and associated contacts. In order to establish an effective bond contact model for these cemented granules in DEM simulation, Delenne et al. [39] first presented an experimental investigation on the mechanical behavior of cemented granules by performing simple tests on a pair of aluminum rods glued together by epoxy resin. This study is attractive and impressive but the effect of normal force on the mechanical behavior of cemented aluminum rods is neglected probably due to technical difficulties. Moreover, the bond material adopted in the experiment is an epoxy resin whose mechanical properties are quite different from a real rock material. Hence, taking these two unsatisfactory aspects into account, Jiang et al. [40,41] conducted a series of experimental tests on a pair of aluminum rods glued together by cement. The cemented rods were prepared by means of a specially designed sample preparation device. Then, the mechanical behavior of the analogous cementation samples was examined in both simple loading and complex loading tests using the newly developed auxiliary loading devices in order to propose a bond contact model for rocks. The main objective of this paper is to perform a distinct element modeling of rock mechanical behavior and fracture propagation. For this purpose, based on the laboratory tests on the analogous cementation sample, a normal force dependent bond contact model was proposed and implemented into a two-dimensional DEM code. A series of laboratory tests on rock, including uniaxial compression, biaxial compression, direct tension as well as Brazilian tension tests

12

p 3

t

rp

Units: mm

Bond thickness: 0.6 Fig. 1. The analogous cementation sample used in experiments (after Jiang et al. [41]).

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

Normal force

Normal force

Normal force

177

Horizontal force Horizontal force

a

c

b Normal force Horizontal force

Normal force

Initial eccentric distance e

Horizontal force Horizontal force

Horizontal force

e

d

Fig. 2. Test schemes designed to characterize mechanical behaviors of the cemented samples: (a) tension test; (b) compression test; (c) shear test under different normal forces; (d) rolling test under different normal forces; and (e) shear-rolling test under different normal forces (after Jiang et al. [40]).

b

kn

Normal bond strength Rt

1 Compressive force

Fs Tangential force

Fn Tensile force

a

Normal displacement

un Residual strength due to friction

Tangential bond strength Rs

ks 1

Residual strength due to friction

Tangential displacement

us

1

ks

Tangential bond strength Rs

c

Rolling resistance

M

Rolling bond strength Rr

km Residual strength due to particle shape and load

1

Rolling rotation

θ

Fig. 3. Mechanical response of bond contact model: (a) normal mechanical response; (b) tangential mechanical response; and (c) rolling mechanical response (after Jiang et al. [40,45]).

2.3. Strength envelopes 2.3.1. Cemented sample with 0.6 mm bond thickness Fig. 4(a) and (b) provides the strength envelope of the cemented sample, with a 0.6 mm bond thickness, expressed on shear strength-normal force and rolling strength-normal force planes, respectively. Note that the black solid lines and red dotted lines shown in Fig. 4 represent peak strength envelopes and residual

strength envelopes, respectively. Fig. 4 shows that the peak shear and rolling strength do clearly depend on the applied normal force. With the increasing of normal force Fn, shear strength Rs and rolling strength Rr increase first, which can be regarded as the frictional part, and then decrease up to compressive strength Rc where the peak strength envelopes intersect with residual strength envelopes which can be regarded as the bond part. So the shear force Rs of the bond can be subdivided into two parts, i.e.

178

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

a

b 24

14

Rolling strength Rr (N.m)

Shear strength Rs (kN)

16

12

Experimental data Eq.(1) Eq.(3)

10 8

Critical normal force

6 4

20

12

Compression strength

Tensile strength 2

4

Compressive strength

0

-Rt 0

Fncr 20 Rc Normal force Fn (kN) 10

30

-10

40

-Rt

0

10 20 Rc Normal force Fn (kN)

30

40

d

c 15

4

Rolling strength Rs

12

Peak shear force Rs

8

Tensile strength

0 -10

Experimental data Eq.(2) Eq.(4)

16

R = μ fR

9 6 3 R = μ fR g (ln(1/f ))

0 0.0

0.3

0.6

3

2

f

R = rβ fR g (ln(1/f)) /6 R

r fR /6

1

fs

0.9

1.2

1.5

0 0.0

stress ratio f

0.3

0.6

0.9

1.2

1.5

stress ratio f

Fig. 4. Strength envelope of the cemented sample with 0.6 mm bond thickness in terms of: (a) shear strength vs. normal force; (b) rolling strength vs. normal force; (c) shear strength envelope; and (d) rolling strength envelope analysis.

frictional part Rsf and bond part Rsb as displayed in Fig. 4(c). In this case it is possible to write Rsf ¼ μb f Rc

ð1Þ


f s 1 Rsb ¼ μb f Rc g s ln f

ð2Þ


1 Rs ¼ Rsf þ Rsb ¼ μb f Rc 1 þ g s ln f

f s !

ð3Þ

where f ¼(Fn þRt)/(Rt þRc) is the stress ratio, μb, gs and fs are the simulated friction coefficient, bond coefficient and envelope shape factor, Rc is the compressive strength, and Rt is the tensile strength. The rolling strength can be treated in the same way Rrf ¼ 16 β b rf Rc

ð4Þ


f r 1 1 Rsb ¼ β b rf Rc g r ln 6 f

ð5Þ


f r ! 1 1 Rs ¼ Rsf þ Rsb ¼ βb rf Rc 1 þ g r ln 6 f

ð6Þ

where βb, gr and fr are the simulated friction coefficient, bond coefficient and envelope shape factor, r is the common radius of two contact particles. As above said, gs and gr determines the percentage of bond part, fs and fr control the envelope shapes, and the values obtained by fitting the experimental data are 0.986, 2.15, 0.761 and 3.055, respectively. After the normal force Fn exceeds compressive strength Rc, the bond is broken and the residual strength increases linearly with

normal force. The residual shear strength Rsr is determined by Mohr–Coulomb criterion; the residual rolling strength Rrr can be described by bond rolling resistance coefficient βb, the common radius r and normal force Fn [45,46]. Then the strength envelopes are re-formulated as follows: Rsr ¼ F n μb

ð7Þ

Rrr ¼ 16 F n βb r

ð8Þ

The predictions given by Eqs. (1)–(8) with these values are also presented in Fig. 4.

2.3.2. Zero-thickness cementation sample In nature, the contact thickness between mineral grains is so thin that the central thickness of two neighboring mineral grains can be regarded as zero. In principle, the bond contact model used for simulating most rock types should be proposed from the laboratory tests on a zero-thickness cementation sample. But the zerothickness cementation sample is much more difficult to prepare in the laboratory due to poor flow ability of cement and such experiments are not yet available. However, the strength envelope of a zero-thickness cementation sample can be obtained through theoretical analysis based on the results of cemented sample with 0.6 mm bond thickness, whose derivation procedure will be introduced in detail elsewhere in the future. If the parameter t in Fig. 1 is set to zero, cementation model becomes a zero-thickness. At this moment, the complete contact between two particles in this model is composed of both an inter-particle contact and an inter-bond contact, and that they transmit the applied normal compression force Fn in a parallel mode. The inter-particle contact force Fnp and

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

inter-bond contact force Fnb can be therefore calculated by

8

In addition, substituting the compressive strength Rc into Eq. (10), the critical normal force Rcbr of a zero-thickness cementation sample is expressed as Rcbr ¼

knb þ knp Rc knb

¼1

40

50

60

70

80

90 100

Rcbr

Eq. (2) Eq. (8) Eq. (10)

60 45 30 15

Critical normal force

Tensile strength 0 -20

-Rt 0

20

40

60

80Rcbr 100

ð15Þ

Fig. 5 presents the complete strength envelope of a zerothickness cementation sample. Three strength relationships are considered: shear strength-normal force, rolling strength-normal force and shear strength-rolling strength. Being shear and rolling strength mainly dependent on the applied normal force Fn, the relationships between shear/rolling strength and normal force consist of three stages (Figs. 5(a) and (b)) When normal force Fn is tensile, only inter-bond contact exists in the zero-thickness cementation sample. The peak strength is therefore completely dependent on bond strength. The peak shear and rolling strength envelopes at this stage are expressed by Eqs. (3) and (6), respectively. When normal force Fn is compressive, but less than critical normal force Rcbr, both the inter-particle friction (rolling resistance) and bond strength contribute to peak shear strength (peak rolling strength), whose values are predicted by Eqs. (11) and (12). If peak strength is reached, the bond breaks, shear force and rolling resistance decrease to their residual value predicted by Eqs. (13) and (14), which are determined by residual inter-particle and inter-bond strength due to friction and rolling resistance. When the normal force Fn is compressive and larger than critical normal force Rcbr, the bond is already broken and peak strength envelopes are coincident with residual strength envelopes, and the corresponding expressions are Eqs. ((13) and (14). Fig. 5(c) presents the shear strength-rolling strength curve with the value on the x-axis representing moment and the value on the y-axis representing shear force. Fig. 5(c) shows that the shear strength-rolling strength envelope is elliptic. Such characteristic shape is also verified by other experimental data [40,41]. The corresponding theoretical equation is expressed by F 2s M 2 þ R2s R2r

30

Normal force Fn (kN)

ed

ð14Þ

20

larg

Rrr ¼ 16 β p rF np þ 16 β b rF nb

10

75

Correspondingly, the residual shear strength Rsr and rolling strength Rrr can be expressed by ð13Þ

Critical normal force

0 -10 -Rt 0

ð12Þ

Rsr ¼ μp F np þ μb F nb

16

Tensile strength

Rolling strength Rr (N.m)

"
 # 1 1 fr 1 Rr ¼ β b rf Rc 1 þ g r ln þ β p rF np 6 f 6

24

ge d

where knp and knb are inter-particle normal contact stiffness and inter-bond normal contact stiffness, respectively. The peak strength of a zero-thickness cementation sample is contributed by bond strength and inter-particle strength when the normal force Fn acting on the contact is compressive. Therefore, the peak shear strength Rs and rolling strength Rr of a zerothickness cementation sample can be expressed by "
 # 1 fs Rs ¼ μb f Rc 1 þ g s ln ð11Þ þ μp F np f

En lar

ð10Þ

Eq.(1) Eq.(7) Eq.(9)

32

En

F nb ¼ F n knb =ðknp þ knb Þ

40

ð9Þ

Shear strength Rs (kN)

F np ¼ F n knp =ðknp þ knb Þ

179

ð16Þ

Shear force Fs Rs

Rr

Moment M

Fig. 5. Strength envelope of cemented sample with zero bond thickness in terms of: (a) shear strength vs. normal force; (b) rolling strength vs. normal force; and (c) shear strength vs. rolling strength.

2.3.3. Comparisons between BPM model and the proposed model A comparison between the BPM model by Potyondy and Cundall [33] and the bond model proposed here has been carried out. The similarities and differences between the two models are summarized in the following. The normal stress and shear stress in the BPM model are derived by the beam theory, simplifying the microscopic bonded particles as a macroscopic beam in which the stresses of the bond are calculated as

σc ¼  τc ¼

Fs A

F n MR þ I A

ð17Þ ð18Þ

where Fn, Fs and M are the normal force, shear force and moment of the contact, A is the contact area, I and R are the moment of

180

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

inertia and radius of the parallel bond. The tensile and shear strengths are a constant value. However, the local effect caused by particles is neglected in the derivation, which violates the Saint-Venant's principle. The rock mechanical model proposed in this study is based on the laboratory tests allowing to avoid the complex analysis of the stress distribution in the bonded particles. Fig. 6 presents the strength envelope of BPM model obtained by a simple transformation of the originally proposed equation. In Fig. 6, σc and τc are the maximum tensile stress and the maximum shear stress, respectively; A, I and R are the area, moment of inertia and radius of the parallel bond, respectively. More details about the model and the parameters can be referred to Potyondy and Cundall [33]. Fig. 6(a) shows that peak shear strength Rs is initially constant, and then increases linearly with normal force Fn; Fig. 6(b) shows that peak rolling strength Rr increases linearly with normal force; Fig. 6(c) shows that the shear strength-rolling strength envelope exhibits a rectangular shape. Accordingly, two bond failure modes are presented: tensile failure and shear failure. In contrast, in the rock model presented in this study, both the shear strength and the rolling strength increase with normal force in a parabolic mode first, and then linearly with normal force (Fig. 5(a) and (b)). The shear strength-rolling strength envelope exhibits an elliptical shape (Fig. 5(c)). Three bond failure modes have been taken into account in this model, namely: tensile failure, shear-rolling failure and compressive failure.

3. Calibration of the bond contact model The bond contact model was implemented into the distinct element code using C þ þ language. While the description of the model is compelling and attractive, it requires calibration of microscopic parameters to match the mechanical behavior of the rock. The experimental data used for calibration in this study comes from [42,43] where a series of laboratory tests were reported on Lac du Bonnet granite, such as uniaxial and biaxial compression tests, direct tension and Brazilian tension tests. As described in the following, the fitting parameters fs, gs, fr and gr are obtained from the experimental results of cemented granules. Fitting the experimental data in Fig. 4(a) and (b), the values of the four fitting parameters are determined. The bond tensile strength Rt is determined by simulating the direct tension test and Brazilian test. The normal and tangential stiffness of particles, kn and ks, and bond compressive strength Rc are determined by simulating a uniaxial compression test. Stiffness parameters kn and ks depend on the elastic modulus and Poisson ratio, and strength parameter Rc depends on the strength of the uniaxial compression test. 3.1. Generation of DEM rock sample To calibrate model properties, a large number of numerical simulations were carried out by adjusting the microscopic parameters until the mechanical properties obtained were consistent with experimental data published in [42,43], according to procedures that will be introduced in the next sections. Microscopic material parameters (Table 1) and particle size distribution were finally obtained. The DEM sample was composed of fifteen classes of particle diameters with a maximum diameter of 2.0 mm, a minimum of 0.5 mm, an average (d50) diameter of 1.3 mm, uniformity coefficient (Cu) of 2.4, and curvature coefficient (Cc) of 1.1. The initial planar void ratio is 0.20, which is representative of a dense sample. The total number of particles in each DEM sample (10,000) ensures macroscopic mechanical responses to remain representative with respect to samples containing more particles. The multilayer undercompaction method, originally proposed by Jiang et al. [48], was employed to generate the DEM rock sample because of its capability at controlling both the homogeneity and density of sample. After the generation, the sample was “consolidated” in 1D under a specific vertical pressure by lowering the top rigid wall with fixed side walls until the equilibrium state. During this process, the microscopic material parameters (Table 1) were assigned to contacts continuously as vertical consolidation pressure increases so that the particles can get into contact with the adjacent particles and bond together. The DEM rock sample aimed at matching Lac du Bonnet granite was finally generated.

Table 1 Microscopic material parameters used in DEM simulations. Parameters

Values 3

Fig. 6. Strength envelope of BPM bond model in terms of: (a) shear strength vs. normal force; (b) rolling strength vs. normal force; and (c) shear strength vs. rolling strength (after Potyondy and Cundall [33]).

Particle density ρs (kg/m ) Initial void ratio e Normal stiffness of particles kn (N/m) Tangential stiffness of particles ks (N/m) Bond tensile strength Rt (N) Bond compressive strength Rc (N) Interparticle friction coefficient μp Interparticle rolling resistance coefficient βp Friction coefficient of bonds μb Rolling resistance coefficient of bonds βb

2700 0.20 1.8  1011 9.47  1010 6.5  104 8.0  107 1.0 1.5 0.5 0.5

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

3.2. Simulations of tension tests The tensile strength of rock is an important mechanical property, and ignoring it may result in engineering failures, especially in underground and rock slope excavations. Direct tension and Brazilian tension tests are often used to determine tensile strength of rock in the laboratory. Thus, these two tests were simulated and the results compared with experimental data in order to calibrate bond tensile strength, Rt.

3.2.1. Direct tension tests In order to obtain the microscopic bond tensile strength Rt, a parametric analysis was carried out by varying values of bond tensile strength Rt while maintaining other parameters constant. A set of direct tension tests was simulated on DEM samples with size similar to those of the specimens used in laboratory. The numerical results reveal that rock tensile strength σt increases linearly with the bond tensile strength, Rt. Model calibration allows us to obtain the stress– strain relationship from direct tension test (Fig. 7(a)). Fig. 7(a) shows the elasto-brittle behavior of the sample and the agreement between experimental [42] and numerical strength σt values. An abrupt stress drop (Fig. 7(a)) is observed at a strain of 0.014%, and this is ascribed to the local bonds breakage during the loading process. The number of tension-induced bond failures increases almost instantaneously when axial stress reaches the peak value. The failure patterns observed in the experiments and in the numerical simulation (Figs. 7(b) and (c)) are in good agreement, and the failure occurs by tension-induced cracks in the intermediate portion of the sample.

a 10

50

Axial stress (MPa)

8

Peak stress from lab

40

experiments (after Martin, [42])

6

30

4

20

2

10

Number of broken bonds

DEM numerical simulation result

Tensile failure 0 0.000

0.005

0.010

b

0.015

0.020

0 0.025

c 45mm

Crack

D=25mm

Fig. 7. Comparison between direct tension test and DEM simulation: (a) stress and number of broken bonds vs. strain; (b) failure pattern observed in laboratory test (after Martin, [42]); and (c) failure pattern observed in DEM numerical simulation.

181

3.2.2. Brazilian tension test Brazilian tension test was used to validate the bond tensile strength Rt obtained by simulation of direct tension test. The diameter of Brazilian sample in the numerical simulation (44.97 mm) is consistent with that of the specimen tested in the lab. The test was run by moving the upper and lower rigid walls towards each other at a constant rate (0.01 m/s). Fig. 8(a) presents the relationships between tensile stress, number of broken bonds and axial strain obtained from the DEM numerical simulation. Fig. 8(a) shows that the stress increases linearly with strain up to a peak value that matches well the experimental data. The number of broken bonds caused by tension and shear-rolling increases linearly from a strain of 0.05% to about 0.16%, when the rate of failure increases abruptly at the peak value of tensile stress which agrees with the experimental results presented in Fig.8(b). Fig. 8(c) illustrates in black the tension-induced bond failure and in red the shear-rolling bond failure. Cracks caused by tension are distributed along the loading direction, while shear-rolling induced cracks are mainly presented at loading points due to stress concentration. The DEM sample was split open by tension-induced cracks, as shown in Fig. 8(d), and in good agreement with the failure pattern observed in the laboratory (Fig. 8(e)). 3.3. Compression tests The stress–strain relationship of brittle rock obtained in compression tests can be classified into five regions, as shown in Fig. 9(a) ([42,43,49,50]): (1) Crack closure region. The initial stage of the stress-strain curve exhibits nonlinear deformability due to closure of pre-existing micro cracks in the specimen, and it depends on the initial crack density and geometry. (2) Elastic deformation region. Once the pre-existing cracks are closed, stress, lateral strain and volumetric strain increase linearly with axial strain. Hence, the elastic deformation of the intact rock specimen dominates this part of stress–strain curves. (3) Crack initiation and stable growth region. This region is characterized by initiation of new cracks and the controlled propagation of the pre-existing cracks. Brace et al. [51] reported that this region begins at a stress level equal to 30–50% of peak stress, and referred to this stress as crack initiation stress σci. In addition, the ending point of this region identifies the transition of volumetric strain from compression to dilatation. (4) Crack damage and unstable growth region. This region was first defined by Bieniawski [52] to identify where the structure of rock specimen changes. The density of cracks increases rapidly and can reach about sevenfold [53]. This prompt increase results in the occurrence of volumetric strain reversal due to lateral strain growing so fast to exceed axial strain. Such change often occurs at an axial stress level of 70–85% of peak stress, and is referred to as crack damage stress σcd. (5) Post peak region. The beginning of this region is marked by the peak stress (σf), which is used to establish the failure strength envelope. At this stage, the macro fracture surface initiates and specimen failure occurs. The three characteristic stress levels, namely, crack initiation stress σci, crack damage stress σcd and peak stress σf, are very important in the stress–strain relationship, but it is very difficult to detect them in laboratory except for the peak stress. 3.3.1. Uniaxial compression tests In order to model the uniaxial compression behavior of Lac du Bonnet granite, other microscopic parameters in the model should be calibrated to match macroscopic response. The identification of

182

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

Fig. 8. Comparison between Brazilian tension test and DEM simulation: (a) stress and number of broken bonds vs. strain; (b) laboratory results (c) distribution of broken bonds classified according to mechanism (black: tensile failure; red: shear and rolling failure); (d) failure pattern observed in DEM numerical simulation; (e) failure pattern observed in laboratory test (after Martin, [42]).

deformability parameters is carried out on DEM samples by modeling uniaxial compression tests. It is found that Young's modulus E increases linearly with normal stiffness kn when fixing the ratio of normal stiffness to shear stiffness, and Poisson ratio υ is only related to the ratio ks/kn. Different combinations of kn and ks/kn were adopted to conduct simulations and the results were compared with the published experiments on Lac du Bonnet granite. The calibrated values for the stiffness are equal to kn ¼ 1.8  1011 N/m and ks ¼9.47  1010 N/m, and they result in simulated macro deformability properties E¼67.53 GPa and υ ¼0.256, which are very close to the average experimental values of E¼69 GPa and υ ¼0.26 [42,43]. The inter-bond friction coefficient μb mainly affects the postpeak mechanical response and it is not clear to what macroscopic parameter it should be calibrated. Thus, μb ¼ 0.5 is used as a reasonable non-zero value. Inter-particle friction coefficient μp ¼1.0 was selected based on the ratio of inter-particle friction coefficient to inter-bond friction coefficient obtained from our experiments on aluminum rods. In addition, the rolling resistance plays a significant role in modeling crack propagation [54], hence the inter-bond and inter-particle rolling resistance coefficients were determined by modeling crack propagation hereinafter and set equal to 0.5 and 1.5, respectively. Another parametric analysis was carried out by varying bond compressive strength Rc while fixing other values. The result reveals that the rock compressive strength σc increases with the bond compressive strength Rc in a logarithmic mode. The bond compressive strength Rc is set to 8.0  107 N by calibration against experimental data. Fig. 9(b) presents the stress-strain relationship and the number of broken bonds with axial strain as obtained from numerical uniaxial compression test. The peak strength σc and elastic modulus E are in good agreement with experimental data shown in Fig. 9(a).

In addition, the macroscopic failure of the rock sample is caused by the initiation and propagation of micro-cracks which can be identified by looking at the history of the number of broken bonds in the DEM simulation (Fig. 9(b)). The failure process of the DEM rock sample can also be classified into five distinct stages as follows: Stage I: No crack occurs and the stress-strain relationship is almost linear which is different from experimental observation. Such discrepancy may arise from the fact that the DEM sample is consolidated before loading which leads to partial closure of initial pores. At the same time the laboratory sample contains a different and particular set of defects. Stage II: No new crack is generated since number of broken bonds remains zero at this stage. Likewise, the elastic deformability plays a dominant role because the relationships between stress, lateral strain, volumetric strain and axial strain are all linear. Stage III: The onset of failure indicates the start of third region and the corresponding stress is referred to as crack initiation stress, σci, in DEM simulations. As a result of crack generation, stress and lateral strain against axial strain depart from linearity. Moreover, the volumetric strain reaches its peak value at the end of this phase due to the increasing rate of lateral strain exceeding that of axial strain. Stage IV: The axial stress level at the volumetric strain reversal marks the beginning of this stage, and it is referred to as crack damage stress, σcd. Besides, it can be observed by the slope of the number of broken bonds-axial strain curve that cracks propagate rapidly, which is in good agreement with the main characteristic of the unstable crack region in real rocks. Stage V: Appearance of the peak stress σf marks the beginning of the post-peak behavior. After the peak, both lateral strain

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

183

f

cd f

ci f

E

VN

V V

Axial stress (MPa) Peak stress

f

200

Crack damage stress

3000 cd

150 ci

Crack initiation stress

2000

100

Tensile failure 50

1000

cc

Number of broken bonds

4000

250

Crack closure stress -0.20

-0.16

-0.12

-0.08

-0.04

Lateral strain (%)

0 0.000.0

Shear and rolling failure 0

0.1

0.2

0.3

0.4

0.5

0.3

0.4

0.5

0.3

Volumetric strain (%)

0.2

Axial compression test

0.1 0.0 -0.1

Compression

Expansion

-0.2 -0.3 -0.4 -0.5 0.0

0.1

0.2

Axial strain (%) Fig. 9. Comparison of stress–strain relationship between experimental and DEM results: (a) stress-strain relationship obtained experimentally from uniaxial compression test on Lac du Bonnet granite (after Martin & Chandler, [43]) and (b) stress–strain diagrams obtained numerically from the DEM uniaxial compression test showing also the number of broken bonds for different failure mechanisms.

and volumetric strain responses display a continuing dilation, while a residual strength determined by friction is reached. Similarly to axial stress, the relationship between broken-bond

number and axial strain tends to stabilize. At this stage, the failure patterns of the DEM sample and laboratory specimen are illustrated in Fig. 10(a) and (d), respectively, which show

184

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

Fig. 10. Failure patterns of rock samples obtained from numerical and experimental uniaxial compression tests: (a) post-peak DEM sample; (b) post-peak contact bond distribution; (c) failure bond distribution classified according to mechanism (black: tensile failure; red: shear and rolling failure); and (d) post-peak rock sample of Lac du Bonnet granite (after Martin [42]). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

900 DEM simulation by the proposed model Laboratory test on Lac du Bonnet grainte (afterMartin, [42]) DEM simulation by the BPM model(after Potyondy and Cundall, [33])

700

83 m

α

2b=8.92 mm

α

3

300

2a=17.83 mm

8 7. =1 2a

m

400

m

106.9 mm

500

106.9 mm

600 7. =1 2a

m

Peak stress σ1 (MPa)

800

200 100 0

0

5

10

15

20

25

30

35

Confining pressure σ3 (MPa)

53.5 mm

53.5 mm

Fig. 11. Strength envelopes obtained from biaxial DEM numerical simulations and laboratory tests.

Fig. 12. Geometry and size of the adopted DEM (a) single-flawed and (b) doubleflawed samples following [44] experimental test preparation.

that both samples develop shear planes rather than the axial splitting plane. In addition, Fig. 10(c) shows that the cracks caused by tension distribute in the sample uniformly, while shear-rolling induced cracks are mainly concentrated in shear planes. However, the experimental stress–strain curve exhibits a brittle behavior, which is missing in the simulated curve. This may result from the different methods adopted to measure strain in the sample. The axial strain was recorded in laboratory by strain gauge attached on the specimen, while in the DEM sample was measured by three measurement circles arranged in the upper, middle and lower part of sample.

BPM and referring to other simulations [34–36], it is observed that the shear strength of BPM model (Fig. 6(a)) is almost independent of normal force until the bond breaks. This makes difficult to model a strongly sloping strength envelope.

3.3.2. Biaxial compression tests The biaxial compression tests on the DEM sample were numerically performed at confining pressures of 1 MPa, 5 MPa, 10 MPa, 20 MPa and 30 MPa in order to calibrate the bond contact model and microscopic parameters. The strength envelopes obtained from DEM tests and experimental tests [42] (Fig. 11) are consistent with each other. To verify the model capabilities with respect to other existing well-known models, some tests have been carried out on the samples using the BPM bond model [33] implemented into PFC [55]. The numerical results by the BPM model are also presented in Fig. 11 and these show that the BPM model does not succeed at modeling the steep slope of the strength envelope for the Lac du Bonnet granite. Comparing the proposed model with

4. Distinct element modeling of cracking in flawed samples It is known that wing cracks initiate at the tip of pre-existing flaws and propagate towards the direction of the major principal stress as load increases [56]. Secondary cracks also start from the tip of pre-existing flaws but grow in two possible directions: (1) coplanar or nearly coplanar to the pre-existing flaw and (2) opposite to the wing crack. The Distinct Element Method (DEM) integrated with the proposed model was employed to simulate crack propagation and rock fracture, whose results were compared with experimental observations on granite samples [44]. Fig. 12 presents the DEM single-flawed and double-flawed samples comparable to those in [44], both in terms of size (53.5  106.9 mm2) and flaw length and bridge distance between the two flaws (17.83 mm and 8.92 mm). In DEM simulation, the pre-existing flaws are generated by deleting the particles that lie in the flaws. In single-flawed samples, as illustrated in Fig. 12(a), the center of flaw is coincident with that of the sample and flaw angle, α, varies from 0o to 90o at 15o increments. The set-up in double-flawed samples is a combination of a horizontal flaw and an underneath inclined flaw, which

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

185

Fig. 13. Cracking patterns in single-flawed prismatic samples of Hwangdeung granite obtained from the uniaxial compression experiments [44] and DEM simulations with different flaw inclination angles: (a) 0o; (b) 30o; (c) 45o; (d) 60o; (e) 75o; and (f) 90o. For each panel (a–f), the experimentally observed features (left), the DEM cracks (center), and the bond type (right) are shown.

is illustrated in Fig. 12(b). Again, the horizontal flaw center is fixed at the midpoint of the sample. The rock bridge represents the link between the center of the horizontal flaw and the top of the inclined flaw. The inclined flaws have an inclination angle, α, of 30–90o at 15o intervals. Because the specimens in laboratory have a 1 mm aperture thickness generated by water-jet system, the preexisting flaws in the DEM simulations were produced by deleting a suitable number of particles. It is worth to say that the same set of properties calibrated on the Lac du Bonnet granite has been used for this set of numerical simulations of a different granite (Hwangdeung granite) [44]. 4.1. Cracking patterns The flawed samples have been loaded under uniaxial compressive condition. Figs. 13 and 14 present the cracking patterns in the single-flawed and double-flawed samples, respectively. The lefthand side of Figs. 13 and 14 show the cracks observed in the

experiments as published in Ref. [44]. The middle and right-hand side images in Figs. 13 and 14 are failure patterns and distributions of bonds for DEM samples, respectively.

4.1.1. Cracking patterns in single-flawed samples The cracking pattern shown in Fig. 13(a) indicates that four cracks initiate at the tip of the pre-existing flaw and propagate along the direction parallel to the major principal stress for sample containing a horizontal pre-existing flaw. For samples containing a pre-existing flaw with inclinations α ¼30–75o, two types of cracks are observed from Fig. 13(b)–(e) in both the experimental specimens and DEM samples: wing cracks and secondary cracks. Both these cracks initiate at the flaw tips, and propagate along the direction of the major principal stress. For α ¼ 90o, Fig. 13(f) shows that only wing cracks initiate at the tip of the flaw and propagate in the direction of the major principal stress. In general, Fig. 13 shows that the crack initiation and propagation are relevant to the

186

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

4

5 2

1

2

1

3

6 4 2

5

3

4

2

6

1

5

6 1

2

4

3

3

6 4

4

2

1 3

1

6

5

6 2

4

1

4 2

5 1

3

5 6 1

2

5

6

3

3

5 2

1

5 4 2 1

6

2

3

5 2

5

6 4

1

3

6 1

5

4

5

4

5 1

6

4

5

1

3

Fig. 14. Cracking patterns in double-flawed prismatic samples of Hwangdeung granite obtained from the uniaxial compression experiments [44] and DEM simulations with different flaw inclination angles: (a) 30o; (b) 45o; (c) 60o; (d) 75o; and (e) 90o. For each panel (a–e), the experimentally observed features (left), the DEM cracks (center), and the bond type (right) are shown. Numbers in figures represent the order of fracture evolution.

inclinations with respect to loading direction, and that cracking patterns obtained from numerical simulations are in good agreement with experimental results, both in terms of fracture length and direction.

4.1.2. Cracking patterns in double-flawed samples Fig. 14 provides cracking patterns in double-flawed samples obtained from experimental and numerical uniaxial compression tests. Two wing cracks are relevant to the coalescence of two flaws in both the experimental specimen and DEM sample. One wing crack develops from the middle of the inclined flaw to the right tip of the horizontal flaw (Fig. 14(a), (1), whereas another runs from the upper tip of the inclined flaw to the middle part of the horizontal flaw (Fig. 14(a), (2). In addition, a secondary crack develops at the lower tip of inclined flaw (Fig. 14(a), (3), whereas other cracks initiate from the middle and tip of the horizontal flaw, and propagate in the upward and downward directions (Fig. 14(a), (4) and (5). The experimental and the numerical observations for samples with α ¼45o (Fig. 14(b)) shows that the coalescence is associated to

the developing of two cracks: one is the link of the upper tip of the inclined flaw and the right tip of horizontal flaw (Fig. 14(b), (1); the other crack propagates from the upper tip of the inclined flaw to the middle of horizontal flaw (Fig. 14(b), (2). In addition, a wing crack emanating from the lower tip of the inclined flaw propagates downward (Fig. 14(b), (3); two cracks initiating from the tips of horizontal flaw propagate downward and upward (Fig. 14(b), (4) and (6). Finally, another crack initiating from the middle part of the horizontal flaw propagates in the direction parallel to the maximum compressive stress (Fig. 14(b), (5). In the case of α ¼60o, one crack links the lower tip of the inclined flaw to the right tip of horizontal flaw both in the experimental and numerical samples (Fig. 14(c), (1), whereas another runs from the upper tip of the inclined flaw to the middle part of the horizontal flaw (Fig. 14(c), (2). Moreover, both the experimental observation and numerical simulation show that a wing crack develops at the lower tip of the inclined flaw (Fig. 14(c), (3), one crack initiates at the left of horizontal flaw (Fig. 14(c), (6) and one at the middle (Fig. 14(c), (5). They propagate in a stable manner and reach the upper boundary of the samples.

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

In the laboratory specimen and DEM sample for α ¼75o, there is one coalesced crack from the lower tip of inclined flaw to the right tip of horizontal flaw (Fig. 14(d), (1). Additionally, Fig. 14(d) and (e) shows that two wing cracks propagate from the left and right tips of horizontal flaw to the upper boundary (Fig. 14(d) and (e), (5) and (6), and one wing crack initiating at the left tip of horizontal flaw propagates to the lower boundary (Fig. 14(d) and (e), (4). In conclusion, the cracking patterns observed in experiments and numerical tests under uniaxial compression loading are in good agreement. For the single-flawed samples, all cracks initiate from the flaw tip. In the double-flawed samples, the cracks may develop at the tip of inclined flaw or likely initiate at the middle of horizontal flaw, or possibly first emerge at the tip of horizontal flaw depending on the flaw inclination and geometry. Although many features of the cracking pattern were captured by the proposed DEM model, the cracking path is somewhat different from experimental observation, especially in single-flawed samples. Such discrepancy may arise from the fact that the use of an assembly of bonded circular particles is difficult to fully resemble the heterogeneous rock specimen, especially for the distribution of the initial irregular micro cracks in the rock specimen. 4.2. Stress analysis Fig. 15 presents the relationships between peak stresses normalized with respect to the maximum experimental value in the

a

1.8 DEM numerical simulation by the proposed model Laboratory test result (after Lee & Jeon, [44]) DEM numerical simulation by BPM model (after Lee & Jeon, [44])

Normalized peak stress

1.6 1.4 1.2 1.0

0.6 0

10

20

30

40

50

60

70

80

90

ο

Flaw inclination angle α λ λ

Normalized crack initiation stress

b

single-flawed samples (Fig. 15(a)), crack initiation stresses normalized with respect to the maximum experimental value in the double-flawed samples (Fig. 15(b)) and flaw inclination angle as obtained from DEM numerical simulations and experimental tests. Fig. 15(a) shows that the normalized peak stresses decrease with increasing α from 0o to 30o, and then increase for values larger than 30o. The results of numerical simulations are slightly larger than those of experiments but much closer than DEM simulations by the BPM model. Moreover, in single-flawed samples, the preexisting flaw for α ¼30o reduces the strength of intact DEM sample most, while flaw for α greater than 75o lowers the strength least. This argument is also supported by another experimental model conducted by Lajtai [57]. Fig. 15(b) shows that the normalized crack initiation stresses increase when flaw angle α is smaller than 45o, then decrease slightly with increasing of α from 45o to 60o, and finally increase again. The trend obtained from our DEM simulation approaches experimental data a bit more than DEM simulation by the BPM model, especially when inclination angles α are in the range of 30–60o. In addition, in the cases of double pre-flawed samples, when the location of horizontal flaw is fixed, the flaw with inclination of 30o initiates most easily, and the opposite occurs for flaw inclination larger than 75o. However, Fig. 15(b) shows that the crack initiation stresses obtained from our numerical simulations demonstrate differences from laboratory tests quantitatively. This may be due to the fact that the identification method of the crack initiation stress in the numerical simulation is different from that in the experiments. In addition, the use of mechanical properties calibrated for the Lac du Bonnet granite can explain the shifting of the numerical results with respect to the experimental ones, although the almost constant shifting suggests the significance of the calibrated values to simulate the behavior of a crystalline granite. Nevertheless, the proposed model succeeds at predicting the role played by the geometry of the flaw assemblage, as demonstrated by the similar trend in the relationships of Fig. 15(a) and (b).

5. Discussion

0.8

0.4

187

1.6 1.4 1.2

DEM numerical simulation by the proposed model Laboratory test result (after Lee & Jeon, [44]) DEM numerical simulation by BPM model (after Lee & Jeon, [44])

1.0 0.8 0.6 0.4 0.2 20

30

40

50

60

70

80

90

ο Flaw inclination angle α τ τ Fig. 15. Peak stress and crack initiation stress versus flaw inclination angle as obtained from DEM uniaxial compressive tests and laboratory tests [44] on Hwangdeung granite samples with different pre-existing flaws: (a) single preexisting flaw and (b) double pre-existing flaws.

One of the main aims of this paper is to introduce a new contact model based on the cemented contact experiments whose results show that the breakage of bonds can be induced by normal, shear and torsion forces independently or in combination. Rock can be represented as a heterogeneous material comprised of cemented grains [3]. This can be a suitable representation both for crystalline as well as for some sedimentary rocks. In this paper, the grains are regarded as the cemented granules whose mechanical behaviors were investigated by carrying out series of loading tests. Based on the experimental results, a model was established in which the grains were circular and the cement existed at and in vicinity of the contacts between the particles. The mechanical response and strength envelope are both based on the circular particles. Such a choice causes the DEM specimens to show slightly anisotropic mechanical behavior, since the multi-layer undercompaction method was used to generate DEM samples. In contrast, non-circular particle shapes may induce a strongly anisotropic mechanical behavior if the multi-layer under-compaction method is used. Nevertheless, the effect of particle shapes will not change the conclusions in this paper considering the more isotropic characteristics of the studied lithologies. Among the contact stress theories for uncemented elastic spheres at contact, the Hertzian contact stress theory employs nonlinear – spring stiffness in normal and tangential direction. In contrast, in this study constant stiffnesses are used in normal, tangential and rolling direction. This can significantly reduce computation time of DEM analyses. Such a slight difference on

188

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

the elasticity of the particles may lead to a slight difference on the elasticity of any assembly of granulates, but will not affect its macroscopic non-linear mechanical behavior resulting from the re-arrangement of the particles. In fact, it is known that the deformation of the rock is mainly induced by the generation and propagation of micro cracks. Hence, the employment of constant stiffnesses in the model may have a slight influence on the elastic behavior but it will not affect the main mechanical behavior of the DEM material as presented in this paper. For the laboratory samples, the initial non-linear stress-strain curves observed in the compression test result from the closure of the pre-existing micro-cracks. However, in this paper all the contacts between the particles are bonded with cementation and no pre-existing micro cracks exist in the DEM sample. Once the DEM sample is loaded and no contact reaches the plastic stage, the DEM sample demonstrates the initial linear rather than non-linear stress–strain behavior. Such initial non-linearity may be numerically simulated by employing nonlinear spring stiffness in normal, tangential, and rolling direction in the contact model. However, since the preexisting micro cracks do not play important role in the mechanical behavior of intact rock before the crack initiation stage, the preexisting micro cracks are not taken into account in this study. The particle size distribution is another controlling variable affecting the results of DEM simulations in combination with initial void ratio of the specimen, which jointly form particle packing patterns. The effect can be finally controlled by a microscopic term relative to the average coordination number. The effect may be quantitatively different but will be qualitatively identical for different models (i.e. the BPM and this model). In this paper, the particle size distribution was chosen to produce a densely packed specimen in which the largest diameter 2 mm is close to the reported value for crystals [3] of Lac du Bonnet granite. This numerical model is validated by capturing the mechanical behavior of Lac du Bonnet granite, and the effect of size distribution will not change the conclusions in this paper. The DEM simulation results show that this model can simulate the mechanical behaviors of a specific rock (i.e. Lac du Bonnet granite), on the aspects of (a) the brittle failure, (b) ratio of uniaxial compression and tensile strength, (c) strength envelope, and (d) propagation of flaws. It is also shown that neither the ratio of uniaxial compression to tensile strength nor the strength envelope can be qualitatively simulated by the BPM. Potyondy and Cundall [33] also reported that the slope of the strength envelope simulated by BPM is too low for the use of circular and spherical grains in the model. By establishing the clumped particle model [4], the strength envelope fitted the experimental results well. The reason why the BPM model shows this kind of behavior comes mainly from the fact shown in Eqs. (9) and (10) that tangential force does not contribute to the failure controlled by σc, and neither normal force Fn nor moment M contributes to the failure controlled by τc. Since macroscopic mechanical behavior is strongly linked to microscopic behavior at contacts, this demonstrates the importance of microscopic constitutive law for bonding of particles in controlling the mechanical behavior of rocks Eventually, it must be stressed that microscopic grain structures are generally 3D and that some tests are characterized by a more axisymmetric geometry (uniaxial tests on Lac du Bonnet granite), whereas other can be more easily simulated in 2D (Brazilian tests on Lac du Bonnet granite and uniaxial compression on pre-flawed prismatic samples of Hwangdeung granite). Although 3D DEM appears to be a promising option for the analyses, the 2D modeling is the only possible choice at this stage. This is because (a) a rock specimen with a large number of grains is needed to investigate the propagation of flaws. This is impossible in 3D analyses due to the present computational capacity of PCs; (b) a 3-D bond contact model is needed in 3D DEM, for which

apparatus and devices will have to be developed at first to carry out 3D tests on the bonded granules idealized by two glued aluminum spheres; (c) in terms of visualizing internal information on the propagation of flaws, 2D DEM is easier than 3D DEM; (d) one important application of DEM is to simulate large-scale boundary-value problems in geotechnical engineering using current computers, in which the size and boundary effects must be reduced to a minimum, which requires an extremely large number of particles. Again, this is only possible for 2D DEM models due to current PCs computing capacity.

6. Conclusions A series of laboratory experiments were conducted on a pair of cemented aluminum rods aimed at stimulating the mechanical behavior and fracturing of rocks via the distinct element method (DEM). Based on the experimental data, a normal force dependent bond contact model was developed and compared with BPM model. By using a DEM code into which the proposed model was embedded, uniaxial and biaxial compression tests, direct tension and Brazilian tension tests were numerically simulated. The results are compared with the published experimental data of Lac du Bonnet granite [42,43] for calibrating the model and identifying microscopic material parameters for such a type of crystalline rocks. With the identified parameters, an investigation of crack propagation and rock fracture was then performed on both single-flawed and double-flawed samples of different flaw angles. The peak stresses in single-flawed samples and crack initiation stresses in double-flawed samples were analyzed and compared with those obtained from laboratory tests performed on a different granite but with similar characteristics and from DEM simulations using BPM model [33]. The following conclusions can be drawn from this study: The proposed bond contact model was successfully calibrated to match experimental data [42,43] on intact Lac du Bonnet granite. The limitations of a gently sloping strength envelope and high ratio between tensile strength and compressive strength frequently encountered in DEM simulation of rocks were tackled and solved. In addition, the recorded evolution of the number and type of bond failures indicate that tension-induced cracks initiate first and the proportion of these cracks is higher. An excellent agreement is observed between our DEM numerical and other reported experimental [44] data on the crack propagation patterns in both DEM single-flawed and double-flawed samples which depend strongly on the flaw orientation and geometry. The DEM numerical results were generally in a good agreement with the experimental observation [44] on the peak stress in single-flawed samples and crack initiation stress in double flawed samples, showing that the single flaw of α ¼ 30o lowers the strength most and the sample containing double flaws with α ¼30o is most easily to initiate, even if values for the microscopic properties were calibrated on a different granite (Lac du Bonnet). This indicates that both the model and the parameter values are highly suitable to validate the behavior of such a type of crystalline rocks. Our approach has several advantages over previous models: the bond contact model is based on experimental observations; the crack initiation stress is identified explicitly by onset of broken bonds; the mechanical behavior of granite was reproduced successfully. Nevertheless, further investigation on the cracking mechanism from the microscopic point of view is necessary. Our future efforts will include: the use of this model to simulate boundary-value problems, such as the instability of rock slopes; performing 3D tests on bonded spheres in the laboratory to obtain the 3D contact model for DEM and to investigate macroscopic and microscopic mechanical behavior of rocks in simple tests.

M. Jiang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 175–189

Acknowledgement The research was funded by the Major Project of Chinese National Programs for Fundamental Research and Development (973 Program) with Grant nos. 2011CB013504 and 2014CB046901, China National Funds for Distinguished Youth Scientists with Grant no. 51025932, and the FP7-PEOPLE-2011-IRSES with Grant agreement no. PIRSES-GA-2011-294976. All of these supports are greatly appreciated. Reference [1] Shao JF, Chau KT, Feng XT. Modeling of anisotropic damage and creep deformation in brittle rocks. Int J Rock Mech Min Sci 2006;43(4):582–92. [2] Chen WZ, Li SC, Zhu WS, Qiu XB. Experimental and numerical research on crack propagation in rock under compression. Chin J Rock Mech Eng 2003;22 (1):18–23 [in Chinese]. [3] Li YP, Chen LZ, Wang YH. Experimental research on pre-cracked marble under compression. Int J Solids Struct 2005;42(9–10):2505–16. [4] Reyes O, Einstein HH. Failure mechanisms of fractured rock-a fracture coalescence model. In: Proceedings of 7th international congress on rock mechanics, Germany, Aachen, vol. 1; 1991. p. 333–40. [5] Bobet A, Einstein HH. Numerical modeling of fracture coalescence in a model rock material. Int J Fract 1998;92(3):221–52. [6] Scavia C. The displacement discontinuity method for the analysis of rock structures: a fracture mechanics. In: Aliabadi MH, editor. Fracture of rock, WTT pressure, computation. Boston: Mechanics Publications; 1999. p. 39–82. [7] Vasarhelyi B, Bobet A. Modeling of crack initiation, propagation and coalescence in uniaxial compression. Rock Mech Rock Eng 2000;33(2):119–39. [8] Tang CA, Yang WT, Fu YF, Xu XH. A new approach to numerical method of modeling geological process and rock engineering problems-continuum to discontinuum and linearity to nonlinearity. Eng Geol 1998;49(3–4):207–14. [9] Tang CA, Wong RHC, Chau KT, Lin P. Analysis of crack coalescence in rock-like materials containing three flaws-Part II: numerical approach. Int J Rock Mech Min Sci 2001;38(7):925–39. [10] Wong RHC, Tang CA, Chau KT, Lin P. Splitting failure in brittle rocks containing pre-existing flaws under uniaxial compression. Eng Fract Mech 2002;69 (17):1853–71. [11] Erdogan F, Sih GC. On the crack extension path in plates under plane loading and transverse shear. ASME J Basic Eng 1963;85:519–27. [12] Hussian MA, Pu EL, Underwood JH. Strain energy release rate for a crack under combined mode I and mode II. Fracture Analysis. ASTM STP 560. American Society for Testing and Materials 1974;560:2–28. [13] Sih GC. Strain–energy–density factor applied to mixed mode crack problems. Int J Fract 1974;10:305–21. [14] Shen B, Stephansson O. Modification of the G-criterion for crack propagation subjected to compression. Eng Fract Mech 1994;47:177–89. [15] Li LC, Tang CA, Zhu WC, Liang ZZ. Numerical analysis of slope stability based on the gravity increase method. Comput Geotech 2009;36(7):1246–58. [16] Potyondy DO. Simulation stress corrosion with a bonded-particle model for rock. Int J Rock Mech Min Sci 2007;44(5):677–91. [17] Cundall PA. A computer rock model for simulating progressive large scale movements in blocky rock systems. In: Proceedings of the symposium of international society of rock mechanics, vol. 1; 1971. Nancy, France. Paper No. П-8. [18] Cundall PA, Strack ODL. A discrete element model for granular assemblies. Geotechnique 1979;29(1):47–65. [19] Jiang MJ, Yu HS, Harris D. Discrete element modeling of deep penetration in granular soils. Int J Numer Anal Methods Geomech 2006;30(4):335–61. [20] Utili S, Crosta GB. Modeling the evolution of nature subject to weathering: 2. Discrete element approach. J Geophys Res 2011;116(F1):1–17. [21] Scholtes L, Donze FV. Modeling progressive failure in fractured rock mass using 3D discrete element method. Int J Rock Mech Min Sci 2011;52:18–30. [22] Utili S, Nova R. DEM analysis of bonded granular geomaterials. Int J Numer Anal Methods Geomech 2008;32:1997–2031. [23] Thompson N, Bennett MR, Petford N. Analyses on granular mass movement mechanics and deformation, with distinct element numerical modeling: implications for large-scale rock and debris avalanches. Acta Geotech 2009;4:233–47. [24] Jiang MJ, Murakami A. Distinct element method analyses of idealized bondedgranulate cut slope. Granul Matter 2012;14:393–410. [25] Crosta GB, Calvetti F, Imposimato , Roddeman D, Frattini P, Agliardi F. Granular flows and numerical modelling of landslides. Report of damocles project; August 2001.

189

[26] Calvetti F, Crosta G, Tatarella M. Numerical simulation of dry granular flows: from the reproduction of small-scale experiments to the prediction of rock avalanches. Riv Ital Geotec 2000;2:21–38. [27] Jiang MJ, Harris D, Yu HS. Kinematic models for non-coaxial granular materials: part II: evaluation. Int J Numer Anal Methods Geomech 2005;29 (7):663–89. [28] Jiang MJ, Harris D, Zhu HH. Future continuum models for granular materials in penetration analyses. Granul Mater 2007;9:97–108. [29] Jiang MJ, Leroueil S, Konrad JM. Insight into shear strength functions of unsaturated granulates by DEM analyses. Comput Methods Appl Mech Eng 2004;31(6):473–89. [30] Jiang MJ, Zhang WC, Sun YG, Utili S. An experimental investigation on loose cemented granular materials via DEM analyses. Granul Mater 2013;15 (1):65–88. [31] Jiang MJ, Yan HB, Zhu HH, Utili S. Modeling shear behavior and strain localization in cemented sands by two-dimensional distinct element method analyses. Comput Geotech 2011;38:14–29. [32] Ghazvinian A, Sarfarazi V, Schubert W. A study of the failure mechanism of planar non-persistent open joints using PFC2D. Rock Mech Rock Eng 2012;45 (5):677–93. [33] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [34] Cho N, Martin CD, Sego DC. A clumped particle model for rock. Int J Rock Mech Min Sci 2007;44(7):997–1010. [35] Wang Y, Tonon F. Modeling Lac du Bonnet granite using discrete element model. Int J Rock Mech Min Sci 2009;46(7):1124–35. [36] Wang Y, Tonon F. Calibration of a discrete element for intact rock up to its peak strength. Int J Numer Anal Methods Geomech 2010;34:447–69. [37] Christian E, Robert S, Peter E. A discrete element model to describe failure of strong rock in uniaxial compression. Granul Mater 2011;13(4):341–64. [38] Lan HX, Martin CD, Hu B. Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. J Geophys Res 2010;115(B01202):1–14. [39] Delenne JY, Youssoufi MSE, Cherblanc F, Benet JC. Mechanical behavior and failure of cohesive granular materials. Int J Numer Anal Methods Geomech 2004;28(15):1577–94. [40] Jiang MJ, Sun YG, Li LQ, Zhu HH. Contact behavior of idealized granules bonded in two different interparticle distances: an experimental investigation. Mech Mater 2012;55:1–15. [41] Jiang MJ, Sun YG, Xiao Y. An experimental investigation on the mechanical behavior between cemented granule. Geotech Test J 2012;35(5):678–90. [42] Martin CD. The strength of massive Lac du Bonnet granite around underground openings. Winnipeg: University of Manitoba; 1993 [Ph.D. thesis]. [43] Martin CD, Chandler NA. The progressive fracture of Lac du Bonnet granite. Int J Rock Mech Min Sci 1994;31(6):643–59. [44] Lee H, Jeon S. An experimental and numerical study of fracture coalescence in pre-cracked specimen under uniaxial compression. Int J Solid Struct 2011;48 (6):979–99. [45] Jiang MJ, Yu HS, Harris D. Bonds rolling resistance and its effect on yielding of bonded granulates by DEM analyses. Int J Numer Anal Methods Geomech 2006;30(8):723–61. [46] Jiang MJ, Yu HS, Harris D. A novel discrete model for granular material incorporating rolling resistance. Comput Geotech 2005;32(5):340–57. [47] Jiang MJ, Yu HS, Leroueil S. A simple and efficient approach to capturing bonding effect in naturally micro-structured sands by discrete element method. Int J Numer Methods 2007;69(6):1158–93. [48] Jiang MJ, Konrad JM, Leroueil S. An efficient technique for generating homogeneous specimens for DEM studies. Comput Geotech 2003;30 (7):579–97. [49] Yang SQ, Dai YH, Han LJ, Jin ZQ. Experimental study on mechanical behavior of brittle marble samples containing different flaws under uniaxial compression. Eng Fract Mech 2009;76(12):1833–45. [50] Tapponnier P, Brace WF. Development of stress-induced microcracks in westerly granite. Int J Rock Mech Min Sci 1976;13:395–430. [51] Brace WF, Paulding BW, Scholz C. Dilatancy in the fracture of crystalline rocks. J Geophys Res 1966;71(16):3939–53. [52] Bieniawski ZT. Mechanism of brittle fracture of rock. Parts І–Ш. I Int J Rock Mech Min Sci 1967;4:395–430. [53] Hallbauer DK, Wagner H, Cook NGW. Some observations concerning the microscopic and mechanical behavior of quartzite specimen in stiff, triaxial compression tests. Int J Rock Mech Min Sci 1972;10:713–26. [54] Wang YC, Mora P. Modeling wing crack extension: implication for the ingredients of discrete element model. Earth Env Sci 2008;165:609–20. [55] Itasca Consulting Group Inc. PFC2D/3D (Particle Flow Code in 2/3 Dimensions), Version 3.1. Minnesota; 2004. [56] Bobet A. The initiation of secondary cracks in compression. Eng Fract Mech 2000;66(2):187–219. [57] Lajtai EZ. Brittle fracture in compression. Int J Fract 1974;10(4):525–36.