- Email: [email protected]

Contents lists available at ScienceDirect

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

Technical Note

Discrete element modeling of anisotropic rock under Brazilian test conditions K. Duan, C.Y. Kwok n Department of Civil Engineering, The University of Hong Kong, Hong Kong

art ic l e i nf o Article history: Received 7 August 2014 Received in revised form 16 April 2015 Accepted 20 April 2015 Keywords: Brazilian test Inherently anisotropic rock DEM Fracture pattern Micro-scale parameters

1. Introduction Anisotropy is one of the most distinct features that must be taken into account in rock engineering as many rocks exposed near the Earth's surface show well-deﬁned fabric elements in terms of bedding, stratiﬁcation, layering, foliation, ﬁssuring or jointing. Such rocks are said to be inherently anisotropic as their properties (physical, mechanical and hydraulic) vary with direction. Rock anisotropy affects different geotechnical operations like rock cutting performance, hydraulic fracture propagation in shale reservoirs, and development of excavation-induced damage around underground structures.1 Therefore, an accurate understanding of anisotropic properties is required for the design and construction of related rock engineering projects. Among the many mechanical parameters, tensile strength is a key one because rocks are in nature much weaker in tension than in compression. Many rock mechanics applications like stability of underground excavations, rate of rock blasting and propagation of hydraulic fractures are highly dependent on the tensile strength. The Brazilian tensile test (diametrical compression of circular discs)2 has been widely adopted in laboratory to determine the tensile strength of rock materials.3–7 The effect of anisotropy on the mechanical responses of rock discs under Brazilian test conditions has been investigated by various approaches. n

Corresponding author. E-mail address: [email protected] (C.Y. Kwok).

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

The theoretical relation between stresses and strains within a disc of anisotropic material under diametrical loading was proposed by Lekhinitskii.8 Based on this solution, Amadei et al.9 presented a comprehensive analytical solution to measure the tensile strength of anisotropic rocks. Chen et al.10 presented an analytical procedure to determine the deformability and tensile strength of anisotropic rocks from the results of Brazilian tests. However, this solution is implicit and numerical techniques are needed to calculate the stress ﬁeld. In the laboratory, the effect of anisotropy on the indirect tensile strength and failure patterns under diametrical loading conditions has been investigated by conducting Brazilian tests on various types of rocks.11–15 In general, the failure of anisotropic rocks under diametrical compression is very complex in terms of fracture mode and direction. Vervoort et al.16 generalized the experimental results of nine different anisotropic rocks into four trends of the Brazilian tensile strength (BTS) which are normalized by the strength for loading perpendicularly to weak planes versus anisotropy angles as illustrated in Fig. 1. ● Trend 1: constant value over the entire anisotropy angles (Fig. 1 (a)); ● Trend 2: constant value between 0° and 45°, followed by a linear decrease (Fig. 1(b)); ● Trend 3: systematic decrease of Brazilian tensile strength over the entire interval (Fig. 1(c)); ● Trend 4: decrease of strength from very low anisotropy angles (between 0° and 30–40°) and followed by a leveling off (Fig. 1(d)).

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

47

Fig. 1. Variation of average failure load, relative to the failure load for loading perpendicular to weak planes (a) trend 1; (b) trend 2; (c) trend 3; and (d) trend 4.16

The fracture types along the weak layers and inside rock matrix were also investigated for these four different trends in Ref. 16. However, possible important parameters (including the strength of weak layers, ratio between tensile strength and cohesion, and anisotropic Young's modulus) and the ways in which they can affect variations of strength and fractures were not investigated in detail. Numerical simulations have been performed to understand the crack initiation and propagation process of anisotropic rocks. The anisotropic mechanical behaviors of Opalinus Clay in Brazilian test and compression test have been studied through combined ﬁnitediscrete element method (FEM/DEM).17,18 In this model, a transversely isotropic elastic constitutive law was implemented to describe the elastic response, while DEM algorithms and non-linear fracture mechanics principles were employed to capture rock fracturing. The fracture behavior of slate rock has been modeled in Ref. 19 based on two-dimensional discrete element (UDEC) where the schistosity is conceptualized as a set of parallel continuous weak layers. Most recently, the particle based DEM was adopted to simulate the behaviors of transversely isotropic rock by employing also continuous smooth joints to represent weak planes.20 However, with careful examination of microstructure of intact anisotropic rocks such as sedimentary or metamorphic rocks, the inherently anisotropy is not necessarily be straight and continuous at microscopic scale.13 Therefore, there is a challenge to develop a more realistic numerical approach to explicitly represent the existence of weak layers in micro-scale. In this study, a new numerical approach is proposed to represent the weak layers in inherently anisotropic rocks more realistically.21 The effect of strength, stiffness and number of weak layers with different loading direction will be examined in detail. The numerical results are compared with previous experimental results quantitatively and failure mode and fracture direction are

investigated. The numerical simulations link the strength anisotropy, deformation behaviors and fracture patterns on sample scale to the micro-properties of weak layers.

2. Numerical methodologies DEM has been successfully used in modeling the behaviors of rocks under different stress conditions.22–27 Moreover, DEM models provide an avenue to investigate the effect of microscopic properties on the macro-scale response, like role of grain interlocking on strength,28 effect of porosity on the deformability and strength,29,30 effect of pore size and pore distribution.31 The advantage of DEM over other continuum methods is its ability to explicitly model the initiation and propagation of cracks from micro-scale to macro-scale without applying complex constitutive laws.32 In this study, DEM models will be generated to represent inherently anisotropic rocks based on the bonded particle model (BPM) and smooth-joint model provided by PFC2D.33 Here, a brief introduction of these two models is provided. 2.1. The bonded particle model (BPM) In PFC models, rock is represented as an assembly of discs (2D) or particles (3D) bonded at their contacts.22 These discrete elements can move with respect to each other and the bonds between them can break once the stress induced by applied load exceeds the corresponding bond strength. Once a bond fails, the stress will be redistributed which may further induce new bond breakage. In this way, meso- or macro-fracture can form as coalescences of micro-cracks. Detailed introduction of the algorithm and related force-displacement relationship for BPM can be found in Ref. 33. The

48

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

k s = (k s/A) + k

s

(3)

μ = μc

(4)

σc = σc

(5)

cb = τc

(6)

φb = 0

(7)

where k¯ n and k¯ s are normal and shear stiffness per unit area (stress/displacement), respectively, λ¯ is radius multiplier, μ is friction coefﬁcient, σc and cb are normal strength and cohesion of the smooth joint and φb is the friction angle. These properties can be overwritten by assigning properties of smooth joint directly. As opposed to the parallel bond model, the shear strength of a smooth joint model is not assigned directly, but is instead determined as follows:

τ = cb + σ tan φb

Fig. 2. Construction of inherently anisotropic rock model.

macro-properties of DEM models for rocks are determined by the related micro-parameters including the size distribution of particles, strength of parallel bond (normal strength σ¯c and shear strength τ¯c ), stiffness of particles (k n and k s ), stiffness of parallel s n bonds (k¯ and k¯ ) and friction coefﬁcient between particles (μ ). c

One problem is that these parameters cannot be measured directly in laboratory, thus an intensive calibration process is necessary. The general procedure based on trial-and-error to choose the micro-parameters has been presented by Itasca.33 Studies have been carried out to explore the relationship between microparameters and macro-properties which can serve as guidance for the selection of parameters for a speciﬁc rock type.30, 34–36 2.2. The smooth-joint model The smooth-joint model was ﬁrst proposed to model the presence of joints in fractured rock masses.37,38 The smooth-joint model simulates the behavior of an interface regardless of the local particle contact orientation along the interface. The behavior of a joint can be modeled by assigning smooth joint models to all contacts between particles that lie on opposite side of the joint. A typical smooth-joint contact in 2D is illustrated in Fig. 2. When the smooth-joint model is created, the parallel bond will be deleted and replaced by the smooth joint contact model. Particles intersected by a smooth-joint can pass through each other by sliding along the pre-deﬁned joint surface rather than move around one another. The properties of the smooth joint will be inherited from the parallel bond in the following ways:

λ¯ = λ¯ pb

k n = (k n/A) + k

(1) n

(2)

(8)

where τ is the shear strength of the smooth joint when the normal stress acting on it is σ . As demonstrated by Eq. (8), the shear strength of smooth joint model is determined by the combination of cohesion (cb ), friction angle (φb ) and normal stress (σ ) acting on it. The smooth-joint model has been successfully applied to the study of fractured rock mass failure,37,39 scale effect,40 representative element volume (REV),41 and anisotropic behavior of jointed rock masses.42,43 Most recently, Park and Min20 performed DEM simulations of transversely isotropic rocks by employing smooth joint model to represent the weak planes. However, all of these applications are based on a series of continuous smoothjoint contacts which represent fully-persistent discontinuities often encountered in fractured rock mass. This study proposes a new numerical approach to explicitly represent the inherent anisotropy by imposing individual smooth-joint contacts on the bonded particle model. The algorithm and procedure for the generation of inherently anisotropic samples will be discussed in the next section.

3. Genesis of anisotropic samples In this study, the simplest form of anisotropy, planar anisotropy, which has a set of parallel planes of weakness, is modeled. For the 2-D model, all the specimen cross sections are assumed to be orientated perpendicular to the strike of the bedding planes. The generation of sample with horizontal weak layers is presented in Fig. 2. The ﬁrst step is to create an isotropic model. To introduce the horizontal anisotropy, any sub-horizontal parallel bonds (those dipping within 20° and þ20°) will be removed and replaced with horizontal smooth-joint contacts (dipping 0°). Similarly, samples with different anisotropy angles can be generated. In this study, seven anisotropy angles (θ ¼0°, 15°, 30°, 45°, 60°, 75°, and 90°) are considered. The sample used for the indirect tensile test is a disc with diameter of 50 mm and unit thickness for the 2D model. The particle size follows a uniform distribution with Rmin ¼0.2 mm and Rmax/Rmin ¼1.66. There are altogether 8182 discs. In the BPM, the response of sample is sensitive to particle size and size distribution as has been demonstrated by previous studies.22,34 The Brazilian strength was found to be affected by the ratio of particle size to Brazilian disc diameter.22 Since the major objective of this study is to investigate the effect of smooth joint properties on the Brazilian tensile test behaviors, we do not aim at investigating the sensitivity of the

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

results to the particle size, thus we adopt same particle size and particle size distribution for all the samples. The angle range of parallel bonds been replaced is 720°. This angle range is believed to be a key factor affecting the degree of anisotropy21 and will be examined in Section 4. As can be noted from previous experimental results, the Brazilian tensile strength of anisotropic rocks varies from 1 to 2 MPa (e.g., Opalinus Clay) to more than 20 MPa (Asan Gneiss).16 DEM model can reproduce mechanical properties of both soft and hard rocks by modifying the micro-parameters. As the main purpose of this study is to model the behaviors of inherently anisotropic rocks, a dimensionless analysis was ﬁrst carried out to investigate the effect of weak layers properties. Therefore, we adopt the micro-parameters calibrated for the Lac du Bonnet Granite22 to generate the isotropic model (E¼ 78.6 GPa, UCS¼ 202.4 MPa) as listed in Table 1. The properties of smooth joint models can be derived from the parallel bonds based on Eqs. (1)–(7). The calculated parameters for smooth joint models are listed in Table 2. After the generation of anisotropic samples, the specimens are loaded by moving the top and bottom platens toward one another at a ﬁxed velocity. During the Brazilian test, the force acting on the platens is monitored and the maximum value is recorded (P). The Brazilian test is terminated when the force reaches 80% of the maximum value. After the test, the Brazilian tensile strength (BTS) can be calculated by

σt =

P πRt

(9)

where P is the maximum load, R is the radius of the Brazilian disc and y is the thickness (1 for the 2D model). However, this equation is obtained by assuming the rock material is isotropic and homogeneous.44 Since most anisotropic rocks fails in both tensile and shear modes,15,19 σt does not always equals to the actual tensile strength. Similar to Ref. 16, the equation will be mainly used for comparison purpose. The BTS of isotropic sample estimated from Eq. (9) is 35 MPa. Secant modulus calculated from 50% of peak tensile strength is adopted to evaluate the stiffness of the Brazilian discs under indirect tensile test. In Section 4, the effect of smooth joint parameters on the BTS, stiffness and crack mechanism is investigated systematically with

Table 1 Microscopic properties for isotropic rock sample22. Particle properties

Value

Bond properties

Value

Ec

62 GPa

Ec

62 GPa

k n/ k s

2.5

2.5

μ

0.5 1.66

k n/ k s σc τc

0.2 mm

λ

1.0

Rmax /Rmin

Rmin ρ

1577 36 MPa 1577 36 MPa

3169 kg/m3

49

the use of DEM. Bear in mind that the BTS and stiffness discussed in Section 4 are normalized by the results of isotropic sample.

4. Simulation results 4.1. Effect of strength of smooth joint The effect of smooth joint strength is investigated by reducing the values of normal strength (σc ) and cohesion (cb ) of the smooth joint with a factor of 1, 1/2, 1/3, 1/5 and 1/10. Other parameters are inherited from the parallel bonds. As can be seen in Fig. 3(a), the normalized BTS is signiﬁcantly affected by the strength of weak layers. When the strength of weak layers is close to that of rock matrix (parallel bonds), the sample is almost isotropic. Thus, the normalized BTS is approximately constant over the entire interval, which is similar to trend 1 proposed by Ref. 16. When the strength of weak layers is reduced by a factor of 1/2, a considerable decrease of BTS can be noticed especially when the anisotropy angle is high (75° and 90°), following trend 2. The BTS decreases constantly over the entire interval when the strength of smooth joint models reduced to 1/3 of the intact part (trend 3). However, when the strength of weak layers is reduced further to and below a certain value (1/5 of intact part), the variation of BTS exhibits the same trend, namely, decreases at low anisotropy angles (0–45°) and is followed by a leveling off which is very similar with the trend 4. Further reduction of the strength of smooth joints does not affect the trend anymore. Reducing the strength of smooth joint does not affect the stiffness before the reduction factor reaches 0.1. When σc : σ¯c = 1: 10, the relative stiffness stays constant at low anisotropy angles (0–30°) and decreases when θ > 45°. In the DEM model, micro-cracks can be classiﬁed into tensile and shear failure of parallel bonds and smooth joints. In this study, the numbers of cracks developed on parallel bond and smooth joint are counted and normalized by the total number of microcracks. When the strength of smooth joints is closed to that of the intact part, all specimens fail mainly by tensile failure of parallel bonds, as indicated by the square solid line in Fig. 3(c). The failure of smooth joints is negligible when the anisotropy angle is low (θ ≤ 30° ) and it starts to occur when the anisotropy angle exceeds 30° but only takes a small percentage (less than 20% even when the loading is applied parallel with weak layers). With the decrease of smooth joint strength, failure of weak layers becomes more dominant. When the strength is reduced by a factor of 1/3 (trend 3), the percentage of smooth joint cracks takes less than 20% percentage at low anisotropy angles (0° and 15°), but it increases dramatically to about 70% when θ ¼90°. Finally, when the strength of smooth joint is reduced to an extremely low value, the failure of weak layers is predominant even at low anisotropy angles (40% when θ ¼ 0°). 4.2. Effect of ratio between normal strength and shear strength

Table 2 Microscopic properties for smooth joint model inherited from parallel bond. ~250000

Normal stiffness, k¯ n [GPa/m] Shear stiffness, k¯ s [GPa/m]

~100000

Friction coefﬁcient, μ Dilation angle, ψ [deg] Tensile strength, σc [MPa] Cohesion, cb [MPa] Friction angle, φb [deg]

0.5 0 157 157 0

From Eq. (8), it can be seen that the shear strength of smooth joint model is dependent on the cohesion (cb ), friction angle (φb ) and normal stress (σ ) acting on it. The effect of ratio between cohesion (cb) and normal strength of smooth joint will be investigated by increasing cb and keeping σc ﬁxed at 15.7 MPa (which is 1/10 of intact part). Again, other properties of smooth joint models are inherited from parallel bonds. When the anisotropy angle is extremely low (0° and 15°) and high (75° and 90°), the effect of normal to shear strength (cohesion) of smooth joints is minor as can be seen in Fig. 4(a) and (b). In contrast, when the anisotropy angle is medium (30–60°), the effect of strength difference is more signiﬁcant. For instance, when the anisotropy

50

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

Fig. 3. Effect of smooth joint strength. Smooth joint strength is reduced by a factor of 1, 1/2, 1/3, 1/5, and 1/10 while keeping σc = cb (angle range: 7 20°). Variation of (a) normalized BTS, (b) normalized stiffness, and (c) percentage of micro-cracks (solid lines represent failure of parallel bond and dash lines represent failure of smooth joint).

Fig. 4. Effect of ratio between normal strength and cohesion of smooth joint. σc : σc = 1: 10; sc:cb ¼ 1:1, 1:2, 1: 3, 1: 5, and 1:10 (angle range: 7 20°). Variation of (a) normalized BTS, (b) normalized stiffness, and (c) percentage of micro-cracks (solid lines represent failure of parallel bond and dash lines represent failure of smooth joint).

angle is 45°, the relative BTS can vary from 0.7 to 1.0 as the ratio between cohesion and normal strength increase from 1 to 10. In addition, increase of relative stiffness between 45° and 90° can be noted in Fig. 4(b) with the increase of cohesion. When the anisotropy angle is low (0° and 15°), more micro-cracks

appear on the parallel bonds, thus the ratio between normal and shear strength (cohesion) of smooth joints does not affect the tensile strength of discs too much (Fig. 4(c)). When the inclination angle reaches 30°, the increase of shear strength of smooth joint prohibits the shear failure and force the failure to happen in parallel bonds. It is

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

more obvious at intermediate angles that the solid circle line is above the solid triangle and solid square lines. This explains the high BTS encountered at θ ¼ 30° and 45° when the shear strength is high.

51

When the inclination angle exceeds 45°, tensile failure of smooth joints becomes dominant, thus the increase of shear strength of smooth joints will not affect the BTS much.

Fig. 5. Effect of smooth joint stiffness.σc : σ¯c = 1 : 10 . The stiffness of smooth joint is reduced by 1, 1/3, 1/5, 1/10, 1/25 while keeping k¯ n = k¯ s (angle range: 7 20°). Variation of (a) normalized BTS, (b) normalized stiffness, and (c) percentage of micro-cracks (solid lines represent failure of parallel bond and dash lines represent failure of smooth joint).

Fig. 6. Effect of ration between smooth joint normal stiffness and shear stiffness. (k s = 1/25k n , kn/k s ¼1:1, 3:1, 5:1, 10:1 and 25:1) (angle range: 7 20°). Variation of (a) normalized BTS, (b) normalized stiffness, and (c) percentage of micro-cracks (solid lines represent failure of parallel bond and dash lines represent failure of smooth joint).

52

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

Fig. 7. Effect of ration between smooth joint normal stiffness and shear stiffness. (k¯ n = 1/25k n , k¯ s/kn ¼1:1, 3:1, 5:1, 10:1 and 25:1) (angle range: 7 20°). Variation of (a) normalized BTS and (b) normalized stiffness.

Fig. 8. Effect of friction angle. Smooth joint stiffness is inherited from parallel bond,σc : σ¯c ¼ 1:10 (angle range: 7 20°). Variation of (a) normalized BTS and (b) normalized stiffness.

Fig. 9. Effect of angle range. σc : σ¯c ¼1:10, other parameters of smooth joint are inherited from parallel bonds. Variation of (a) normalized BTS, (b) anisotropy ratio of BTS, (c) normalized stiffness, and (d) anisotropy ratio of stiffness.

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

53

Table 3 Micro-parameters calibrated for different types of rock. Microparameters Particle Parallel bond

Smooth-joint

Ec [GPa] E¯c [GPa) σ¯c [MPa] τ¯c [MPa] Angle range Normal stiffness, k¯ n [GPa/m] Shear stiffness, k¯ s [GPa/m] Tensile strength, σc [MPa] Cohesion, cb [MPa] Friction angle [φ ]

Postaer Sandstone

Boryeong Shale

Leubsdorfer Gneiss

Mosel Slate

22

33

66

110

22

33

66

110

177 4 177 4 7 20° 55,000

62 714 62 714 730° 21,000

857 20 857 20 7 35° 80,000

1007 23 1007 23 7 55° 34,000

55,000

21,000

80,000

34,000

17 17 0

5 50 0

10 40 30

7 7 0

Table 4 Young's modulus perpendicular and parallel to weak planes, determined on cylindrical specimens by uniaxial loading test. Rock

Postaer Sandstone Boryeong Shale Leubsdorfer Gneiss Mosel Slate

Experimental results16

Numerical results

E_0 (GPa)

E_90 (GPa)

E_0 (GPa)

E_90 (GPa)

– 23.5 62.2 49.0

27.9 41.4 85.0 101.8

27.8 24.3 61.7 49.8

28.5 40.8 85.5 100.9

4.3. Effect of smooth joint stiffness The effect of deformability of weak layers is studied in this part by reducing the stiffness of smooth joint models by different factors. Two scenarios are designed here: one is by reducing the shear and normal stiffness together from the original value of 250,000 GPa/m; the other is by increasing one stiffness value while keeping the other ﬁxed at 10,000 GPa/m. Both σc and cb are 15.7 MPa.

It is shown in Fig. 5(a) that the reduction of weak layer stiffness signiﬁcantly decreases the BTS at low anisotropy angles. However, the effect becomes negligible when θ ¼90° Similar effect can be noted on the variation of normalized stiffness in Fig. 5(b). This is due to the fact that the loading is mainly carried by the rock matrix. Thus, the reduction of smooth joint stiffness does not lower the BTS and stiffness much. In fact, it causes more failures of parallel bonds instead of smooth joints for all range of anisotropy angle as illustrated in Fig. 5(c). This analysis is further investigated by looking at the effect of normal and shear stiffness separately. In each case, one stiffness is assigned a low value as 10,000 GPa/m (1/25 of intact part) while the other stiffness is increased by 3, 5, 10 and 25 times. As shown in Fig. 6(a) and (b), both the normalized BTS and stiffness at low anisotropy angles decrease with the decrease of normal stiffness of smooth joint. The trends are very similar to the results of decreasing both normal and shear stiffness together (Fig. 5(a) and (b)). Possible explanation to this observation is the ease of closing of smooth joints when normal stiffness is low, which causes the

Fig. 10. Comparison of BTS from experiments and DEM simulations as a function of anisotropy angles. Experimental results are from Ref. 16. (a) trend 1, Postaer Sandstone; (b) trend 2, Boryeong Shale; (c) trend 3, Leubsdorfer Gneiss; and (d) trend 4, Mosel Slate.

54

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

Fig. 11. Tensile stress versus vertical displacement and increment of micro-cracks from simulated Boryeoung Shale with different anisotropy angles: (a) θ¼ 0°; (b) θ¼ 45°; (c) θ¼ 90°.

stress concentration on the parallel bonds. This explanation is supported by Fig. 6(c) where more breakages of parallel bonds are encountered at low anisotropy angles when the normal stiffness is low. However, the increases of shear stiffness of smooth joint do not vary the normalized BTS and stiffness much as shown in Fig. 7 (a) and (b). 4.4. Effect of friction angle The results of normalized BTS and stiffness as a function of anisotropy angles with different friction angles are compared in Fig. 8. It is clearly that the effect of friction angle on the mechanical response with different anisotropy angles is not as obvious as other parameters. As can be derived from Eq. (8), higher friction angle of smooth joints will result in higher shear strength. When θ ¼0° and 15°, the failure is dominated by the breakage of parallel bonds, so the change of friction angle in smooth joints will not affect the BTS. On the other hand, when θ ¼90° the failure is dominated by the tensile failure of smooth joints; therefore the change of friction angle does not affect the BTS. Only when the anisotropy angle is medium, the BTS and stiffness increase with higher friction angle but still very slightly. 4.5. Effect of angle range The angle range of parallel bonds being replaced by smooth joint models is thought to be a key factor controlling the degree of

anisotropy.21 The larger the angle range the more smooth joints will be included which represent heavily bedded rock. In this study, samples with six different angle ranges (7 10°, 715°, 720°, 725°, 730° and 735°) are generated with σc and cb ¼ 15.7 MPa. Simulation results are illustrated in Fig. 9. It is obviously illustrated in Fig. 9(a) that the larger the angle range, the more reduction of normalized BTS and stiffness especially towards high anisotropy angles. As presented in Fig. 9(b), the anisotropy ratio increases from 1.34 to 2.21 as the angle range changes from 7 10° to 735°. Meanwhile, the anisotropy ratio of stiffness increases from 1.07 to 1.35. Note that angle range of 7 20° is the same case as the trend 4 in Fig. 3(a). Tuning the angle range can allow us to control the anisotropy angle at which the BTS would level off. Therefore, this parameter can be calibrated to represent rocks with different degrees of anisotropy.

5. Calibration and validation Based on the results from Section 4, the following calibration procedures should be followed to reduce iterations: 1. The angle range which dictates the anisotropy ratio is ﬁrst selected. A reasonable value to start with is 710°. 2. The stiffness of parallel bond is calibrated to match Young's modulus when θ ¼ 90° as the effect of smooth joint stiffness the least in this direction. Both the parallel bond strength and smooth joint strength are set to a large value at this stage.

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

55

Fig. 12. Fracture patterns of simulated Boryeong Shale with different anisotropy angles and comparison with experimental results of Mancos Shale.45 (Blue lines: smooth joint contact model; red lines: failure of smooth joint; and black lines: failure of parallel bond). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

3. The stiffness of the smooth joint is reduced to match Young's modulus when θ ¼0°. The normal and shear stiffnesses are set equal to reduce free parameters. During this step, the Young’s modulus when θ ¼90° matched in step 2 may be also reduced. Therefore, a few iterations between steps 2 and 3 may be required. 4. The strength of parallel bond is calibrated to reproduce the BTS when θ ¼ 0°. θ ¼0° is chosen as the effect of smooth joint strength at this direction is minimal. Normal and shear strength of parallel bond are set equal for simplicity. 5. The strength of smooth joint (σc , cb and φ ) can be tuned to match different trends of BTS: for trend 1, the strength of smooth joint should be close to that of parallel bonds; for trend 2, cohesion of smooth joint should be reduced slightly (0.5–1 of intact part) while the normal strength should be reduced gradually to match the decrease of BTS when θ 445°; for trend 3, the strength of smooth joint can be reduced by a factor between 3 and 5; for trend 4, very low strength should be assigned to the smooth joint model (less than 0.1 of intact part). If the degree of anisotropy cannot be matched in step 5, the angle range should be increased and steps1–5 have to be repeated. Following the procedures discussed above, the numerical models are calibrated to represent four typical rocks exhibiting different trends, namely, Postaer Sandstone (trend 1), Boryeong Shale (trend 2), Leubsdorfer Gneiss (trend 3) and Mosel Slate (trend 4). The corresponding micro-parameters are listed in Table 3. Young's moduli perpendicular and parallel to weak layers determined from uniaxial compression test on cylinder samples are listed in Table 4. Good

agreement can be found between the numerical and experimental results in terms of stiffness and tensile strength as compared in Fig. 10. Further analyses are conducted to validate the numerical model by looking at the tensile stress–displacement curves and evolution of fractures. Fig. 11 illustrates the tensile stress calculated from Eq. (9) as a function of the vertical displacement of platen for simulated Boryeong Shale when θ ¼0°, 45° and 90°. As expected, the tensile stress increases linearly before peak and followed by a sudden drop upon brittle failure. Difference can be noted between the increments of micro-cracks, namely, the dominant failure mode switches from tensile failure of parallel bond when θ ¼ 0° (Fig. 11(a)) to tensile failure of smooth joint when θ ¼90° (Fig. 11(c)). The evolution of fractures at 80% pre-peak, peak, and 80% postpeak stages are illustrated in Fig. 12 together with the ﬁnal fracture patterns of Mancos Shale obtained in the laboratory.45 In general, the ﬁnal simulated macroscopic fracture patterns are consistent with those observed in the laboratory. When θ ¼ 0° the simulated failure mode is formed as tensile splitting of rock matrix. A few tensile failures of parallel bonds occur in the pre-peak stage. Meso-fracture starts to appear when the tensile stress reaches peak. Finally, a straight fracture forms along the vertical loading path in a brittle manner. It is worth noting that tensile failures of smooth joint originate from the edge of the disc in the post-peak stage, which is also observed in laboratory as secondary fractures.45,46 When θ ¼45°, the inﬂuence of smooth joints can be noted as cracks mainly initiate in terms of failure of smooth joints at pre-peak stage. Inclined fractures form at peak stage along the direction of weak layers. Then, several tensile fractures initiate in rock matrix and a major curved macro-fracture appears in a

56

K. Duan, C.Y. Kwok / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 46–56

same pattern with that observed from experiment sample. When θ ¼90°, independent tensile failures of smooth joints originate in the middle of specimen at pre-peak stage and grow towards the specimen boundaries with increase of loading. Parallel bonds between these failed weak layers start to connect at peak stage and ﬁnally penetrate the whole sample at post-peak stage. Again, this failure mode agrees well with those observed in laboratory on shales.12,45

6. Conclusions The behaviors of inherently anisotropic rocks under Brazilian test condition are investigated in this study based on DEM simulations. The effect of weak layers on the indirect tensile strength, stiffness, and fracture patterns is studied systematically. Simulation results are compared to previous experimental results and good agreement can be found. The DEM model can provide some insights on the relationship between the BTS and micro-parameters for different types of rock. By imposing individual smooth joints into the bonded particle model, the numerical approach proposed in this study can explicitly represent the presence of weak layers normally encountered in inherently anisotropic rocks. The behavior of the numerical model changes from nearly isotropic to heavily anisotropic with the reduction of smooth joint strength. Also, the ratio between shear strength and normal strength affects the failure pattern especially when the anisotropy angle is between 30° and 0°. The stiffness of weak layers will affect the strength especially when loading is applied perpendicular to them. Moreover, the normal stiffness plays a more signiﬁcant role. The angle range of parallel bonds being replaced by weak layers is a key parameter controlling the degree of anisotropy. The friction angle of weak layers plays a minor role on the mechanical behaviors. This numerical approach provides a new way to study the behaviors of inherently anisotropic rocks. This numerical model can be further used to investigate properties of inherently anisotropic rocks like fracture toughness and strength criterion under different stress conditions. Also, it can be employed in modeling engineering projects related to anisotropic rocks like borehole stability and propagation of hydraulic fracturing in shale gas exploitation in the future.

Acknowledgment The research was funded by the National Natural Science Foundation of China (NSFC) (Grant no. 51428902).

References 1. Everitt RA, Lajtai E. The inﬂuence of rock fabric on excavation damage in the Lac du Bonnett granite. Int J Rock Mech Min Sci 2004;41:1277–303. 2. ISRM. Suggested methods for determining tensile strength of rock materials. Int J Rock Mech Min Sci Geomech Abstr 1978;15:99–103. 3. Li D, Wong LNY. The Brazilian disc test for rock mechanics applications: review and new insights. Rock Mech Rock Eng 2013;46:269–87. 4. Andreev GE. A review of the Brazilian test for rock tensile strength determination. Part I: calculation formula. Min Sci Technol 1991;13:445–56. 5. Andreev GE. A review of the Brazilian test for rock tensile strength determination. Part II: contact conditions. Min Sci Technol 1991;13:457–65. 6. Steen BV, Vervoort A, Napier JAL. Observed and simulated fracture pattern in diametrically loaded discs of rock material. Int J Fract 2005;131:35–52. 7. Newman DA, Bennett DG. The effect of specimen size and stress rate for the Brazilian test—a statistical analysis. Rock Mech Rock Eng 1990;23:123–34. 8. Lekhnitskii SG. Anisotropic plates. Defense Technical Information Center; 1968. 9. Amadei B, Rogers J, Goodman R. Elastic constants and tensile strength of anisotropic rocks. In: Proceedings of the 5th International congress on rock mechanics; 1983. 10. Chen C-S, Pan E, Amadei B. Determination of deformability and tensile strength of anisotropic rock using Brazilian tests. Int J Rock Mech Min Sci 1998;35:43–61. 11. Barla G, Innaurato N. Indirect tensile testing of anisotropic rocks. Rock Mech 1973;5:215–30.

12. Cho JW, Kim H, Jeon S, Min K-B. Deformation and strength anisotropy of Asan gneiss, Boryeong shale, and Yeoncheon schist. Int J Rock Mech Min Sci 2012;50:158–69. 13. Tavallali A, Vervoort A. Effect of layer orientation on the failure of layered sandstone under Brazilian test conditions. Int J Rock Mech Min Sci 2010;47:313–22. 14. Tavallali A, Vervoort A. Failure of layered sandstone under Brazilian test conditions: effect of micro-scale parameters on macro-scale behaviour. Rock Mech Rock Eng 2010;43:641–53. 15. Debecker B, Vervoort A. Experimental observation of fracture patterns in layered slate. Int J Fract 2009;159:51–62. 16. Vervoort A, Min K-B, Konietzky H, et al. Failure of transversely isotropic rock under Brazilian test conditions. Int J Rock Mech Min Sci 2014;70:343–52. 17. Lisjak A, Tatone BS, Grasselli G, Vietor T. Numerical modelling of the anisotropic mechanical behaviour of opalinus clay at the laboratory-scale using FEM/DEM. Rock Mech Rock Eng 2014;47:187–206. 18. Lisjak A, Grasselli G, Vietor T. Continuum–discontinuum analysis of failure mechanisms around unsupported circular excavations in anisotropic clay shales. Int J Rock Mech Min Sci, 65; 2014. p. 96–115. 19. Debecker B, Vervoort A. Two-dimensional discrete element simulations of the fracture behaviour of slate. Int J Rock Mech Min Sci 2013;61:161–70. 20. Park B, Min K-B. Discrete element modeling of transversely isotropic rock. In: Proceedings of the 47th US rock mechanics symposium; 2013. 21. Kwok CY, Duan K, Tham LG. Numerical simulation of strength and deformation behavior of inherently anisotropic rocks. In: Proceedings of the 48th US rock mechanics symposium, Minneapolis; 2014. 22. Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329–64. 23. Wang B, Chen Y, Wong Tf A. discrete element model for the development of compaction localization in granular rock. J Geophys Res Solid Earth 2008;113 (B3):B03202. 24. Schöpfer Childs C. The orientation and dilatancy of shear bands in a bonded particle model for rock. Int J Rock Mech Min Sci 2013;57:75–88. 25. Cho N, Martin CD, Sego DC. Development of a shear zone in brittle rock subjected to direct shear. Int J Rock Mech Min Sci 2008;45:1335–46. 26. Park J, Song JJ. Numerical simulation of a direct shear test on a rock joint using a bonded-particle model. Int J Rock Mech Min Sci 2009;46:1315–28. 27. Duan K, Kwok CY, Tham LG. Micromechanical analysis of the failure process of brittle rock. Int J Numer Anal Meth Geomech 2015;39:618–34. 28. Scholtès L, Donzé F-VA. DEM model for soft and hard rocks: Role of grain interlocking on strength. J Mech Phys Solids 2013;61:352–69. 29. Schöpfer MPJ, Abe S, Childs C, Walsh JJ. The impact of porosity and crack density on the elasticity, strength and friction of cohesive granular materials: insights from DEM modelling. Int J Rock Mech Min Sci 2009;46:250–61. 30. Weng M-C, Li H-H. Relationship between the deformation characteristics and microscopic properties of sandstone explored by the bonded-particle model. Int J Rock Mech Min Sci 2012;56:34–43. 31. Fakhimi A, Alavi Gharahbagh E. Discrete element analysis of the effect of pore size and pore distribution on the mechanical behavior of rock. Int J Rock Mech Min Sci 2011;48:77–85. 32. Cundall P. A discontinuous future for numerical modelling in geomechanics? Proc ICE Geotech Eng 2001;149:41–7. 33. Itasca. PFC2D:particle ﬂow code in 2 dimensions. 4.0 ed.. Minneapolis: Itasca; 2010. 34. Koyama T, Jing L. Effects of model scale and particle size on micro-mechanical properties and failure processes of rocks—a particle mechanics approach. Eng Anal Bound Elem 2007;31:458–72. 35. Yoon J. Application of experimental design and optimization to PFC model calibration in uniaxial compression simulation. Int J Rock Mech Min Sci 2007;44:871–89. 36. Fakhimi A, Villegas T. Application of dimensional analysis in calibration of a discrete element model for rock deformation and fracture. Rock Mech Rock Eng 2007;40:193–211. 37. Ivars DM, Pierce ME, Darcel C, et al. The synthetic rock mass approach for jointed rock mass modelling. Int J Rock Mech Min Sci 2011;48:219–44. 38. Cundall P, Potyondy D, Lee C. Micromechanics-based models for fracture and breakout around the mine-by tunnel. In: Martino JB, Martin CD, editors. Proceedings of International conference on deep geological disposal of radioactive waste, Winnipeg. Toronto: Canadian Nuclear Society; 1996. p. 113–22. 39. Scholtès L, Donzé F-V. Modelling progressive failure in fractured rock masses using a 3D discrete element method. Int J Rock Mech Min Sci 2012;52:18–30. 40. Zhang Q, Zhu HH, Zhang LY, Ding XB. Study of scale effect on intact rock strength using particle ﬂow modeling. Int J Rock Mech Min Sci 2011;48:1320–8. 41. Esmaieli K, Hadjigeorgiou J, Grenon M. Estimating geometrical and mechanical REV based on synthetic rock mass models at Brunswick Mine. Int J Rock Mech Min Sci 2010;47:915–26. 42. Bahaaddini M, Sharrock G, Hebblewhite B. Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a non-persistent jointed rock mass under uniaxial compression. Comput Geotech 2013;49:206–25. 43. Chiu C-C, Wang T-T, Weng M-C, Huang T-H. Modeling the anisotropic behavior of jointed rock mass using a modiﬁed smooth-joint model. Int J Rock Mech Min Sci 2013;62:14–22. 44. Jaeger JC, Cook NG, Zimmerman R. Fundamentals of rock mechanics. 4th ed.. Oxford: Wiley; 2009. 45. Simpson NDJ, Stroisz A, Bauer A, Vervoort A, Holt RM. Failure mechanics of anisotropic shale during Brazilian tests. In: Proceedings of the 48th US rock mechanics symposium, Minneapolis, Paper ARMA 2014-7399. 46. Tan X, Konietzky H, Frühwirt T, Dan D. Brazilian tests on transversely isotropic rocks: laboratory testing and numerical simulations. Rock Mech Rock Eng 2014:1–11.