A comprehensive parametric study of grain-based models for rock failure process simulation

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

Contents lists available at ScienceDirect

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

A comprehensive parametric study of grain-based models for rock failure process simulation

T

Xin Wanga, Ming Caia,b,c,

⁎

a

MIRARCO – Mining Innovation, Laurentian University, Sudbury, Canada P3E 2C6 Bharti School of Engineering, Laurentian University, Sudbury, Canada P3E 2C6 c Key Laboratory of Ministry of Education for Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110004, China b

ARTICLE INFO

ABSTRACT

Keywords: Grain-based model 3DEC-GBM Parametric study Rock failure modeling

Grain and contact are the two key components in 3DEC-GBMs (grain-based models) because they control the micro-mechanical behavior and consequently the macro-mechanical behavior of rock. Three types of grains – rigid, elastic, and breakable grains – are considered in this study to explore the influence of grain properties on the mechanical behavior of rock. Equivalent 2D plane strain problems are solved using 3DEC to reduce computation time. A non-uniform distribution of grain size is used to generate numerical models using Neper, a software package for polycrystal generation and meshing. A comprehensive parametric study of properties of contact and grain of the three types of grains is conducted and calibration procedures are suggested for numerical modeling. The modeling results indicate that both the contact and the grain properties affect the macroscopic mechanical behavior of synthetic rock specimens. All the three types of grains produce reasonable results for pre-peak deformation (E and υ) and macroscopic strength (UCS, σt, C and ϕ) parameters. Each grain type has advantages and disadvantages in numerical modeling. Rigid grains cannot produce volumetric strain properly and a near-zero residual friction angle of contact is needed to capture post-peak deformation behavior due to the grain interlocking effect. Elastic grains also need a very low residual friction angle of contact and brittle rock deformation behavior cannot be captured. The complex calibration procedure and long computation time are the main issues of breakable grains, although it simulates well the volumetric strain, does not require a very low contact residual friction angle, and captures brittle rock behaviors better than other two types of grains. The findings from this study are useful for rock failure process simulation using the 3DEC-GBM modeling approach.

1. Introduction Brittle rock failure is a gradual deformation process that involves microcrack closure, initiation, propagation and coalescence, which eventually leads to macro-fracture formation in rock.1–6 It is widely accepted that grain-scale microstructures control the micro-mechanical behavior of the grain itself and hence the macro-mechanical behavior of rock.1–9 As can be seen from Fig. 1, different mineral types with random polygonal shapes and various mineral sizes are identified at grain scale in a granite sample. A proper representation of grains (such as shape, size, and mineral types) and contacts (such as stiffness anisotropy and strength variation) at the grain scale is essential for capturing brittle rock failure process using numerical models. Numerical modeling has emerged as a powerful tool to study mechanical behaviors of rock, rock-like materials and jointed rock masses

⁎

at both the laboratory and the field scales. Considerable efforts have been made to model macroscopic mechanical behaviors of rock based on grain-scale consideration of grain characteristics such as grain shape, grain size distributions, and grain contact heterogeneity. Indirect (continuum) and direct (discontinuum and hybrid continuum-discontinuum) numerical methods can be used to conduct rock failure process modeling.11 In the continuum modeling approach, rock damage is represented by a constitutive relation with a failure criterion in the form of elastic-plastic, elastic-brittle, or strain-softening models. Finite Element Method (FEM), Finite Difference Method (FDM) and Boundary Element Method (BEM) are the main continuum modeling methods. PLAXIS 2D & 3D,12 RS2 and RS3,13 ABAQUS,14 FLAC 2D and 3D, Map3D16 and EXAMINE3D13 are the most popular continuum analysis software packages. Cundall17 pointed out that phenomenological and inappropriate representation of progressive failure of rock are the two

Corresponding author at: MIRARCO – Mining Innovation, Laurentian University, Sudbury, Canada P3E 2C6. E-mail address: [email protected] (M. Cai).

https://doi.org/10.1016/j.ijrmms.2019.01.008 Received 5 July 2018; Received in revised form 24 January 2019; Accepted 26 January 2019 1365-1609/ © 2019 Published by Elsevier Ltd.

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 1. Grain scale microphotographs of (a) LdB granite and (b) Forsmark granite. Four minerals are denoted by Qtz (Quartz), Plag (plagioclase), Kfsp (potassium feldspar), and Bt (biotite). Image size is 2.78 × 2.10 mm2.10.

major limitations of the indirect modeling scheme. Compared with the indirect modeling scheme, the direct modeling scheme represents a rock as an assembly of discrete particles or blocks and rock failure progress is modeled directly without using pre-defined constitutive models. Discrete Element Method (DEM), hybrid, and combined Finite Element and Discrete Element Method (FEM/DEM), Discontinuous Deformation Analysis (DDA), Numerical Manifold Method (NMM) are the main direct modeling methods. Particle Flow Code (PFC2D & 3D), Universal Distinct Element Code (UDEC & 3DEC)15 and YADE18 are the most popular DEM codes. ELFEN,19 Y-code,20–22 and IRAZU23 are the most popular FEM/DEM (or FDEM) codes. The direct modeling method is widely used to model mechanical behaviors of rock under various loading conditions such as uniaxial compression, triaxial compression, direct and indirect tensions, shearing, and thermal loading using the commercial or open source codes mentioned above.

procedure and the advantages and disadvantages of using different types of grains in numerical modeling using GBMs. In this study, a brief introduction of GBM and Voronoi tessellation in Neper is provided in Section 2. Micro- and macro-mechanical properties of Lac du Bonnet (LdB) granite and the generation of 3DEC-GBMs are presented in Section 3. Next, a comprehensive parametric study of micro-properties of contacts and grains is conducted, followed by presenting calibration procedures for 3DEC-GMBs modeling. Finally, numerical modeling results of calibrated 3DEC-GBMs with different types of grains are presented and discussed.

In numerical modeling, different grain shapes are used to consider rock heterogeneity at the grain scale. For instance, disk and sphere,11 square and cube,26 triangular and tetrahedron1,4,7,21,27 and polygonal and polyhedron3,23,29,37,38 shapes have been used. As can be seen in Fig. 1, polygonal and polyhedron grain structures seem to be more realistic representations of grain structures.3 Voronoi tessellation is one of the widely used techniques to generate polygonal or polyhedron grains for microstructure modeling of materials.3,7,38–41 DEM-Voronoi models have been implemented in UDEC to model rock as a pack of polygonal-shaped grains, which are bonded along contacts.42 In addition, the grain-based model (GBM) has been made available in UDEC to better represent rock as a polygonal grain structure.3,7,43 However, external tools are needed to construct GBMs in 3DEC. For instance, Ghazvinian et al.38 constructed 3D Voronoi GBM models using Neper44 and simulated brittle rock failure in 3DEC. Gao and Stead45 proposed a modified Trigon logic method to model brittle rock failure in both 2D and 3D. The rock specimen is assembled by cutting a conventional polygonal or polyhedral block into smaller triangular or tetrahedral blocks. A polyhedron grain in a GBM model can be rigid or deformable, elastic or plastic. For simplicity and computation efficiency, polygonal grains are usually considered as unbreakable (rigid or elastic) (i.e., Farahmand and Diederichs7 and Ghazvinian et al.38) and the overall rock behavior is controlled by the contact's behavior for rigid grains. In GBM modeling, micro-properties of grains and contacts are not the same for different types of grains. Each grain consideration has advantages and disadvantages in terms of computation time and complexity of model calibration. To the authors’ best knowledge, there is no comprehensive study that thoroughly discusses the model calibration

The DEM was developed by Cundall46 to simulate movement and deformation of irregularly shaped particles in geotechnical problems. One of the most important applications of the DEM is simulating rock deformation behaviors by treating a rock as an assemblage of circular or spherical particles that are cemented at the contact points. This is called bonded particle model (BPM) and it has been implemented in PFC2D & 3D and used to study rock deformation behaviors at laboratory and field scales.11,47 The grain-based model (GBM) was originally developed in PFC2D by Potyondy4 to address the issue of the unrealistic assumption for rigid disk-shaped grains. Before the GBM modeling approach was developed, clustering and clumping techniques were developed to address the issue of grain shape. Grains are breakable in clustered particle models (ClsPM) (Potyondy and Cundall11). Grains are assumed as rigid and unbreakable in clumped particle models (ClmPM) from Cho et al.30 Because grains can fail (i.e., grain crushing) in reality, the assumptions of rigid and unbreakable grains and that only contact failure controls rock failure are not realistic. Cho et al.30 summarized three main disadvantages of the conventional PFC modeling approach: high tensile (σT) to uniaxial compressive strength (UCS) ratio, poor strength envelope (low confined strength), and a very low friction angle of the macroscopic strength parameter. Scholtès L and Donzé48 proposed a method to implicitly consider grain interlocking for solving the issues of low UCS/σT ratio and poor strength envelope in DEM modeling. In the PFC-GBM, a rock is modeled as a synthetic material with fully packed deformable, breakable, polygonal grains and breakable grain contacts. Details of the grain structure generation procedure are described by Potyondy4. Gao et al.1 mentioned that an intrinsic deficiency of the PFC-GBM modeling approach is the issue of porosity and it is difficult to model low porosity rocks due to the circular shape of the particles.

2. Brief Introduction of GBM and Neper 2.1. GBM

i.e.,1,4,8,11,20,21,24–36

61

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 2. Comparison of grain shapes and fracture development paths in 3DEC-GBM. (a) Tetrahedron grains generated from 3DEC; (b) polyhedron grains generated from Voronoi tessellation in Neper.9.

GBM has been developed in UDEC using the UDEC-Voronoi modeling approach.1,3,7,45 A rock is modeled as a pack of polygonal-shaped interlocked blocks of convex topology that are bonded along the contacts. The UDEC-GBM modeling approach can overcome some limitations of the PFC-GBM modeling approach and better capture Poisson's ratio, σT/UCS ratio, and failure envelope due to the full contact between grains and a better interlocking resulted from polygons of different shapes and sizes. On the other hand, if grain breakage is not considered, the strong interlocking induced by the polygonal-shaped blocks can cause undesired model response, particularly in the post-peak deformation stage. In some studies,7,30 very low contact friction angles have to be used to off-set the effect of block interlocking in order to match the laboratory test results. Rock also can be modeled as a 3D assemblage of grains that are bonded along faces of grains using 3DEC-GBM. Grain faces are modeled as discontinuities acting as distinct boundaries interacting between blocks. In this manner, each grain can detach, rotate, and slide against one another. In 3DEC-GBMs, tetrahedral and Voronoi tessellations can be used to generate grains (as shown in Fig. 2). Few studies have adopted the 3DEC-GBM modeling approach to simulate rock failure due to the high demand on computing power. A weak coal was simulated using a Trigon logic tessellation in 3DEC by Gao and Stead45 Brittle rock damage was simulated using polyhedron grains in 3DEC by Ghazvinian et al.38 Brittle rock failure with the consideration of inter- and intra-grain contact failures was studied using a 3DEC-GBM modeling approach.9 Brittle rock fracturing processes were investigated using both tetrahedral and polyhedron grains by Azocar.34 Due to different morphology and arrangement between tetrahedral (smooth surface) and polyhedron (rough surface) grains, different modeling results are expected. For example, a relatively smooth failure path is expected due to more shear contact failures when tetrahedral grains are used (Fig. 2(a)). On the other hand, a rough failure path can be expected due to tensile failures of dominant contacts when polyhedron grains are used (Fig. 2(b)).

generate tessellations: one is from cell morphological properties, such as grain size distribution, grain shape distribution and detailed grain data obtained from laboratory observation; the other is from multi-scale tessellations. Multi-scale tessellations are realized by subdividing the cells of a primary tessellation into a secondary tessellation and combining them into one tessellation in the end. More details about the other two modules of Neper and the application of Neper can be referred to some references.49–51 In this study, Neper is used to generate grains based on grain size distribution. 3. Lac du Bonnet (LdB) Granite and 3DEC-GBM Generation 3.1. Micro- and macro-mechanical properties of LdB granite The well-studied and documented Lac du Bonnet (LdB) granite, which is the main rock type at AECL's Underground Research Laboratory (URL) in Manitoba, Canada, is chosen to conduct the parametric study. The mineral composition of LdB granite is 42% feldspar, 29% quartz, 22% plagioclase and 6% biotite. The overall mean size of the four minerals is about 2.0 mm from the grain size distribution. It should be mentioned that a proper GBM should produce a statistical grain size distribution similar to that from real rock. Based on the published mineral experiment data and other research works (e.g., Bass52 and Lan et al.3), micro-material properties for the four minerals of LdB granite are summarized in Table 1. Fig. 3 shows the nonlinear Hoek-Brown and the linear Mohr-Coulomb strength envelopes associated with the laboratory compression test data. For calibration purpose, the macro-mechanical properties of LdB granite are summarized in Table 2. 3.2. 3DEC-GBM generation It is important to obtain a reasonable representation of the grain Table 1 Micro-material properties of the four minerals of LdB granite (from Bass52 and Lan et al.3).

2.2. Brief introduction of Neper GBMs have been adopted in various DEM codes and a proper generation of grain structures is required in numerical modeling. Neper is a robust open source package for polycrystal generation and meshing in 2D and 3D.44 It can be compiled and run on any Unix-like operating system. The tool is mainly composed of three modules: generation, meshing and visualization. Neper generates polycrystals as Voronoi tessellations with convex-shaped cells in a space domain (i.e., cuboidal-, cylindrical- and spherical-shaped domains). There are two ways to

Mineral type

ρ (kg/m3)

E (GPa)

υ

Grain size (mm)

Abundance (%)

Feldspar Quartz Plagioclase Biotite/chlorite Average

2560 2650 2630 3050 2605

69.8 94.5 88.1 33.8 78.1

0.28 0.08 0.26 0.36 0.22

3.00 1.50 3.50 0.75 2.50

42 29 22 6

Notes: ρ means density; E means Young's modulus; υ means Poisson's ratio. 62

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

computation time and to simplify the models, 3DEC is used to study plane strain problems (similar method is adopted by Guo et al.35 and Wang and Cai9). The displacements are fixed in the Y direction to realize the plane strain state and the thickness of all models are 1 mm along the Y direction (Fig. 4(b)). 3.3. Model set-up For UCS test modeling, axial load is applied via a constant vertical velocity (Vz) on the top of the model while the bottom of the model is fixed in the vertical direction (Fig. 4(b)). The vertical velocity is applied directly on the grid points of the top surface of the model. This is an ideal loading condition (frictionless loading) that does not cause the end effect induced by loading platens. It should be noted that the model is clamped by two rigid loading plates at both ends for the case of rigid grains, which is not shown in Fig. 4(b). It is known that the loading velocity has a large influence on the modeling results. The propagation speed of the loading disturbance, which is a dynamic process that is implemented by a time-stepping algorithm in 3DEC, depends on the physical properties of the block system. A sensitivity analysis of the loading velocity on the modeling results was undertaken and a proper loading velocity was determined to ensure that the model remains in quasi-static equilibrium. It is found that for the given specimen geometry and material properties, a constant velocity of 0.1 m/s is a proper loading velocity to conduct the UCS test modeling. For the direct tension test simulation, a smaller vertical velocity (0.01 m/s) is applied on the top of the model in the reversed Z direction. For the triaxial compression test simulation, confinements are applied on both the X-Z and the Y-Z planes (Fig. 4(b)). It should be mentioned that confinement on the X-Z plane has no influence on the overall model response in the modeling because the displacement is fixed in the Y direction (to ensure the plane strain condition). Axial and lateral strains are monitored by multiple pairs of monitoring grid points and the strains are calculated by averaging grid displacements using a FISH script (FISH is a programming language embedded within 3DEC that enables users to define new variables and functions). As shown in Fig. 4(b), seven and thirteen pairs of grid points are arranged 2 mm away from the model boundaries to monitor the axial and lateral strains, respectively. It should be noted that because the model is simplified based on the plane strain condition, the Poisson's ratio calculated from the recorded strains refers to the plane strain Poisson's ratio in this study. For the case of rigid grains, the axial stress is calculated based on the vertical force and the surface area of the rigid loading platens. For the cases of elastic and breakable grains, the axial stress is measured using two methods in order to check the loading velocity: one is to monitor the reaction forces of sub-contacts on the top surface of the model and the other is to average the zone stress σzz for all blocks using a FISH script. Trial test runs were performed to check the validation of the methods. The stress–strain curves derived from the two methods overlap for the adopted loading velocity, meaning that the loading velocity is appropriate for the numerical modeling.

Fig. 3. Laboratory compression test data and strength envelopes of the nonlinear Hoek-Brown and the linear Mohr-Coulomb failure criteria (test data from Martin53).

microstructures with topological and statistical properties similar to that of real rock,4 which is important for considering geometric and material heterogeneities at the grain-scale in 3DEC-GBM modeling. It is widely accepted that generating a non-uniform grain size distribution is deemed imperative but difficult. In UDEC, only a near-uniform grain size distribution is produced by the Voronoi generator. In 3DEC, there is no available built-in generator and a third-party software is needed to generate the grains. The open-source software Neper is used to generate grain structures here. Based on the grain size distribution of LdB granite, different distributions such as lognormal, normal and Weibull distributions can be used to fit the grain size distribution. In Neper, a normal distribution (average grain size μ = 2 mm, standard deviation σ = 1 mm with cutoffs at 0.6 and 6 mm to eliminate very small and large grains) is used to generate the grain structures within a model dimension of 50 mm × 100 mm × 1 mm, which contains 1715 grains with equivalent grain sizes from 0.6 to 6 mm, as shown in Fig. 4(a). The cutoffs are necessary because very small grains increase the total numbers of grains and contacts and hence the computation time significantly. Once the grains are generated in Neper, a file containing the information of the grains is saved in a format (e.g., blocks.3ddat) that can be imported to 3DEC directly. It is impractical to simulate a full 3D model as that used in the laboratory test with all the grain details to mimic reality. To reduce the Table 2 Macroscopic mechanical properties of LdB granite.53. Property

Value

Deformation properties

Stress levels

Strength properties

E (GPa)

υ

σci (MPa)

σcd (MPa)

ϕ* (°)

C* (MPa)

σT (MPa)

UCS (MPa)

69 ± 5.8

0.26 ± 0.04

82 (41% UCS)

172 (86% UCS)

55.6

30.9

9.3 ± 1.3

200 ± 22

Notes: σci and σcd denotes crack initiation stress and crack damage stress, respectively; σT means tensile strength; UCS means uniaxial compressive strength; it should be mentioned that C* (cohesion) and ϕ* (friction angle) of the Mohr–Coulomb failure criterion are estimated from the laboratory triaxial compression test data with a confinement pressure range of 0–40 MPa for linear failure criteria (the maximum confinement pressure was 60 MPa in the laboratory triaxial compression tests). 63

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 4. (a) 3DEC-GBM model generated in Neper using grain size distribution of LdB granite; (b) loading conditions and monitoring points for UCS and triaxial compression test simulations. Note: the thickness is 1 mm for all modeling cases.

4. Comprehensive parametric study

process modeling. As mentioned in Section 1, a polyhedron grain can be rigid or deformable, elastic or plastic. The use of each type of grain has its advantages and disadvantages in terms of computation time and model calibration. To deepen our understanding of the role of grains in 3DECGBM and to shed light on the advantages and disadvantages of the types of grains used in numerical modeling, three types of grains are considered herein. As can be seen in Fig. 5, grains are treated as rigid (Case 1), elastic (Case 2) and breakable (Case 3) grains in the modeling. For simplicity, Case 1, Case 2, and Case 3 are denoted as C1, C2, and C3, respectively hereafter. Three grain cases (C1, C2, and C3) are realized through various

4.1. Parametric study considerations As shown in Fig. 5, the grains (blocks) and the contacts (joints) are the two key components that control the micro- and macro-mechanical behaviors of a 3DEC-GBM synthetic specimen. A common practice of model simplification is setting grains as a rigid or an elastic material.7,37,38 It should be mentioned that in most GBM modeling studies, the grain itself is assumed elastic and only grain contact failure is allowed. It is known that grains can also fail in shear or in tension. Therefore, grains should be breakable for an improved rock failure

Fig. 5. Three cases of grain considerations for parameter calibration in 3DEC-GBM. 64

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

considerations of the blocks such as rigid and deformable blocks. Deformable blocks can be assigned different material models such as linear elastic, Mohr-Coulomb, anisotropic elastic, Hoek-Brown and strain-softening plasticity models. For C1, blocks are rigid and the overall macroscopic behavior of a specimen is determined by the contact behavior. For C2, blocks are deformable and a linear elastic isotropic constitutive model is assigned to the blocks. In this model, stress increments are derived based on strain increments according to Hooke's law:15

= 2G

ij

ij

2 G 3

+ K

k , k ij

to their residual values. In addition, contact dilation can also occur when the contact is in the shear state. More details about the governing constitutive relations of the contacts can be found in the 3DEC manual.15 As shown in Fig. 5, micro-properties of contacts determine the overall macroscopic behavior of rock for C1. For C2 and C3, both the micro-properties of grains and the contacts affect the macroscopic behavior of rock. In other words, contact properties are important and are required for the three cases. For the Coulomb slip joint model, microparameters of deformation (normal and shear stiffness) and strength (contact cohesion, friction angle and tensile strength) are needed for model calibration. For grains in C1, only microscopic deformation parameters are needed to run the models. For grains in C2, microscopic deformation parameters are needed and the macro-deformation behavior is influenced by both the grains and the contacts. The overall strength behavior is controlled only by the contacts in C2. For grains in C3, both the micro-properties of grains and contacts contribute to the overall mechanical behavior of a specimen. To obtain a deep understanding the role of micro-parameters of grains and contacts on the macro-mechanical behavior of the SRMs and to develop calibration procedures for each simulation case, a comprehensive parametric study is conducted. The objective of the parametric study is to establish a representative range of grain-scale micro-parameters that can reproduce the desired macro-scale properties measured in laboratory tests. Based on the parameters listed in Tables 1 and 2 and the numbers in literature,1,7,54–56 Table 3 presents the reference grainscale input parameters for the grains and the contacts in the parametric study. It should be mentioned that only one type of grain and one type of contact are considered in the parametric study for each case. For simplicity, the residual cohesion and tensile strength of the grains are assumed to be 5 and 0 MPa for C3. The residual cohesion and tensile strength of the grain contacts are also assumed to be zero, which are normal practice in GBM modeling.

(1)

where ij and ij are the stress and strain increment tensors, respectively. K and G are the bulk and shear modulus, respectively, and ij is the Kroenecker delta. The new stress values ( ijN ) are determined from Eq. (2): N ij

=

ij+

(2)

ij

For C3, a strain-softening model is used to make the blocks deformable and breakable. The strain-softening model is based on the Mohr-Coulomb model with non-associated shear and associated tension flow rules. In 3DEC, the yield and potential functions, plastic flow rules and stress corrections the strain-softening model are identical to those of the Mohr-Coulomb model. In the strain-softening model, model parameters such as cohesion, friction angle, and tensile strength are softened after the onset of plastic yield, while these parameters are assumed to be constant in the Mohr-Coulomb model. These model parameters (i.e., cohesion, friction angle, and tensile strength) can be defined by piecewise-linear softening functions of the characteristic plastic shear strain. More details about the strain-softening model can be found in the 3DEC manual.15 In addition, the Coulomb slip contact model is used for the contacts for all the three cases. For deformation behavior of contacts, stress–displacement relations are assumed to be linear in the elastic range and are governed by the contact's normal and shear stiffness:15 n, c sh, c

=

JKn un, c

(3)

=

JKs ush, c

(4)

4.2. Influence of JKn/JKs on deformation properties According to previous studies,11,24,30,37,38,45 the Poisson's ratio is governed by the normal (JKn) to shear (JKs) stiffness ratio of contacts. The Young's modulus (E) and the Poisson's ratio (υ) of the numerical specimens are normalized to the macroscopic deformation parameters from the laboratory tests (E = 69 GPa and υ = 0.26, Table 2). Fig. 6

where Δσn,c and Δσsh,c are the effective normal and shear stress increments, respectively, JKn and JKs are the contact's normal and shear stiffness, respectively, and Δun,c and Δush,c are the normal and shear displacement increments, respectively. Tensile and shear states are the two main considerations for the strength of contacts. In the tensile sate, a limiting tensile strength (i.e., JσT) is assumed for the intact joint, which has no previous slip or separation. The contact breaks in tension once its tensile strength is reached to the limiting tensile strength; then σn,c is dropped to zero and the contact is marked as a tensile crack. Similarly, in the shear state, shear stress is determined by the following considerations. Thus, if

JC +

sh, c

n tanJ

=

sh, max

Table 3 Reference micro-parameters for grains and contacts in the parametric study. Grain block properties Deformation parameters (C1 C2 C3)

(5)

Strength parameters (C3)

then sh, c

e = JKs ush ,c

(6)

sh, c

(7)

sh, max

then, the σsh,c is dropped to the residual shear stress: sh, c

= JCr +

n tanJ r

=

sh, r

ρ (kg/m3)

2605

E (GPa)

78.13

υ Peak grain strength

0.22 Cp (MPa) ϕp (°) σTp (MPa) ψp (°) Cr (MPa) ϕr (°) σTr (MPa) ψr (°)

62 65 29 26 5 31 0 22

εcps (%) εfps (%)

0.1 0.15

Grain strength when strain is greater than the characteristic plastic strain Plastic strain

or else, if

(8)

where JC and Jϕ are the peak cohesion and friction angle of the contact, e ush , c is the elastic component of the incremental shear displacement, JCr and Jϕr are the residual cohesion and friction angle of the contact. A contact slips and forms a shear crack once the maximum shear stress is reached and the corresponding cohesion and friction angle are reduced

Contact properties JKn (GPa/m) JKs (GPa/m) JKn/JKs Peak contact strength

45000

2.5 JCp (MPa) Jϕp (°) JσTp (MPa)

60 60 20

Residual contact strength

JCr (MPa) Jϕr (°) JσTr (MPa)

0 31 0

18000

(C1 C2 C3)

Notes: JKn and JKs denote normal and shear stiffness of grain contacts, respectively; C and JC denote cohesions of grains and contacts, respectively; ϕ and Jϕ denote friction angles of grains and contacts, respectively; σT and JσT denote tensile strengths of grains and contacts, respectively; ψ means dilation angle; εps means characteristic plastic strain. 65

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 6. Effect of contact JKn/JKs ratio on E and υ for three cases: C1 with JKn = 45,000 GPa/m (a), C2 with JKn = 250,000 GPa/m (b) and C3 with JKn = 380,000 GPa/m (c).

shows the influence of the JKn/JKs ratio on the normalized E and υ when JKn is kept at 45,000 (C1), 250,000 (C2) and 380,000 (C3) GPa/ m. For all the three grain cases, it is found that the JKn /JKs ratio influences both E and υ but within different extents. As the JKn/JKs ratio increases from 1 to 10, E decreases while υ increases. As the JKn/JKs ratio increases by keeping JKn constant, the model becomes softer due to the decrease of JKs, which leads to more lateral deformations and hence large Poisson's ratios. For the case of rigid grains (C1), a relatively low JKn (45,000 GPa/m) is required to match the macroscopic E and υ compared with the values used in C2 (250,000 GPa/m) and C3 (380,000 GPa/m). Among the three grain cases, the macroscopic E and υ of C1 show a high sensitivity to the JKn/JKs ratio, which demonstrates that only the contact's deformation properties determine the overall macroscopic deformation properties. For C2 and C3, a higher JKn is required to match the macroscopic deformation parameters due to the influence of the deformability of both the grains and the contacts on deformation.

for C2 and C3. It is seen that E increases as JKn (also JKs) increases. A high increasing rate of E is found as JKn increases from 10,000 to 240,000 GPa/m for C2 and from 80,000 to 400,000 GPa/m for C3. As the contact stiffness further increases, the macroscopic E is approaching the Young's modulus of the micro-grains (78 GPa, see Table 3). On the other hand, the influence of JKn on υ is not notable compared with its influence on E. As the contact stiffness increases, the normalized macroscopic υ increases from 0.85 to 1.11 for C2 (Fig. 7(b)) and from 0.97 to 1.05 for C3 (Fig. 7(c)). Comparing the results shown in Fig. 7 and Fig. 6, it is seen that the υ value depends mostly on the JKn/JKs ratio, which indicates that the Poisson's ratio should be the first macroproperty parameter to be calibrated. 4.4. Influence of JCp/JσTp on macroscopic σT As mentioned in previous studies, the strengths of grain contact control directly the overall strength of the model and the post-peak deformation behavior.11,24,45 Two scenarios are considered to investigate the effect of the JCp/JσTp ratio on the tensile strength of a specimen using direct tensile loading. The first scenario is fixing the contact's peak cohesion JCp at 100 (C1) and 60 (C2 and C3) MPa and varying the contact's tensile strength JσTp from 10 to 50 MPa (C1) and 6–60 MPa (C2 and C3). The second scenario is fixing JσTp at 21.5 (C1) and 14 (C2 and C3) MPa and changing JCp from 43 to 215 MPa (C1) and from 14 to 140 MPa (C2 and C3), which results in JCp/JσTp ratios ranging from 2 to 10 for C1 and 1–10 for C2 and C3. Fig. 8 shows the modeling results with the macroscopic tensile strengths normalized to the laboratory tensile strength (9.3 MPa, see Table 3). Similar trends are found for all the three cases. For the cases with identical contact peak cohesion (JCp = 100 MPa for C1 and 60 MPa for C2 and C3), the macroscopic direct tensile strength decreases as the JCp/JσTp ratio increases. As expected, the macroscopic

4.3. Influence of JKn and JKs on deformation properties The Young's modulus (E) is closely related to contact's normal and shear stiffness, as noted in some works.37,38,45 Fig. 7 presents the influence of JKn on E and υ when the JKn/ JKs ratio is kept at 4 (C1) and 2.52 (C2 and C3). In general, E increases as JKn (also JKs) increases for all the three cases although there is a difference between C1 and C2 (C3). The model stiffness increases due to the increase of JKn and JKs. In addition, different trends of variation of υ can be seen in the three simulation cases. For C1 (rigid grains), the macroscopic E increases linearly as JKn increases from 10,000 to 120,000 GPa/m, while the macroscopic υ is independent of JKn and JKs (Fig. 7(a)); similar results were reported by Kazerani and Zhao.37 In addition, similar trends are found for E and υ 66

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 7. Effect of contact JKn on E and υ for three cases: C1 (a), C2 (b) and C3 (c).

tensile strength is strongly related to the grain contact's tensile strength, which directly controls the overall tensile strength of the specimen. For the cases with the same contact tensile strength (JσTp = 21.5 and 14 MPa), it is found that the JCp/JσTp ratio has a small influence on the macroscopic tensile strength for all the three cases. Compared with the grain contact's cohesion, the grain contact's tensile strength plays a key role in determining the macroscopic tensile strength. Hence, the macroscopic tensile strength should be the first macroscopic strength parameter to be calibrated after the deformation parameters are calibrated. It should be mentioned that the JCp/JσTp ratio also has an influence on the characteristic stresses and the post-peak deformation behavior, which will be discussed in the next section.

important stress thresholds have been identified based on laboratory compression test results.28,57–60 These stress thresholds are: crack closure stress σcc, crack initiation stress σci, crack damage stress σcd and peak stress σf.28,53 New microcracks initiate at σci and the number of cracks increases rapidly when σcd is reached. Based on laboratory test data, it is found that σci and σcd are approximately 0.3–0.6 and 0.7–1.0 times of the peak strength, respectively.28,53,57–59,61,62 σci and σcd of LdB granite are about 0.4–0.45 and 0.7–0.8 times of peak strength, respectively.53 Based on the stress–strain curves, volumetric strain–crack volumetric strain curves and AE evolution, different approaches are proposed to determine the characteristic stresses.28,53,62–64 A method combining the cumulative number of tensile and shear cracks, volumetric strain and crack volumetric strain is used to estimate the characteristic stresses in this study, details about the method can be found in Wang, Cai9. It should be noted that only the cumulative number of

4.5. Influence of JCp/JσTp ratio on characteristic stresses (σci, σcd, σf) When considering the evolution of rock failure process, several

Fig. 8. Effect of JCp/JσTp ratio on macroscopic tensile strength: (a) C1; (b) C2 and C3. 67

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 9. Effect of the JCp/JσTp ratio on stress threshold σci, σcd, and σf for simulation cases C1, C2, and C3.

tensile and shear cracks are used to estimate the stress thresholds due to the poor volumetric strain obtained in C1. The influence of the JCp /JσTp ratio on the characteristic stresses is presented in Fig. 9. The three stress thresholds are normalized to the corresponding values obtained from the laboratory tests for the LdB granite (i.e., σf, σcd, and σci are 200, 172, and 82 MPa, respectively, see Table 2). It is seen from Fig. 9(a), (c) and (e) that for a fixed grain contact cohesion of JCp = 100 MPa, the three stress thresholds decrease as the JCp /JσTp ratio increases from 2 to 10 for the three simulation cases. Compared the extent of σci variation with these of σcd and σf, it is found that the contact's tensile strength (JσTp) has a large influence on σci and a relatively small impact on σcd and σf for the three simulation cases. In addition, when JσTp is fixed at 21.5 MPa for C1 (Figs. 9(b)) and 15 MPa for C2 (Fig. 9(d)) and C3 (Fig. 9(f)), the characteristic stresses, in general, increase as the JCp/JσTp ratio increases but at different extents. When the contact's tensile strength is fixed, σf and σcd depend largely on

the contact's cohesion (especially for C2 and C3), while σci shows a small dependency on the contact's cohesion. Furthermore, the characteristic stresses of C3 show small variations compared with these of C1 and C2. As can be seen from the results shown in Fig. 9, the JCp/JσTp ratio has a large influence on the damages of the grain contacts for the three simulation cases. Hence, a proper determination of the JCp/JσTp ratio is required to match these important characteristic stresses in 3DEC-GBM modeling. 4.6. Influence of microscopic strength properties on macroscopic C and In numerical modeling, a calibrated 3DEC-GBM model should match not only the macro-deformation properties (i.e., E and υ) but also the characteristic stresses (i.e., σf, σcd, σci, σT) and the macroscopic strength parameters (i.e., C and ϕ). A suite of triaxial compression test simulations are conducted to derive the strength envelopes and 68

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 10. Effect of grain contact cohesion (a, b), peak friction angle (c, d) and residual friction angle (e, f) on normalized macroscopic C and ϕ for rigid grains (C1) and elastic grains (C2).

investigate how the grain and grain contact's microscopic properties affect the macroscopic C and ϕ, which are the Mohr–Coulomb strength parameters for representing strength envelopes. As mentioned above, three cases (C1, C2 and C3) with different grains are considered to reveal the role of grains and grain contacts on the macro-mechanical behavior of specimens.

(10°), macroscopic C increases with the increase of JCp (from 22 to 66 MPa); on the other hand, macroscopic ϕ shows a decreasing trend as JCp increases (small variations compared with that of C, Fig. 10(a)). Furthermore, when the grain contact's JCp and Jϕr are kept at 44 MPa and 10°, both C and ϕ show an increasing trend due to the increase of Jϕp (Fig. 10(c)). In addition, as shown in Fig. 10(e), as Jϕr increases from 5° to 25°, macroscopic C increases while ϕ deceases. For rigid grains, the rock's macroscopic cohesion and friction angle are functions of the grain contact's cohesion, peak and residual friction angles. Very small residual friction angle (i.e., 5°) needs to be used to off-set the effect of block interlocking in order to match the macroscopic C and ϕ.

4.6.1. Rigid grain (C1) For C1, the macro-mechanical behavior is determined by the properties of the grain contacts. Once JKn, JKs and JσTp are determined, the influence of other three contact properties (JCp, Jϕp and Jϕr) on macroscopic C and ϕ are investigated in this section. Fig. 10 presents the variations of macroscopic C and ϕ with different grain contact's properties JCp, Jϕp and Jϕr. Macroscopic C and ϕ are normalized to the corresponding strength properties of LdB granite (Table 2). It is found from Fig. 10(a) that for fixed Jϕp (40°) and Jϕr

4.6.2. Elastic grain (C2) Similar to the C1 case, the influence of contact properties JCp, Jϕp, and Jϕr on macroscopic C and ϕ are investigated for C2. Fig. 10(b), (d) and (f) show the variations of macroscopic C and ϕ with different grain 69

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

contact properties of JCp, Jϕp and Jϕr. For fixed Jϕp (20°) and Jϕr (5°), as shown in Fig. 10(b), C increases quickly with the increases of JCp (from 15 to 75 MPa), while the increase of JCp has a limited impact on ϕ. It can be understood that macroscopic C is strongly influenced by the contact's cohesion. Furthermore, under different Jϕp (from 20° to 40°), different trends of C variations are found when compared the trends of C2 with that of C1 for fixed JCp (40 MPa) and Jϕr (5°). ϕ shows a clear increasing trend while C presents a slightly decreasing trend with the increase of Jϕp (Fig. 10(d)), which means that ϕ is closely related to Jϕp. In addition, as shown in Fig. 10(e), similar results are found for both C1 and C2. C increases while ϕ deceases as Jϕr increases from 5° to 25°. Again, a small residual friction angle needs to be used to calibrate the macroscopic C and ϕ. For C2, grains are deformable but not breakable. The deformation behavior is therefore controlled by both the micro-properties of grains and contacts. Similar to the C1 case, the overall strength in C2 is determined by the contact's strength. Compared with C1, the main advantage of using elastic grains is that volumetric strain can be properly captured, which will be further discussed in Section 6.2.

dominant role in controlling the overall strength of a specimen. In other words, the strength of grains controls overall rock strength when the contact's strength is sufficiently high (i.e., C and ϕ show small variations when Jϕp > 55° in Fig. 11(b)). Fig. 12(a) and (b) show the variation of C and ϕ with different grain properties of Cp and ϕp, respectively. It is seen that grain contact's cohesion influence the macroscopic C and ϕ; both increase as Cp increases, although the rate of increase of ϕ is smaller than that of C. Furthermore, ϕ increases with the increases of ϕp (from 55° to 75°), while C shows a decreasing trend as ϕp increases (Fig. 12(b)). This means that grain's peak friction angle influences both C and ϕ. For deformable and breakable grains (C3), it is seen from Figs. 11 and 12 that C and ϕ are functions of the strength parameters of both the grains and the contacts. The model complexity increases significantly when a proper consideration of contact's residual friction angle is needed. 5. Calibration procedure Based on the parametric study conducted in Section 4 for the three simulation cases, it is found that each grain-scale micro-parameter has a different role in affecting the macro-mechanical properties of 3DECGBMs. Numerical models should be calibrated prior to perform modeling. The baseline of a calibrated model is that the modeling results should match the mechanical properties obtained from laboratory tests. However, the input grain-scale properties cannot be selected randomly; they should be chosen from a constrained “open system.” For instance, the deformation parameters of grains should respect the actual deformation properties of individual minerals. The strengths of grains and grain contacts should be assigned based on reasonable simplifications and the values are in representative ranges, although there are no solid methods to measure these micro-parameters directly. Hence, calibration procedures are developed for different considerations of grains (C1, C2 and C3) based on the parametric study conducted in Section 4. As shown in Fig. 13, the calibration targets include not only the deformation properties but also characteristic stresses and macroscopic strength parameters. It should be pointed out that all the property parameters of grains and grain contacts should be selected carefully and they should be in acceptable ranges. The calibration procedure is listed in the following.

4.6.3. Deformable and breakable grain (C3) Based on the results of C1 and C2, it is found very low residual friction angle (i.e., 5°) of contacts needs to be used, which is too low for residual friction angles of rock. The purpose of considering C3 is to address the issue of very low residual friction angle of contacts used in C1 and C2 simulations. For C3, a strain-softening model is assigned to grains to simulate breakable grains. In this consideration, the overall deformation and strength behaviors are controlled by both the grains and the contacts, which in turn increase the number of micro-parameters and hence the model complexity. Based on reference values in literature,1,7,54–56 some assumptions are made in the following to simplify the parametric study of C3:

• For contacts: a. Peak strength: JCp and Jϕp are the two main micro-properties that need to be investigated; JσTp is assumed to be 14.5 MPa. b. Residual strength: Jϕr is assumed to be 31°, which is equal to the average basic friction angle of granite; JCr and JσTr are assumed to be zero.

•

(A) Deformation property calibration (for C1, C2 and C3)

For grains:

(1) Generate grain structures in Neper using the grain size distribution method and import the block data file into 3DEC and assign initial reference values of micro-properties to the grains and grain contacts. (2) Change contact's normal to shear stiffness ratio (JKn/JKs) to calibrate the Poisson's ratio (υ) using UCS tests. (3) Keep the calibrated JKn/JKs ratio and adjust the contact's normal and shear stiffness to obtain the Young's modulus (E) using UCS tests.

a. Peak strength: Cp and ϕp are the two strength properties of grains that need to be studied; σTp is assumed to be 29 MPa (two times of JσTp); ϕp is assumed to be 26°. b. Residual strength: Cr and ϕr are assumed to be 5 MPa and 31°, respectively. σTr is assumed to be zero; ψr is assumed to be 22°. Considering the brittle failure of LdB granite, the characteristic plastic strain (εps) for cohesion and friction angle of grains are assumed to be 0.1% and 0.15%, respectively.

(B) Strength property calibration

The influences of contact properties JCp and Jϕp on macroscopic C and ϕ are investigated when the grain strength is fixed. Fig. 11(a) and (b) show the variations of macroscopic C and ϕ with different grain contact properties of JCp and Jϕp, respectively. C increases with the increases of JCp (from 21 to 56 MPa), while ϕ shows a decreasing trend as JCp increases (Fig. 11(a)). This means that the contact's cohesion influences C and ϕ. In addition, under different Jϕp (from 35° to 75°), different trends are found for C and ϕ (Fig. 11(b)). As Jϕp increases from 35° to 55°, C increases firstly and then shows a small variation when Jϕp is greater than 55°. In the meantime, Jϕp has a small influence on ϕ, which shows a small variation as Jϕp increases, different than the trends in C1 and C2. As mentioned above, both the strengths of grains and contacts contribute to the macroscopic strength in C3, which means that the lower strength component (either grain or contact) has a

(4) Vary the contact's tensile strength (JσTp) to achieve the target macroscopic tensile strength (σT) using either direct tension or Brazilian tests (for C1, C2 and C3). (5) For C1 and C2, change the contact's residual friction angle (Jϕr) with the calibrated tensile strength (JσTp) of contacts to obtain a range of contact's residual friction angle. For C3, Jϕr is assumed equal to Jϕb. (6) Change the contact cohesion (JCp) with the calibrated JσTp to obtain a range of contact cohesion based on the laboratory characteristic stresses (σci, σcd and σf) using UCS tests (for C1, C2 and C3). (7) For C1 and C2, adjust the contact's friction angle (Jϕp) and cohesion JCp to calibrate the macroscopic cohesion (C) and friction angle (ϕ) 70

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 11. Effect of grain contact's cohesion (a) and contact's peak friction angle (b) on macroscopic C and ϕ for the case of breakable grains (C3).

Fig. 12. Effect of grain cohesion (a) and grain peak friction angle (b) on macroscopic C and ϕ for the case of breakable grains (C3).

using triaxial compression tests. For C3, adjust Jϕp and JCp in reasonable ranges to match the macroscopic C and ϕ. (8) For C3, following Step (7), grain strength parameters (Cp and ϕp) are also adjusted to obtain macroscopic C and ϕ using triaxial compression tests. (9) Steps (7) and (8) are iterated several times to obtain the desired cohesion and friction angle because both the grains and the contacts are allowed to fail.

6.1. UCS test simulation results Fig. 14 shows the mechanical behavior of the calibrated synthetic LdB granite model under uniaxial compression for C1, C2 and C3 using the calibrated grain and contact properties shown in Table 4. A summary of the laboratory test results and the simulation results of the three cases are given in Table 5. Based on the errors shown in Table 5, it is seen that a good agreement is obtained for most of the target properties between the laboratory test and the numerical simulation results. From the axial stress–strain curves shown in Fig. 14, the three cases produce similar deformation responses prior to reaching the crack initiation stress (i.e., 80 MPa). Furthermore, from the lateral stress–strain curves, it is seen that all curves of the three cases are overlapped (i.e., < 50 MPa) and a large overlapped portion is found between the curves of C2 and C3 (i.e., < 140 MPa), which further indicates that the Poisson's ratios are the same for the three cases. It means that the deformation behavior of synthetic LdB can be modeled successfully for all the three cases. Similar UCS values are obtained for C1 (206 MPa), C2 (202 MPa) and C3 (204 MPa) (see Fig. 14 and Table 5) with errors less than 3%. In addition, the two characteristic stresses σci and σcd are also reasonably captured for the three cases, except that the crack initiation stress of C1 is 20.4% smaller than the test value. It should be mentioned again that the characteristic stresses of C1 are estimated based on the accumulated tensile and shear cracks, while the stress levels of C2 and C3 are based on both the accumulated tensile and shear cracks and the crack volumetric strains. For C1, σcd cannot be estimated based on the reversal point in the crack volumetric strain curves, as shown by the red dash line in Fig. 14. In addition, as shown in Fig. 14, large variations in the

6. Simulation results from calibrated models Based on the calibration procedure presented in the previous section (Fig. 13), 3DEC-GBM models are calibrated and the calibrated microparameters of the grains and the contacts of the three cases (C1, C2, and C3) are summarized in Table 4. As pointed out by Potyondy and Cundall11 one-mineral idealization of rock is simpler than the multi-mineral idealization and minor differences in the deformation behaviors of the models are exhibited although the grain structure (shape and size distribution) is the sole heterogeneity source. It should be mentioned that the calibrated model can be viewed as a simplified grain-based model due to the fact that only one contact properties are used in the calibration. Otherwise, ten contact types are needed to calibrate a model if four minerals are considered in the model, which makes the calibration procedure more complex and the uncertainties of the contact properties associated with the grains more pronounced. The purpose of this study is to explore different grain considerations on numerical modeling of rock failure instead of analyzing the effect of contact strength heterogeneities on macroscopic mechanical behavior; therefore, no attempt is made to differentiate the contact types. 71

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 13. Calibration procedures of grain and grain contact parameters for 3DEC-GBMs.

stress–strain curves are observed when the axial stress surpasses the crack initiation (i.e., > 100 MPa), especially for the post-peak axial stress–strain curves. It is seen that C3 shows a more brittle behavior than these of C1 and C2 in the post-peak deformation stage, which indicates that the consideration of breakable grain is necessary in modeling deformation response of brittle rocks. Furthermore, the strengths of C1 and C2 are determined purely by the strengths of contacts and very low contact's residual friction angle is needed to offset the strong interlocking effect of grains. Insufficient brittleness is still

observed in Fig. 14, even though very small contact's residual friction angles (i.e., 5°) are used in the modeling, which in turn indicates that an adjustment of the input parameter of contact is required to match the laboratory test results. Fig. 14 also shows the cumulative numbers of tensile and shear cracks for the three modeling cases. For all the three cases, it is found that the cumulative number of tensile crack is higher than that of shear crack and tensile cracks appear earlier than the shear cracks (Fig. 14), which means that the dominant failure mechanism in the compression

Table 4 Calibrated grain and contact property parameters of LdB granite for three study cases. Grain properties Deformation parameters (C1 C2 C3) Strength parameters (C3)

Contact properties ρ (kg/m3) E (GPa) υ Peak grain strength

Grain strength when strain is greater than the characteristic plastic strain Plastic strain

2605 78.13 0.22 Cp (MPa)

62

ϕp (°)

70

σTp (MPa)

29

ψp (°) Cr (MPa) ϕr (°) σTr (MPa) ψr (°) εcps (%) εfps (%)

26 5 31 0 22 0.1 0.1

72

JKn (GPa/m) JKs (GPa/m) JKn/JKs Peak contact strength

45000 (C1), 250000 (C2), 380000 (C3) 9000 (C1), 99210 (C2), 156000 (C3) 5 (C1), 2.52 (C2), 2.52 (C3) JCp (MPa) 53 (C1), 51.5 (C2), 28 (C3) Jϕp (°) 37 (C1), 31.5 (C2), 56 (C3) JσTp (MPa) 21.5 (C1), 14.5 (C2 and C3)

Residual contact strength

JCr (MPa) Jϕr (°) JσTr (MPa)

0 (C1, C2, C3) 5 (C1 and C2), 31 (C3) 0 (C1, C2, C3)

Note: contact's residual cohesion and tensile strength are assumed to be zero for the three cases

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

Fig. 15. Macroscopic fracture patterns at 0.5% axial strain for all three grain cases with different amplifications (× 2 and × 5).

volumetric strain, and the crack volumetric strain are calculated for C1, C2 and C3. In the volumetric strain curves, positive and negative strains indicate volumetric contraction and dilation, respectively. The total volumetric strain has two contribution sources: the deformations of the grain itself and the contact. Not surprisingly, the rigid grain case C1 shows very different volumetric behavior when compared with these of C2 and C3. Very small volumetric strain is observed and volumetric dilation occurs easily. This is because that the volumetric strain is solely contributed from the sliding, opening and rotation of the contacts, with no contribution from the rigid grains. On the contrary, proper volumetric strain curves are obtained for both C2 and C3 due to the consideration of deformable grains, which have a large contribution to the total volumetric strain. Because poor volumetric strain is captured in C1, the characteristic stresses cannot be estimated from the crack volumetric strain. That is why only the cumulative numbers of tensile and shear cracks are used to derive the characteristic stresses in C1 (rigid grain case). Macroscopic fracture patterns at 0.5% axial strain for the three simulation cases are presented in Fig. 15. Several vertical splitting fractures are formed near the left and the right model boundaries and many minor sub-vertical splitting fractures (parallel or sub-parallel to the loading direction) are also developed within the models. As mentioned in Section 6.1, tensile cracking is the dominant failure mechanism in the uniaxial compression modeling. From the modeling results, it is observed that the main orientation of tensile cracks is perpendicular or sub-perpendicular to the loading direction, which would contribute to the final formation of vertical splitting fractures. For C1 and C2, the overall strengths are controlled by the contacts and hence relatively large contact deformations are observed. Relatively small contact deformations are observed for C3 when compared with C1 and C2 because the overall strength of C3 is governed by the strengths of the grains and the contacts. For C3, the grains tend to fail after the stress passes σcd and

Fig. 14. Modeling curves of stress–strain, associated total volumetric and calculated crack strain–axial strain and cumulative number of tensile and shear cracks for the calibrated UCS tests for C1, C2, and C3.* : Equations are from Martin53.

modeling is tensile cracking. Furthermore, it is observed that cases C3 and C1 produce the highest and the smallest numbers of tensile and shear cracks, respectively. This further indicates that contacts fail at different times during numerical loading when different types of grains are used in the modeling. 6.2. Volumetric strain and final fracture patterns As shown in Fig. 14, the elastic volumetric strain, the total

Table 5 Property comparison between laboratory test results and numerical modeling results. Properties Deformation Stress levels Strength

E (GPa) υ σci (MPa) σcd (MPa) C (MPa) ϕ (°) σT (MPa) UCS (MPa)

Lab. Re.

C1

Er. (%)

C2

Er. (%)

C3

Er. (%)

69.0 0.26 82 172.0 30.9 55.6 9.3 200.0

71.2 0.264 65.3 * 176.2 * 32.2 56.4 10.43 206.2

+ 3.2 + 1.5 − 20.4 + 2.4 + 4.2 + 1.4 + 12.2 + 3.1

68.6 0.254 84.2 178.4 31.1 55.9 9.44 201.7

− 0.6 − 2.3 + 2.7 + 3.7 + 0.6 + 0.5 + 1.5 + 0.9

69.6 0.258 87.6 186.5 31.6 55.2 9.36 204.3

+ 0.7 − 0.8 + 6.8 + 8.4 + 2.3 − 0.7 + 0.6 + 2.2

Notes: Lab. Re. denotes the Laboratory Results; Er. denotes Error; * means the stress levels of C1 are estimated only based on the accumulate number of tensile and shear cracks. Stress levels are estimated based on both the accumulate number of tensile and shear cracks and the crack volumetric strain for C2 and C3. 73

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

of specimens. 6.3. Direct tension test simulation results Direct tension tests are simulated for the three cases with a vertical loading velocity of 0.01 m/s, which is reduced ten times compared with that used in the UCS test simulations. Fig. 16 presents tensile stress–axial strain curves and the final fracture patterns for the three simulation cases. The direct tensile strengths are 10.43, 9.44, and 9.36 MPa for C1, C2, and C3, respectively (see Fig. 16 and Table 5). All the three cases show a brittle tensile splitting behavior with a rapid post-peak stress drop. Almost all the cracks are tensile cracks and there are only a very small number of shear cracks (not included in this figure), which indicates that tensile cracking is the dominant failure mechanism in direct tension. As expected, the tensile strength of the grain contact determines the model's macroscopic tensile strength, which had been discussed in Section 4.4. The macroscopic failure mode, shown as inserts in Fig. 16 (amplified 40 times), is tensile splitting (in terms of pure contact opening mode) along the major fracture and fracture branches developed within the specimen. For the three simulation cases, the final macroscopic fracture is not located in the model's center height, which may be caused by the non-uniform distribution of grain size. Furthermore, the macroscopic tensile fracture, although located in the same part of the specimen (i.e., lower part of the specimen), shows slightly different failure paths for the three simulation cases.

Fig. 16. Tensile stress–axial strain curves and fracture patterns for the three cases in direct tension tests.

in the post-peak deformation stage because the deformation of the grains is governed by the strain-softening model. Hence, the overall failure pattern in C3 is a combination of grain and contact failures. Although axial splitting is the dominant macroscopic fracture mode for the three simulation cases, the final failure modes of the three cases show a combination of tensile and shear fractures due to the coalescence of factures in the post-peak deformation stage, especially for C3. This suggests that the consideration of grain's deformation model has a large influence on the cracking patterns and thus the final failure mode

Fig. 17. Stress–strain curves under different confining pressure for C1 (a), C2 (b), and C3 (c); (d) summarizes the triaxial compression test modeling results for the three simulation cases. 74

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

6.4. Triaxial compression test simulation results

Table 6 Advantages and disadvantages of three types of grains in UDEC and 3DEC-GBM.

Fig. 17 presents the deviatoric stress–axial strain curves of the triaxial compression tests for the three simulation cases. The synthetic LdB granite specimens are subjected to confining pressures ranging from 0 to 20 MPa. It is seen that the strength increases as the confining pressure increases for all cases. Furthermore, different axial strains are required to obtain the peak confined strength. C1 and C2 require relatively large axial strains (i.e., > 1%) to reach the peak confined strength compared with C3 (i.e., < 0.6%). This is caused by the different grain considerations. The strengths of C1 and C2 are solely determined by the contact's strength and the crack nucleation and propagation are suppressed by the applied confinement. The overall strength for the breakable grain case C3 is determined by both the strengths of the contacts and the grains, and a shear band is formed due to grain failure under high confinement. Although the calibrated Mohr–Coulomb strength parameters (i.e., C and ϕ) are in good agreement with the laboratory test results for the three simulation cases (Table 5), the post-peak deformation behaviors are different (Fig. 17(a), (b) and (c)). Fig. 17(d) presents a summary of the triaxial stress–strain curves for C1, C2 and C3. It is seen that the post-peak deformation behaviors are different for the three simulation cases although the confined strengths are similar. C3 shows a brittle behavior and C1 and C2 show less brittle behaviors. This indicates that different types of grains affect the mechanical behavior of synthetic rocks, especially at the post-peak deformation stages under confinement.

Comparison item

C1 (Rigid)

C2 (Elastic)

C3 (Breakable)

Numbers of input parameters Computation time Calibration procedure Input contact friction angle Deformation parameters (E, υ) Macroscopic strength parameters (UCS, σt, C, and ϕ) Volumetric strain Characteristic stresses Capturing brittleness of axial stress–strain curve Post-peak behavior under high confinement

Few Short Simple Very low Reasonable

Few Long Simple Very low Reasonable

Many Very long Very complex Reasonable Reasonable

Reasonable

Reasonable

Reasonable

Poor Poor Poor

Good Reasonable Poor

Good Reasonable Reasonable

Strainhardening

Strainhardening

Strain-softening

Notes: C1 is the rigid case without zoning; C2 is the deformable grain case by assigning the elastic model; C3 is the breakable grain case by assigning a yield criterion or with sub-grain considerations; the comparison are relative among the three cases of grain considerations.

of the stress–strain curve. Based on the advantages and disadvantages discussed above, it is suggested that the grain type should be selected based on the modeling purpose and computation capacity. 7. Conclusions

6.5. Discussion on selection of grain type

Grain and grain contact, which control the micro-mechanical behaviors and consequently the macro-mechanical behavior of a synthetic rock specimen, are the two key components in UDEC- and 3DEC-GBM simulations. In this study, grain structures are generated in Neper using a normal distribution of grain size with cutoffs and then imported into 3DEC for deformation behavior modeling. A comprehensive parametric study is conducted based on the laboratory test results of Lac du Bonnet granite considering three types of grains: rigid, elastic, and breakable grains. Based on the results of the parametric study, calibration procedures are recommended to determine the micro-mechanical parameters of contacts and grains for the three types of grains. The influence of grain type on the mechanical response of the models is investigated using 3DEC-GBM numerical models by conducting uniaxial compression, direct tension and triaxial compression tests. The advantages and disadvantages of using the three types of grains are also discussed. The modeling results show that the micro-properties of the contacts and the grains influence the macroscopic deformation and strength properties of the synthetic rock specimens differently. Similar calibration procedures are recommended for the cases of rigid and elastic grains, while a more complex calibration procedure is suggested for the case of breakable grains. It is found that all the three simulation cases give good results for the pre-peak deformation (E and υ) and the macroscopic strength (UCS, σt, C and ϕ) parameters, but there are large differences in terms of post-peak deformation behavior and the realism of input parameters. In addition, it is found that the use of the three types of grains has advantages and disadvantages in modeling. For rigid grains, the computation time is less but the volumetric strain cannot be modeled accurately and the contact's residual friction angle has to be near zero in order to match the macroscopic post-peak deformation behavior due to the interlocking effect. For elastic grains, the main issues are that very low contact's residual friction angle has to be used and it is difficult to model the brittle behavior (i.e., sudden loose of the strength capacity post-peak). In addition, the complex calibration procedure and excessive computation time bring challenges to modeling using the breakable grains but it resolves the issues of inaccurate volumetric strain, low contact residual friction angle and strain-hardening in the simulation results. It should be mentioned that there are limitations in the modeling

In this study, a comprehensive parametric study is conducted considering three types of grains using the 3DEC-GBM modeling approach. For C1, C2 and C3, grain and grain contact's strength parameters have a large influence on the macroscopic behavior of the model. In general, micro-property selection is in an “open system”, which means that various combinations of micro-properties of grains and grain contacts can produce a “calibrated” model, which can represent the majority of the macroscopic rock properties in terms of deformability, strength, and failure pattern. However, all the grain-scale properties should be chosen in a “constrained open system”, and the selected properties should respect actual grain conditions and are in a reasonable range. There are advantages and disadvantages for each type of grain considering the number of model input parameters, computation time, model calibration complexity and modeling results. Based on the results of the comprehensive parametric study presented in Section 4 and other studies using UDEC-GBM and 3DEC-GBM from some papers,i.e.,3,7,38,43,45 Table 6 presents a brief summary of the advantages and disadvantages of the three types of grains in UDEC- and 3DEC-GBM numerical modeling. It is seen from Table 6 that all the three types of grains give good modeling results for the pre-peak deformation behavior (E and υ) and the macroscopic strength parameters (UCS, σt, C and ϕ) after calibrating micro-parameters of contacts and grains. However, each grain type has limitations in terms of the realism of the input parameters, the number of input grain-scale parameters, the ability to capture brittle rock failure, and computation time, etc. For instance, volumetric stain cannot be captured properly for C1 (rigid grain) and a very low contact's residual friction angle is needed to off-set the block interlocking effect for C1 (rigid grain) and C2 (elastic grain). Strain-hardening behavior under high confinement is found for C1 and C2. For the breakable grain case C3, because there are many grain-scale input parameters, the calibration is complex and time-consuming, and the demand on computation power is high. In many studies, rigid or elastic grains are normally used to conduct GBM modeling to simplify the calibration procedure and reduce computation time but one has to sacrifice the realism of the contact's residual friction angle and brittleness 75

International Journal of Rock Mechanics and Mining Sciences 115 (2019) 60–76

X. Wang, M. Cai

due to time-consuming model calibration and the high demand on computational power; hence single mineral and contact properties are considered and equivalent 2D plane strain simulation, instead of full 3D modeling, is conducted. With these limitations in mind, the results of this study can be used to guide the selection of grain type and the calibration process in GBM modeling. When the parallelization of 3DEC becomes available, the work can be extended to full 3D GBM modeling while considering the heterogeneities of grains and contacts.

29. 30. 31. 32. 33.

Acknowledgements

34.

The authors would like to acknowledge financial support from The Natural Science and Engineering Research Council of Canada (NSERC) (CRDPJ 470490-14), The Nuclear Waste Management Organization of Canada (NWMO), and Rio Tinto. The authors also would like to thank Dr. Romain Quey for providing the open-source software Neper (http:// neper.source-forge.net/) and help in using the software.

35. 36. 37. 38.

References

39.

1. Gao F, Stead D, Elmo D. Numerical simulation of microstructure of brittle rock using a grain-breakable distinct element grain-based model. Comput Geotech. 2016;78:203–217. 2. Wong TF, Wong RHC, Chau KT, Tang CA. Microcrack statistics, Weibull distribution and micromechancial modeling of compressive failure in rock. Mech Mater. 2006;38:664–681. 3. Lan H, Martin CD, Hu B. Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. J Geophys Res: Solid Earth. 2010;115(B1). 4. Potyondy D. A grain-based model for rock: approaching the true microstructure. Proc Rock Mech Nord Ctries. 2010. 5. Bahrani N, Kaiser PK, Valley B. Distinct element method simulation of an analogue for a highly interlocked, non-persistently jointed rockmass. Int J Rock Mech Min Sci. 2014;71:117–130. 6. Bewick R, Kaiser P, Bawden W, Bahrani N. DEM simulation of direct shear: 1. Rupture Constant Norm Stress Bound Cond Rock Mech Rock Eng. 2014;47(5):1647–1671. 7. Farahmand K, Diederichs MS. A calibrated synthetic rock mass (SRM) model for simulating crack growth in granitic rock considering grain scale heterogeneity of polycrystalline rock. In: Proceedings of the 49th US Rock Mechancis / Geomechanics Symposium. San Francisco, CA, USA. Paper 15-430. 2015:14. 8. Peng J, Wong LNY, Teh CI. Influence of grain size heterogeneity on strength and microcracking behavior of crystalline rocks. J Geophys Res: Solid Earth. 2017;122(2):1054–1073. 9. Wang X, Cai M. Modeling of brittle rock failure considering inter-and intra-grain contact failures. Comput Geotech. 2018;101:224–244. 10. Lim SS, Martin CD, Åkesson U. In-situ stress and microcracking in granite cores with depth. Eng Geol. 2012;147:1–13. 11. Potyondy D, Cundall P. A bonded-particle model for rock. Int J Rock Mech Min Sci. 2004;41(8):1329–1364. 12. PLAXIS 2D&3D. Accessed on July 4, 〈https://www.plaxis.com/product〉; 2018. 13. Rocscience. Accessed on July 4, 〈https://www.rocscience.com/rocscience/ products〉; 2018. 14. ABAQUS. Accessed on July 4, 〈https://www.3ds.com/productsservices/simulia/ products/abaqus/〉; 2018. 15. Itasca Consulting Group Inc., 3DEC and UDEC. Accessed on July 4, 〈https://www. itas-cacg.com/software〉; 2018. 16. Map3D. Accessed on July 4, 〈https://www.map3d.com/〉; 2018. 17. Cundall PA. A discontinuous future for numerical modelling in geomechanics? Proc Inst Civil Eng-Geotech Eng. 2001;149(1):41–47. 18. Šmilauer V, et al. Yade Documentation. 2nd ed. The Yade Project; 2015〈http://yadedem.org/doc/〉. 19. Rockfiled. Elfen (Finite Element-Discrete Element). Accessed on July 4; 2018. 〈www. rockfi-eldglobal.com〉. 20. Xiang J, Munjiza A, Latham J-P, Guises R. On the validation of DEM and FEM/DEM models in 2D and 3D. Eng Comput. 2009;26(6):673–687. 21. Mahabadi O, Grasselli G, Munjiza A. Y-GUI: a graphical user interface and pre-processor for the combined finite-discrete element code, Y2D, incorporating material heterogeneity. Comput Geosci. 2010;36(2):241–252. 22. Munjiza AA. The combined finite-discrete element method. John Wiley & Sons; 2004. 23. 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(1):187–206. 24. Diederichs MS. Instability of hard rockmasses, the role of tensile damage and relaxation. Ph D Univ Waterloo. 2000:566. 25. Hazzard JF, Young RP. Simulating acoustic emissions in bonded-particle models of rock. Int J Rock Mech Min Sci. 2000;37(5):867–872. 26. Tang C, Liu H, Lee P, Tsui Y, Tham L. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—part I: effect of heterogeneity. Int J Rock Mech Min Sci. 2000;37(4):555–569. 27. Cai M, Kaiser PK. Numerical simulation of the Brazilian test and the tensile strength of anisotropic rocks and rocks with pre-existing cracks. Int J Rock Mech Min Sci. 2004;41(Supplement 1):478–483. 28. Cai M, Kaiser PK, Tasaka Y, Maejima T, Morioka H, Minami M. Generalized crack

40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64.

76

initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int J Rock Mech Min Sci. 2004;41(5):833–847. Christianson M, Board M, Rigby D. UDEC simulation of triaxial testing of lithophysal tuff. Golden Rocks 2006. In: Proceedings of the 41st US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association; 2006. Na Cho, Martin C, Sego D. A clumped particle model for rock. Int J Rock Mech Min Sci. 2007;44(7):997–1010. 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(2):219–244. Potyondy D. A flat-jointed bonded-particle material for hard rock. 46th US Rock mechanics/geomechanics symposium. Am Rock Mech Assoc. 2012. Rong G, Liu G, Hou D, Zhou C-b. Effect of particle shape on mechanical behaviors of rocks: a numerical study using clumped particle model. Sci World J. 2013;2013. Azocar KD. Investigating the mesh dependency and upscaling of 3D grain-based models for the simulation of brittle fracture processes in low-porosity crystalline rock. Geological Sciences and Geological Engineering. Canada: MA.S. Queen's University; 2016:150. Guo S, Qi S, Zou Y, Zheng B. Numerical Studies on the Failure Process of Heterogeneous Brittle Rocks or Rock-Like Materials under Uniaxial Compression. Materials. 2017;10(4):378. Liu G, Cai M, Huang M. Mechanical properties of brittle rock governed by microgeometric heterogeneity. Comput Geotech. 2018. Kazerani T, Zhao J. Micromechanical parameters in bonded particle method for modelling of brittle material failure. Int J Numer Anal Methods Geomech. 2010;34(18):1877–1895. Ghazvinian E, Diederichs M, Quey R. 3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing. J Rock Mech Geotech Eng. 2014;6(6):506–521. Nygårds M, Gudmundson P. Three-dimensional periodic Voronoi grain models and micromechanical FE-simulations of a two-phase steel. Comput Mater Sci. 2002;24(4):513–519. Zhang K, Wu M, Feng R. Simulation of microplasticity-induced deformation in uniaxially strained ceramics by 3-D Voronoi polycrystal modeling. Int J Plast. 2005;21(4):801–834. Li H, Li K, Subhash G, Kecskes LJ, Dowding RJ. Micromechanical modeling of tungsten-based bulk metallic glass matrix composites. Mater Sci Eng: A. 2006;429(1):115–123. Damjanac B, Board M, Lin M, Kicker D, Leem J. Mechanical degradation of emplacement drifts at Yucca mountain—A modeling case study: Part II: Lithophysal rock. Int J Rock Mech Min Sci. 2007;44(3):368–399. Park JW, Park C, Song JW, Park ES, Song JJ. Polygonal grain‐based distinct element modeling for mechanical behavior of brittle rock. Int J Numer Anal Methods Geomech. 2017;41(6):880–898. Neper. Accessed on 4 July; 2018. 〈http://neper.sourceforge.net〉. Gao F, Stead D. The application of a modified Voronoi logic to brittle fracture modelling at the laboratory and field scale. Int J Rock Mech Min Sci. 2014;68:1–14. Cundall PA. A computer model for simulating progressive, large scale movement in blocky rock systems. Symp ISRM, Nancy, Fr, Proc. 1971;2:129–136. Potyondy DO. The bonded-particle model as a tool for rock mechanics research and application: current trends and future directions. Geosystem Eng. 2015;18(1):1–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(2):352–369. Quey R, Dawson P, Barbe F. Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing. Comput Methods Appl Mech Eng. 2011;200(17):1729–1745. Barbe F, Quey R. A numerical modelling of 3D polycrystal-to-polycrystal diffusive phase transformations involving crystal plasticity. Int J Plast. 2011;27(6):823–840. Quey R, Renversade L. Optimal polyhedral description of 3D polycrystals: method and application to statistical and synchrotron X-ray diffraction data. Comput Methods Appl Mech Eng. 2018;330:308–333. Bass DJ. Elasticity of minerals, glasses and melts. Mineral Phys Crystallogr. 1995;2:45–63. Martin CD. The strength of Massive Lac du Bonnet Granite around underground openings. Civil and Geological Engineering. Ph.D. University of Manitoba, Winnipeg, Canada; 1993:278. Barton NR. The shear strength of rock and rock joints. Int J Rock Mech Min Sci Geomech Abstr. 1976;13(10):1–24. Barton NR, Choubey V. The shear strength of rock joints in theory and practice. Rock Mech. 1977;10(1–2):1–54. Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci. 1997;34(8):1165–1186. Bieniawski Z. Mechanism of brittle fracture of rock: part II—experimental studies. Int J Rock Mech Min Sci & Geomech Abstr. 4. Elsevier; 1967 [407IN13409IN1541913408IN14418IN18423]. Bieniawski ZT. Mechanism of brittle fracture of rock: part I—theory of the fracture process. Int J Rock Mech Min Sci & Geomech Abstr. 4. Elsevier; 1967 [395IN1140511404IN12406]. Martin C, Chandler N. The progressive fracture of Lac du Bonnet granite. Int J Rock Mech Min Sci & Geomech Abstr. 31. Elsevier; 1994:643–659. Eberhardt E, Stead D, Stimpson B. Quantifying progressive pre-peak brittle fracture damage in rock during uniaxial compression. Int J Rock Mech Min Sci. 1999;36(3):361–380. Brace W, Paulding B, Scholz C. Dilatancy in the fracture of crystalline rocks. J Geophys Res. 1966;71(16):3939–3953. Eberhardt E, Stead D, Stimpson B, Read R. Identifying crack initiation and propagation thresholds in brittle rock. Can Geotech J. 1998;35(2):222–233. Nicksiar M, Martin C. Evaluation of methods for determining crack initiation in compression tests on low-porosity rocks. Rock Mech Rock Eng. 2012;45(4):607–617. Zhao X, Cai M, Wang J, Li P, Ma L. Objective determination of crack initiation stress of brittle rocks under compression using AE measurement. Rock Mech Rock Eng. 2015;48(6):2473–2484.