The role of transgranular capability in grain-based modelling of crystalline rocks

Computers and Geotechnics 110 (2019) 161–183

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

The role of transgranular capability in grain-based modelling of crystalline rocks

T

X.F. Lia,b,c, H.B. Lia,b, , J. Zhaoc ⁎

a

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China University of Chinese Academy of Sciences, Beijing 10049, China c Department of Civil Engineering, Monash University, Clayton, VIC 3800, Australia b

ARTICLE INFO

ABSTRACT

Keywords: Grain crushing Intergranular and transgranular cracks Grain-based model Discrete element method (DEM) simulation Heterogeneity

Considering that many numerical methods have been developed to study the mechanical properties of rocks based on micro heterogeneities such as mineral size, morphology, composition and boundary defects, an overview of grain-based modelling is presented in this study to compare the advantages and limitations of different grain-mimicking methods. Intergranular and transgranular fracturing models and difficulties in parameter calibrations are discussed. Four typical grain models, UDEC-GBM, GB-FDEM, cluster and clump, are used to represent the two main families: block-based grains and particle-based grains with and without considering grain crushing. Subsequently, how grain crushing affects crack stresses, Hoek-Brown strength parameters mi, localized shearing and cracking transformation is simulated by these grain models. The simulated results indicate that the approaches capable of grain breakage lead to more consistent results in comparison with laboratory tests, for example, evident dilatancy in uniaxial compressive loading, high strength ratios, nonlinear failure strength envelopes and shear bands under high confinement. The role of transgranular capability is significantly important in the simulation of rock deformation using grain models.

1. Introduction Rock is a kind of microscopically heterogeneous material with widely distributed defects, including pores, microcracks, grain boundaries and cleavages. Microcracks may be the main source of heterogeneity in rocks of low porosity, for example, igneous rocks with strong initial cohesive behaviours and calcite-cemented sandstone [82]. A definition of different microcracks is clarified by Kranz [46], and many researchers distinguish micro-fissures as intergranular cracks generated on grain boundaries and transgranular cracks crossing through mineral grains [96,50]. Microcracks are commonly accepted to govern physical properties such as failure strength, elastic wave velocity, fracture toughness and permeability [82,3,50]. From the viewpoint of microscopic fracturing, there are six mechanisms present in granitic rocks under compression stress as shown in Fig. 1. Sliding-induced intergranular cracks are widely observed at low stress levels, especially when the stress exceeds the crack initiation stress. With a comparison of numerical modelling and laboratory tests, Li et al. [50] suggested using the emergence of intergranular cracks as the indicator of crack initiation. When the stress is approximately equal to the crack damage stress

or short-term failure strength, microcracks begin to propagate and coalesce and are accompanied by a proliferation of intragrain cracks, such as short infiltrating tensile cracks within a single mineral, long propagating cleavage that crosses through grains and sliding-induced shear cracks (as shown in Fig. 1b–d). In triaxial compression testing on a compact rock, intensive acoustic emission events are found to cluster along the elongate volume that forms a through-going shear band across the sample [65]. The shear localization is pronounced at higher pressure, and the underlying mechanism is attributed to micro-bulking of mineral grains resulting from parallel or clustered transgranular cracks, as shown in Fig. 1e and f [76,50]. Apart from compression testing, microscopic observation of granite in tension test has revealed that more than 90% of tensile cracks are present as transgranular fractures and that the remaining 10% follow the intergranular fractures [27]. Li et al. [57–58] also found that the proportion of transgranular cracks increases from 47% to 69% when the loading condition is changed from a quasi-static to a high strain rate state (∼335/s). Many laboratory observations have indicated that the role of transgranular fracturing behaviours is significant in the fracturing process of brittle rocks, and at the same time the transgranular fracturing dominates the

⁎ Corresponding author at: State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China. E-mail addresses: [email protected] (X.F. Li), [email protected] (H.B. Li), [email protected] (J. Zhao).

https://doi.org/10.1016/j.compgeo.2019.02.018 Received 9 November 2018; Received in revised form 14 January 2019; Accepted 14 February 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 1. Idealized crack types revealed in granites under vertical compressive stress. Type a is the intergranular (IG):tensile crack on grain boundaries as a result of closed crack sliding; type b is the transgranular (TG):tensile crack induced by the local concentrated tensile stress at the tip of the sliding crack; type c is the long propagating TG:tensile/shear crack across grains; type d is the shear crack across a single grain; type e is the parallel opening cracks induced by grain shearing; and type f is the shear band-induced large-scale grain pulverization (modified after [96,37,82,50]).

macro properties of materials in a wide range from homogeneous rocks to crystalline rocks [96,101,76,50–51,58]. Numerical methods provide a means to quantitatively investigate fracturing behaviours at the grain scale under conditions that are difficult to achieve in the laboratory. The widely used continuous method is the finite element method (FEM) [95,13], followed by the discrete element method (DEM) [23,86,53] and a combination of the two, such as FDEM [78,61] and FEM/DEM [32]. DEMs have been applied to reproduce the fracturing of rocks in the context of quasi-static loading conditions, including the effects of initial pores on the fracture mechanism [7] and micro heterogeneity on rock strength [59], damage processes and shearing banding [24] and the role of micropolar flow on the brittle-ductile transition [109]. Further applications of DEMs have migrated to impact engineering areas, including dynamic fragmentation and grain pulverization [58], softening in granular materials under extreme dynamic loads [48], dynamic fragmentation distribution of rockfalls [108] and fractal crushing of granular materials [38]. DEMs have the advantages of characterizing the separation, overlapping and sliding of elements by simulating fracturing with the debonding of cohesive contacts. Considering the micro heterogeneity present in granitic rocks, attempts to reproduce grain structures have been explored. These studies include continuum-based frameworks, such as the modified finite difference method (FDM) [98], fracturing-capable FEM [106], smooth particle method (SPH) [67] and distinct schemes, for example, the UDEC-grain boundary model [47,28], 3DEC-Voronoi model [99], three-dimensional polygon-based DEM [55], PFC-Grain model [10–11], Irazu-GBM [1] and multiscale DEM [50–51,58]. From the viewpoint of experimental results, the grain-based model has distinguishing features summarized as: (a) heterogeneity in the mechanical properties of mineral grains can be reproduced; (b) the geometric heterogeneity induced by grain boundaries can be revealed; and (c) grains can be crushed under certain stress levels. The term ‘microcrack’ is classified as intragranular, transgranular and intergranular cracks depending on the fracturing mechanism and the location where the cracks form [46]. Intergranular cracks are abbreviated as ‘IG cracks’ to refer to fracturing behaviours occurred on

grain boundaries. The terminologies of intragranular and transgranular cracks are different in the area of geology because the former indicates cracks crossing one grain and the latter indicates cracks crossing through several grains. In this study, we focus on the definition of grain fracturing, therefore, these grain-crossing cracks are collectively called transgranular fracturing. This contribution specifically discusses the importance of transgranular cracking in grain-based modelling and utilizes four different modelling schemes to show how grain crushing influences the crack stresses, Hoek-Brown strength parameter, localized shearing and cracking transition. In section 2, this paper contains an overview of the grain-based methods for numerical modelling of geomaterials and then describes the algorithm of grain-breakable/unbreakable modelling. Subsequently, the calibration procedure of grain-based models is discussed, and the simulation of rock fracturing under different stress states is performed. Finally, the results are presented, including the effect of intergranular cracking on the fracturing process, damage localization and fracturing mechanism. 2. Overview of grain-based modelling algorithm Considering the influence of mineral structure on rock deformation, the grain modelling method is widely developed based on either the continuum framework or the discontinuum framework. From the viewpoint of grain shape, these methods can be categorized as follows: (a) square-based grains, which are commonly implemented in FEM and FDM such as FLAC2D,3D and rock failure process analysis (RFPA) [49,63,18,98,106], (b) triangle-based grains, for example, the FDEM model implemented in Y-Geo [70]; (c) particle-based grains using the particle DEM, a method that is widely employed for geomaterials such as rocks and soil by using the commercial code, PFC2D/3D, or open source codes, YADE (https://yade-dem.org/doc/), LIGGGHTS, and ESyS-Particle; (d) node-based grains implemented in SPH [67], and (e) polygon-based grains in combination with UDEC and 3DEC (as listed in DEC (Table 1)). Square-based grain methodology: Early studies on rock

162

Square-based method, stochastic grain distribution Square-based method, random grain shape Square-based method, imagebased grain modelling

FDM; FLAC2D

FEM; RFPA2D

FDM; FLAC2D

Uniaxial compression test; Brazilian disc test, Wedge chipping Uniaxial compression test

Brazilian disc test

163

FDEM, Irazu-grain

FDEM, GB-FDEM

DEM. 3DEC-trigon model

DEM, UDEC

NMM DEM; 3DEC-Voronoi

Uniaxial compression test

Compression test, Brazilian disc test

Uniaxial compression test, Brazilian disc test

Uniaxial compression test

Uniaxial compression test Uniaxial compression test

Voronoi-shaped and Trigonshaped grains Voronoi tessellation Voronoi tessellation

Trigon-shaped grains

Actual grain shape

Dual-layer model, Voronoi tessellation of mineral grains, triangle-based element intragrain Polygon-based grain

DEM; UDEC-grain model

DEM; UDEC-grain model

Voronoi tessellation

DEM; UDEC-grain model

Uniaxial compression test, Brazilian disc test; fracturing test Brazilian disc test

Uniaxial compression test, Brazilian disc test

Voronoi tessellation

Meshfree; SPH

Dynamic loading

Actual 3D CT method Node-based method

Triangle-based method, actual grain shape

FDEM; X-ray MicroCT

Brazilian disc test

Triangle-based method, actual grain shape

Polygon-based grain

FDEM; X-ray MicroCT

Discrete element approach, irregular lattice model with rigid-body-spring network

Rate dependence of concrete

Grain methodology

Dynamic behaviours of RC structure Mixed mode fracture test and plate impact Actual material heterogeneity

Method

Application

Uncertainty induced by grain boundary modes Fabric-guided micro fracturing

Elastic minerals, material is homogenous

Microfracturing is achieved in 3D DEM modelling

Inter/intra behaviours on crack fracturing of rocks

The role of grain crushing on rock static fracturing

The role of grain heterogeneity on rock tensile strength

Dynamic failure and strain rate mechanism Stress heterogeneity induces brittle extensile fracturing of rocks

Thermal treatment on grain crushing Rock heterogeneity on rock tensile behaviours

Failure process of granite

Micro heterogeneity and microstructures on the tensile behaviours Heterogeneity and mineral properties on UCS and BTS

Actual micro heterogeneity is reproduced and micro parameters are calibrated by micro indentation

Load-carrying capacity of RC concrete

Rate dependence is a behaviour of viscosity of concrete

Results

Elastic minerals

Elastic minerals

Grain crushable

Grain crushable

Grain crushable

Element plasticity

Elastic minerals, grain crushing is unavailable

Damage degradation

Element plasticity

Deformable grains, fracturing is represented by plasticity Element degradation in grains

Deformable grains

Deformable grains with intergranular capability

Rigid grains

Grain crushing ability

(a) 3D grain model (b) No grain crushing

(a) Digital image-based modelling (b) Intergranular cracking is ignorable (a) Three-dimensional model (a) Stochastic grain distribution (b) Phenomenological failure (a) More realistic grain morphology (b) Grain is unbreakable (c) 2D model (a) Mineralogy is revealed (b) Grain can fail (c) Digital image-based model (a) Two-dimensional grain model (b) Grain parameters are difficult to obtain (c) High computation consumption (a) Grain can be broken (b) FDEM increases computation efficiency (c) Commercial code (a) Realistic reproduction of grain morphology (b) Grain is crushable (c) FDEM increases computation efficiency (a) Grain is unbreakable (b) Grain shape is overlooked (c) Mineralogy is overlooked

(a) Stochastic grain distribution (b) No grain boundary (c) Damage-based grain crushing

(a) Actual grain distribution (b) Grain crushing is available (c) Two-dimensional model;(d) mesh dependent (a) Two-dimensional model (b) Mesh size influences actual fracturing (a) No grain boundaries (b) Two-dimensional model

(a) 3D grain-based model

(a) Grain element is rigid (b) Two-dimensional model (c) Stochastic grain distribution

Advantages/disadvantages

[64] [30]

[74]

[28]

(continued on next page)

[50–51]

[1]

[40,29]

[94]

[47,79,20,19,31,81]

[17] [67]

[112] [18]

[49,63]

[98]

[68],

[69,70]

[39,90]

[44]

[45]

Literature

Table 1 Overview of grain-based modelling of geomaterials including application area, methodology, grain shape modelling, grain crushing ability, research results, advantages/disadvantages and references.

X.F. Li, et al.

Computers and Geotechnics 110 (2019) 161–183

Computers and Geotechnics 110 (2019) 161–183

[54,56,58]

[4,10–11,34–35,85,84,5,83,64]

heterogeneity concentrated on material heterogeneity rather than the fact of inter- or transgranular cracking; therefore, the fracturing behaviours are phenomenologically represented by element degradation indicators such as damage variables or plasticity flow. As shown in Fig. 2c, the grain model consists of a series of mapped mesh elements that are regarded as representative elementary volumes (REVs). These square elements are homogenous and isotropic and can also be used to construct irregularly shaped grains [107]. Two main codes used in the literature are FDM method and RFPA. Generally, the square-based grain method takes advantage of the continuum algorithm, for example, high computation efficiency and high accuracy, and indirectly addresses rock fracturing as element degradation or deactivation. With the aid of digital image processing (DIP) technique, micro heterogeneity was simulated in RFPA2D, and a statistical model in a Weibull-like form was developed to describe the homogenization of rocks [63]. This algorithm is different from the method proposed by Li et al. [49], which aimed to simulate the mineral grains directly. With a focus on the actual microstructures of rock materials, Chen et al. [18] applied DIP technique to FLAC2D and simulated rock fracturing in Brazilian disc testing. This method was extended to FLAC3D [17] and RFPA2D [106]. However, the continuum square-based method has shortcomings: (a) the intergranular cracking on the grain boundary is not revealed and (b) the fracturing of grains is dependent on the element damage state or plasticity rather on realistic separation and sliding behaviours among grains. The node-based SPH method is another grain model developed in the continuum context. SPH nodes are discretised by Voronoi diagram-based domains, and mineral heterogeneity is achieved by setting different material properties [67]. The failure of grains is similar to that in square-based method and is based on the smearing of REVs by degradation of element stiffness and strength. Block-based grain methodology is widely used in DEMs, which model fracturing as debondings of neighbouring blocks. In comparison, continuum methods attempt to characterize rock mass behaviours from an average equivalence viewpoint using the representative elementary volume technique. DEMs aim to explicitly reproduce realistic discontinuities and directly simulate the separation and sliding of elements. The block method is subdivided into a polygon/polyhedronbased grain method [50] and a triangle-based grain method [29] according to the grain morphology. The mineral grains can be modelled as deformable bodies [47,79,20,69,70,55] or rigid bodies [45,44]. In the rigid-body-spring network (RBSN) scheme, the contacts between Voronoi-blocks are regarded as springs, which provide average stresses on the corresponding lattice elements, as shown in Fig. 3. A rheological model with a visco-plastic damage model is used to consider the ratedependent failure features [45,44,39]. Considering the truth that mineral grains are deformable and crushable in experiments, deformable elements are commonly used for rock grains. Shown in Fig. 2e and f are the typical models of triangle-based grains and polygon-based grains, respectively. Fig. 2 illustrates a method based on the DIP technique. This result shows that the irregular grain morphology can be fully reproduced with the aid of DIP technique. Another method that can be used to simulate grain heterogeneity is the stochastic distribution [47,69]. Polygonal grains provide more actual representation of the grain morphology [29] but result in rougher surfaces on macro faults, which lead to unrealistically high values of the macro friction angle [74]. This outcome occurs because (a) the polygon-based grain model restrains shear sliding of grains under high deviatoric stress due to asperities and rough failure paths of grain boundaries and (b) asperities restrain the rotation of minerals. Even so, the polygon-based method is still popular since the Voronoi tessellation model was developed in commercial code UDEC [66]. For example, Lan et al. [47] used Voronoi tessellation to simulate the microstructure of brittle rocks and found that the crack initiation stress is governed by the microheterogeneity of grains, which leads to tensile stress concentration on grain boundaries. Nicksiar and Martin [79] employed the elastic mineral model to investigate the effect of grain heterogeneity on crack initiation. Park et al.

Strain rate mechanism and microfracturing transition Grain crushable

Particle clustered grain DEM; PFC-grain model

DEM; multiscale DEM

Compression test, direct shear test

SHPB test

Element plasticity, grain is crushable Dual-layer modelling, grain crushing is allowed Dual-layer Voronoi grains DEM; 3DEC-Voronoi Uniaxial compression test

Actual grain shape, dual-layer grains, particle-based method

Confinement mechanism is highly related to grain-scale fracturing

(a) Grain crushing is available (b) Stochastic grain distribution (c) Parameter calibration is difficult (d) Smooth joint model is used for grain boundary (a) Grain crushing is available (b) Realistic reproduction of grain shapes

[99]

[52,55]

(a) 3D grain model (b) No grain crushing (a) 3D stochastic grain distribution Grain structure on rock strengths Elastic minerals Polygon-based grains DEM; 3PDEM Compression and tension test

Table 1 (continued)

Grain methodology Method Application

Grain crushing ability

Results

Advantages/disadvantages

Literature

X.F. Li, et al.

164

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 2. A typical demonstration of realistic reproduction of grain heterogeneity on the basis of digital image processing methods (DIPs). (a). The original image of grain distribution; (b) digital image processing in combination with binary segmentation; reproductions of (c) square-based grains; (d) particle-based grains; (e) triangle-based grains; and (f) polygon-based grains (modified after [51]).

[81] considered the diverse mechanical properties of grain boundary contacts between different minerals in GBM-UDEC and obtained an encouraging strength ratio of tensile to compressive strength in the range of 1/20 to 1/10. In these studies, mineral grains are treated as deformable elements, while the overall deformation is attributed to deformation of mineral grains and contacts (as shown in Fig. 3). To extend this method to a three-dimensional model, Voronoi tessellation in combination with the Neper code [89] is implemented in 3DEC [30]. This novel method achieved 3D construction of minerals but failed to consider different minerals in brittle rocks. A similar problem is also encountered in the study of the 3PDEM method, which treats the grains as deformable/rigid polyhedra [55]. The two-scale Voronoi tessellation scheme was incorporated into 3DEC-GBM, therefore, inter- and intragrain cracking can be simulated [99]. Generally, the polygon-based grain method is incorporated with the Voronoi tessellation technique in UDEC, and the polyhedron-based grain can be generated in Neper code [89]. A similar block-based method in the family of grain models is that regarding the fundamental grain as a triangular element. The trianglebased grain method is primarily used in UDEC because of the simplicity of element generation [40]. Triangular grains contribute to smoother pathways that encourage shear sliding [30] and result in a 45% decrease in the uniaxial compressive strength (UCS) compared with Voronoi tessellation [74]. This result occurs because the triangular blocks display an increase in freedom and reduce the possibility of interlocking grains. This method was then extended to 3D models incorporating 3DEC code [28]. Given the importance of grain crushing in rocks under compressive stress, a grain-breakable model was developed in the context of Voronoi tessellation [29]. Dual-layer fracturing was

also achieved in 3D modelling with the aid of two-scale tessellation in Neper code [99]. Dual-layer DEMs have the advantage of characterizing the transgranular fracturing behaviours of rocks and can also be performed for grain heterogeneity modelling [29]. However, the intrinsic drawbacks of DEMs, including the rounding error, the low efficiency in contact detection and the difficulty in microparameter calibration, limit the promoted application of dual-layer DEMs on the in-situ scale. Inheriting both the advantages of the continuum method and the discontinuum method, the hybrid FDEM is widely used in rock engineering [77,61]. The Grasselli group and related researchers extended this method to develop a realistic grain shape-based method in open source code and commercial code such as Y-Geo and Irazu-GBM, [70,68,1]. As shown in Fig. 3, the logic of grain-based FDEMs is similar to that of the dual-layer DEM proposed by Gao et al. [29]. The mineral grains are constructed by an assembly of triangular meshes while the grain boundary is rather rough. The early grain-based FDEM in Y-Geo focused on rock heterogeneity modelling rather than characterizing the role of inter- and intragrain fracturing in rocks [1]. Taking advantage of the more realistic representation of grain structures by the Voronoi tessellation technique, polygon diagrams were adopted to replace the grain phases in conventional grain-based FDEM [1]. Fracturing behaviours were explicitly modelled based on the principles of nonlinear elastic fracture mechanics in Irazu-GBM, and the grain boundaries were differentially treated by reducing the mechanical parameters [1]. Similarly, Li et al. [50–51] proposed a grain-based model in FDEM, and several improvements are achieved by this method: for example, (a) grain boundaries are treated as initial discrete fracture networks; (b) the Mohr-Coulomb criterion is used for intergranular cracks; and (c) the

165

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 3. The polygon-based grains implemented in the DEM and FDEM methods. These methods can be roughly divided into rigid grains, deformable grains with grain-unbreakable ability and dual-layer elements.

grain heterogeneity is realistically reproduced in combination with the DIP technique. GB-FDEM also reduces the number of uncertain parameters used in grain-based models and simplifies the calibration procedure of DEMs [51]. Particle-based grain methodology is a large family in the area of rock modelling and provides some novel ideas in the development of grain-based models. The initial grain shape in particle-based DEMs is circular (2D) or spherical (3D), which is applicable to granular materials such as sand and soil [88]. The bonded particle models (BPM), such as the parallel bond model (PBM), contact bond model (CBM) and softening contact model, were proposed to mimic the cohesion loss, which is observed as acoustic emissions from brittle rocks in the laboratory [86]. In BPM, each particle is regarded as a mineral grain, and therefore, the grain deformability is overlooked as well as the grain morphology. Numerous studies have verified that the circular morphology of grains leads to low interlocking between grains, which reduces restraints on the rotation of particles [86,22,53]. Even if particle rotation is resisted by the construction of parallel bonds in PBM, which is used to mimic the mechanical behaviours of epoxy cement, as shown in Fig. 4b, the smooth grain surface still fails to provide enough resistance to rotation. Due to the circular shape of the grains, the applied compressive stress is easily transformed to local tensile stress because of

the large moment induced by grain rotation [30]. The BPM have three main shortcomings: (a) an unrealistically low strength ratio of UCS to Brazilian tensile strength (BTS) in comparison with laboratory results [22]; (b) a low macro friction angle and (c) approximately linear correlation of failure strength with confining pressures [86]. To solve these intrinsic problems, a kind of notional particle model has been developed to increase the rotation resistance [87]. As shown in Fig. 4a and b, the greatest improvement in the flat joint model (FJM) is considering the progressive damage to bonds. The contact is discretised into subsegments and is kept bonded until all segments are broken. Therefore, the rotational motion of particles is limited to a low level, which has a finite influence on the tensile fracturing of contacts. The FJM has success in mimicking the actual response of brittle rocks, especially the strength ratio and macro friction angle [103]. The drawbacks of this model are that: (a) it lacks an underlying mechanism for the sub-segment discretisation and (b) the simulated microcracks are pervasively distributed, which fails to agree with laboratory results. Irregular grain aggregates have been imposed to increase the grain interlocking and explicitly simulate a rough grain surface on a realistic scale [21,22,53]. The aggregate can be regarded as a deformable cluster or a rigid clump according to the freedom of an individual particle. A range of nonsmooth, non-spherical and non-convex agglomerates can be modelled

166

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

(a).

(c).

Fig. 4. An example of the flat joint contact model proposed in PFC for grain modelling. (a) Typical model used in PFC code (modified after [87]); (b) an example image of an epoxy-cemented glass bead sample with SEM photomicrographs (modified after [36]); (c) a schematic of the cluster model; and (d) a schematic of the clump model.

(b).

(d).

by assembling disks together to create deformable clusters [80]. The cluster model is very similar to BPM except that the mechanical properties of grain boundary contacts are weakened. Fig. 4c and d illustrate the schematics of the cluster and clump used in grain-based modelling. Each particle within the grain has its own freedom of motion, including translation and rotation. Therefore, the agglomerate can be translated, deformed and even crushed in the cluster model [54,58]. The irregular grain boundaries increase the grain interlocking and decrease the degree of rotation-induced tensile fracturing in the bonded model [53]. This improvement overcomes intrinsic drawbacks inherited from circular particles and highly increases the utility of particle DEMs in bonded materials. Numerical results indicate that the strength ratio (σucs/σt) can be increased from 4.76 for PBM to 11.1 for the cluster model [53]. The clump model has distinguishing characteristics that treat the entity as a rigid body, and the motion results from the translation and rotation of the whole aggregate [22,6]. Cho et al. [22] indicated that the clump-geometry model induces a significant improvement of the strength ratio to more than 14.3. Moreover, the clump model allows large overlapping between grains, and therefore, the initial porosity of specimens can be reduced to the actual level. Following the logic of the cluster model, Li et al. [54,58] comprehensively proposed an actual grain model using the DIP technique and then investigated the dynamic behaviours of granitic rock with attention to inter- and intragranular cracking. Deformable and crushable minerals seem to be more comparable with actual conditions and explain more mechanisms of grain-scale behaviours such as strain localization and grain pulverization [58]. Two main grain boundary models, namely, the smooth joint model (SJM) and PBM, are mainly used in clusterbased GBM. The smooth joint model is imposed to model a smooth

interface regardless of the local particle topology [73]. This model is widely used for initial joints that have planar surfaces and is used in GBM because the original code package provided by Itasca adopted this boundary model. Li et al. [58] realized that the grain boundary is not ideally planar in reality and therefore utilized PBM as grain boundaries in GBM. This modification is proven to be valid for two main reasons: (a) the bonded grain boundary has a large realistic rotation restraint that can be better simulated by BPM and (b) the planar approximation in SJM loses its effect when large deformation occurs [75]. 3. Grain-breakable models and fracturing behaviours As outlined above, many grain shapes have been adopted in the framework of grain-based modelling. This study intends to impose the block-based method and particle-based method to compare the differences between common grain methods. Benefiting from the advantages of the discontinuum context, grain-breakable and grain-unbreakable models are employed to investigate the mechanism of grain crushing on microfracturing for rocks. Therefore, four main methods, including a block-based grain-unbreakable model implemented in UDEC-GBM, a GB-FDEM allowing for grain breaking, a clump model in particle DEM and a cluster model, are used in this study. For simplicity, these models are abbreviated to UDEC-GBM, GB-FDEM, clump and cluster in the following sections. 3.1. Grain crushing and transgranular fracturing Grain crushing is a common failure phenomenon in granular materials under compressive stress [105]. For cemented materials, the

167

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 5. Photomicrographs of transgranular cracks and grain crushing encountered in brittle rocks. (a) Stress-induced microcracks in Kfsp and Qtz grains in granitic rocks (after [60]); (b) and (c) are crushed K-feldspar and quartz in rocks under high confining stress (modified after [43]); (d) coalescence of microcracks in shear bands of Westerly granite at a confining stress of 250 MPa (after [102]), and (e) localized shear zone in San Marcos gabbro [101].

term ‘grain crushing’ means grain pulverizations induced by assembled microcracks within grains, as shown in Fig. 5. In the Lac du Bonnet granite and Forsmark granite samples corded by means of solid disking and ring disking (as shown in Fig. 5a), the proportion of grain-crossing cracks ranges from 8% to 69% [60]. From the results for the Verzasca Gneiss after compression with high-level stress, the original grains in the gouge zone form wide and elongated arrays of fractured particles [43]. At the same time, the quartz and feldspar are fractured, while the feldspar is more intensively pulverized than quartz, as shown in Fig. 5b and c. The grain size distribution in fault zones is analysed, and the fractal value of the cracked grains is approximately 1.64. Early studies on compression of the Westerly granite as shown in Fig. 5d and the San Marcos gabbro in Fig. 5e found that the microcracks tend to form

inclined zones, which are defined as macro shear bands under higher compression. Most fracturing patterns of these microcracks are regarded as transgranular fracturing [82]. 3.2. Inter and transgranular fracturing models The failure of rock material is simulated by different contact models as shown in Fig. 6. In FDEM, the cohesive element is imposed to capture nonlinear fracturing behaviours at the tip of the crack, which is known as the fracture process zone (FPZ). The entity is constructed by an assembly of triangular elements, and the four-node crack element is introduced to present the fracturing process, as shown in Fig. 6b. GBFDEM originates from GB-DEM, which treats the grain boundary as an

168

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 6. Fracture behaviours of the element boundary used in FDEM. (a) FPZ model induced by local tensile stress (modified after [62]); (b) four-node element used for a smearing element during the yielding process; dual-layer constitutive models used in GB-FDEM in (c) the normal direction and (d) the tangential direction (modified after [51]); (e) FPZ model of crack element used in Y-GUI for mode I fracturing; (f) mode II sliding; and (g) mixed mode fracturing (modified after [62]).

169

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

initial distinct contact with reduced mechanical properties. Therefore, two different fracturing models exist in GB-FDEM, which is different from the conventional FDEM and Irazu-GBM models [50]. As shown in Fig. 6c and d, two intergranular crack models are simulated by the linear brittle fracture model: mode I, intergranular opening model, and mode II, intergranular shear sliding, and two transgranular crack models are simulated by FPZs as mode I, transgranular opening model, and mode II, transgranular shear sliding. In UDEC-GBM, only the intergranular crack models are employed for intergrain behaviours, which are described as: IG mode I fracturing: An intergranular tensile crack is simulated by the linear bond model of the Mohr-Coulomb criterion. The magnitude of the normal stress of contact is expressed as n:

=

n

ds =

t

1

GIIC =

s

1,

n tan IG + c s n tan IG

s

1,

<

where σc is the compressive strength and ϕIG is the frictional angle on the grain boundary. σs is the shear stress that can be computed from s:

=

s

(4)

ks, IG us

TG mode I fracturing: Transgranular tensile fracturing is characterized by the FPZ model in the form of a linear softening model in GB-FDEM [50]. This model is slightly different from the cohesion model imposed in FDEM, which uses a heuristic scaling function to represent the progressive loss of cohesion, as shown in Fig. 6e–f. In the normal component of the cohesive element, the bonding stress σTG,n can be computed as

kn, TG n, TG, n

= kn, TG (1

n

d ) n, 0.0,

<

p, n

p, n

n

>

n

r, n

(5)

r ,n

TG mode II fracturing: Mode II in transgranular cracks is simulated by a slipping model similar to that shown in Fig. 6d. The bonding tangential stress σTG,s is reduced based on the linear damage model, which can be denoted as

ks, TG s , TG,s

= ks, TG (1 TG , n

tan

<

s

d) s,

p, s

p, s

TG ,

s s

>

r ,s r ,s

p, s

n

,

p, s

s

r, n r, s

(7)

r,n

r,s p, s

p, n

(

TG, n d n

TG, s (s )

s, r ) d s

(8)

Although grain-based models have the advantage of describing the microheterogeneity of rocks, the difficulty of laboratory measurements on the microproperties of minerals limits the promotion of these methods. In DEMs, calibration is an essential step to match the microproperties to macro responses observed in the laboratory [74]. Most parameters are micro properties of rocks, which are difficult to obtain from laboratory tests, especially when the heterogeneity increases the mineral components in a GBM. To date, no standard calibration procedure has been proposed for GBMs because of the complexity of the interactions among grain geometry, size, morphology and material properties [50]. Back analysis based on the logic of trial and error is widely used in calibration. In addition, microindentation is an appropriate approach to obtain the micro mechanical properties of minerals directly. Table 2 summarizes the micro parameters used in different grain-based methods. Most models have more than five uncertain parameters, and this number increases exponentially as the mineral components increase. In this table, the term ‘uncertain parameters’ indicates that the parameters cannot be obtained directly in the laboratory; other macro parameters such as Young’s modulus, Poisson’s ratio and the fracture energy of minerals, which can be obtained in micro/ nanoindentation tests, are regarded as certain parameters. Even if the elastic properties can be approximately estimated from the micro indentation test, the difficulty of the testing operation still leads to obstacles in simulating the grain-scale behaviours of rocks with certainty. As outlined in the literature [51], two kinds of micro properties are distinguished in grain models: (i) continuum parameters and (ii) discontinuum parameters. Continuum parameters, such as the grain size, mineral density, Young’s modulus, Poisson’s ratio and fracture energy, can be directly measured in experimental testing or estimated from the literature [10]. To reduce the number of uncertain parameters, the calibration procedure involves only discontinuum parameters, such as the failure strength of minerals and stiffness of grain boundaries. In GBFDEM, the authors explored the possible bridges between micro parameters and macro behaviours and proposed a standard calibration procedure [58]. Unfortunately, this method is applied only for blockbased methods with straightforward fracturing features. The calibration procedures for particle-based methods, especially for clusters, have been reported in the literature [58]. Table 3 lists the micro parameters of the four models used in this study. The minerals in UDEC-GBM inherit the elastic parameters of GB-FDEM, while the fracturing behaviours are not considered because grains are unbreakable in UDECGBM. Similarly, the micro parameters of the grain boundaries in the clump models are the same as those in the cluster models.

(3)

c

r,s

p, n

4. Calibration procedure for grain-based models

c

s

p, s

,

The cohesion model of tensile/shear cracking is illustrated in Fig. 6c, and GIC, GIIC are energy release rates obtained from laboratory tests.

where σt is the tensile strength of the boundary contact. In UDEC-GBM, the tensile fracturing of an intergranular contact is the same as the criterion proposed above. IG mode II fracturing: Mode II failure on the grain boundary is modelled by a slipping model. The sliding process can be divided into two parts: the loss of cohesion and the frictional sliding after bond breakage. Therefore, the mode II intergranular cracking criterion can be written as

Fs =

s

GIC =

(2)

n

p, n

The displacements corresponding to peak strengths are = tTG /kn. TG = tTG lave/ E , p,s = sTG d / G . The residual displacements can be computed from

where kn,IG is the normal and tangential stiffness of the contacts and Δun is the elastic portion of the normal and displacement increments. In GBFDEM, the stiffness can be simplified as kn,IG = (Ei + Ej)/dijEiEj, where Ei and Ej are the elastic moduli of parent elements I and j, which are connected by contact ij, and dij is the distance between the centroids of two elements. The tensile yield criterion is expressed as

Ft =

p, n

r,n

p, n

(1)

+ kn, IG un

n

dn =

(6)

where d is the damage variable of the elementary boundary denoted as

170

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Table 2 Calibration parameters used in different GBMs. Method

Corresponding code

Micro parameters

Uncertain parameters

Calibration method

Square-based GBM

Flac, RFPA

5

Micro/nanoindentation testing; trial and error

Polygon-based RBSN

RNSN

3

Trial and error

Block-GBM, grain elastic

UDEC, 3DEC

5

Block-GBM, grain plastic

UDEC, 3DEC

Mineral: E, ν, σt, Grain boundary: Mineral: rigid Grain boundary: Mineral: E, ν Grain boundary: Mineral: E, ν, σt, Grain boundary:

UDEC-GBM, Grain-breakable

UDEC

10

Mineral: Micro/nanoindentation testing Grain boundary: trial and error Mineral: Micro/nanoindentation testing; trial and error Grain boundary: trial and error Micro/nanoindentation testing; trial and error

FDEM

Y-GUI

6

Micro/nanoindentation testing; trial and error

GB-FDEM

GB-FDEM

10

Irazu-GBM

Irazu

Cluster

PFC, Yade and so on

10

Mineral: Micro/nanoindentation testing; trial and error Grain boundary: trial and error Mineral: Micro/nanoindentation testing; trial and error Grain boundary: trial and error Micro/nanoindentation testing; trial and error

Clump

PFC, Yade and so on

5

Trial and error

c, ϕ continuous E, ν, τ, N, GI kn, ks, σIG,t, σIG,c, ϕIG c, ϕ kn, ks, σIG,t, σIG,c, ϕIG

Mineral: kTG,n, kTG,s, σTG,t, σTG,c, ϕTG Grain boundary: kIG,n, kIG,s, σIG,t, σIG,c, ϕIG Mineral: E, ν, σt, c, ϕ, GI, GII, pn, pt, pf Grain boundary: unified Mineral: E, ν, σt, c, ϕ, GI, GII Grain boundary: σIG,t, σIG,c, ϕIG Mineral: E, ν, σt, c, ϕ, GI, GII, pn, pt, pf Grain boundary: σIG,t, σIG,c, ϕIG, GIG,I, GIG,II, pIG,n, pIG,t, pIG,f Mineral: kTG,n, kTG,s, σTG,t, σTG,c, ϕTG Grain boundary: kIG,n, kIG,s, σIG,t, σIG,c, ϕIG Mineral: Rigid Grain boundary: kIG,n, kIG,s, σIG,t, σIG,c, ϕIG

8

14

The symbols used above are Young’s modulus, E (GPa); Poisson’s ratio, ν; tensile strength, σt (MPa); cohesion, c (MPa); friction angle, ϕ (degree); viscosity coefficients, τ, N; fracture energy, GI and GII (J m−2); normal stiffness, kn (GPa/m); tangential stiffness, ks (GPa/m); intergranular tensile strength, σIG,t (MPa); intergranular cohesion, σIG,c (MPa); intergranular friction angle, ϕIG (degree); normal, shear and fracture contact penalties, pn (GPa/m), pt (GPa/m), and pf (GPa). Table 3 Suggested micro parameters of different grain models used in this study. Mechanical parameter

Qtz GB-FDEM

Grain size, mm Density, g/cm3 Young’s modulus, GPa Poisson’s ratio Tensile strength, MPa Cohesive strength, MPa Mode I fracture energy, J/m2 Mode II fracture energy, J/m2 Friction coefficient Porosity, % Area ratio, %

2.65 0.12 61.0 185.0 1421.5 3553.7 0.51 0.0 26

Fsp Cluster 2.1–3.6 0.91 94.0 0.32 63.0 ± 6.0 118.0 ± 25.0 0.33 0.5 57

Ma

IG boundary

GB-FDEM

Cluster

GB-FDEM

Cluster

GB-FDEM

2.6

2.9–6.1 2.92 54.3

3.05

1.8–3.1 3.31 37.2

kn = 0.2E/zmin, ks = 0.2G/zmin

42.0 ± 4.0 70.0 ± 12.0

16.0 58.0

28 52

0.71 0.5

1.64

0.8

0.36 52.0 173.0 842.3 1853.5 0.80 0.0 17

54.0 ± 5 86.0 ± 16 0.46 0.5

36.0 142.0 1142.5 2856.5 0.75 0.0

Cluster

Note: E is the average Young’s modulus of adjacent grains; G is the average shear modulus of adjacent grains. zmin is the minimum zone edge length. For particle models, the average Young’s modulus E is E = E1E2/(E1 + E2), and zmin is the distance between the centroids of two particles.

Using empirical micro parameters of a kind of granitic rock [50,58] as listed in Table 3, four grain models are created to investigate the role of transgranular fracturing on rock deformations. The models have identical geometric sizes with 10.0 cm height and 5.0 cm width and are compressed by two plates with a thickness of 1.0 cm. Mineral discretisation is performed using DIP, and the mineral compositions are quartz (Qtz, 26%), feldspar (Fsp, 57%) and mica (Ma, 17%). Quartz has a size range from 2.1 to 3.6 mm, and the corresponding ranges of feldspar and mica are 2.9–6.1 mm and 1.8–3.1 mm, respectively. To reduce the mesh dependency of the FDEM and cluster models, the minimum number of meshes in each grain of GB-FDEM should be more than 20 [51], and the average number of particles in an individual cluster is 50.0.

5. Results and discussions 5.1. The effect of transgranular behaviours on crack stresses Laboratory studies of rock deformation have shown that the stressstrain curves of rocks undergo five regimes: (i) crack closure; (ii) linear elasticity; (iii) stable crack development after crack initiation; (iv) unstable crack development after crack coalescence; and (v) post-peak failure [12,72]. Four characteristic stresses, including the crack closure stress (CC stress), crack initiation stress (CI stress), crack damage stress (CD stress) and peak failure strength, are defined to represent the failure states under different stress levels for rocks. For example, the CC stress is usually observed in high-porosity rocks and is applied for the

171

100

σp 0.4

σci 0.0

50 0 0.0

0.2

0.4

-0.4 0.8

0.6

(b).

250 200

0 0.0

0.8

Axial stress, σ1 IG: tensile crack TG: tensile crack TG: shear crack IG: shear crack εv

σci

100

0.4

σcd

σp 0.0

50 0 0.0

0.2

0.0

σci 0.2

0.4

-0.4 0.8

0.6

Axial strain [%]

0.4

-0.4 0.6

(d). 250 200

Stress, (MPa)

150

0.4

100 50

Crack density; Volumetric strain [%]

Stress, (MPa)

200

σp

σcd

150

Axial strain [%]

(c). 250

0.8

Axial stress, σ1 IG: tensile crack IG: shear crack εv

0.8

Axial stress, σ1 IG: tensile crack IG: shear crack εv

0.4

150 100

σci

σcd

σp 0.0

50 0 0.0

Axial strain [%]

0.2

-0.4

0.4

Crack density; Volumetric strain [%]

150

0.8

Stress, (MPa)

Stress, (MPa)

200

Axial stress, σ1 IG: tensile crack TG: tensile crack TG: shear crack IG: shear crack σ cd εv

Crack density; Volumetric strain [%]

(a). 250

Crack density; Volumetric strain [%]

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Axial strain [%]

(e).

GB-FDEM

Cluster

UDEC-GBM

Clump

Crack density per cm2

(f).

0.25 0.20 0.15 0.10 0.05 0.00

GB-FDEM

Cluster

UDEC-GBM

Clump

Fig. 7. Numerical results for the four kinds of grain-based models, including (a) GB-FDEM; (b) UDEC-GBM; (c) cluster model; and (d) clump model. The fragment state at final state is presented in (e) in which the colour of particle-based models means different debris; (f) is the crack density distribution of the models used above, and the grid area used in this study is 1 cm × 1 cm.

172

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

pronounced dilatancy in rock samples, which is believed to be the main cause of significant dilatancy in grain-unbreakable models. The high interlocking force induced by initial rough surfaces on the grain boundary in the clump model increases the difficulty of tensile crack generation. This feature is the reason why the dilatancy in the clump model is much less significant than that in UDEC-GBM, even if these methods are similar in the modelling of grain crushing behaviours. As discussed in the study by Li et al. [50], micro fracture generation follows a sequence of IG:tensile crack → IG:shear cracks → TG:tenisle → TG:shear cracks in grain-breakable models. Intergranular cracks are formed due to local tensile stress concentrations [47]. IG:shear cracks are usually generated accompanying IG:tensile cracks and account for no more than 10% of the total micro cracks. As shown in Fig. 7a and c, the proportions of IG:shear cracks in the GB-FDEM and cluster models are 10.1% and 1.4%, respectively. An IG:shear crack is a result of selfsimilar propagation under the compression-shear stress state. This fracturing phenomenon disappears in particle-based models because grain sliding is transformed to particle rotation in PBM. TG:tensile fracturing is evident in both GB-FDEM (∼27.4%) and the cluster model (∼77.4%), as shown in Fig. 7a and c. A high proportion of TG:tensile cracks in the cluster model is attributed to two main reasons: (a) the intergranular contact number is much lower than that of transgranular contact and (b) the grain boundary sliding induces intragrain fracturing. CD stress is identified as the initiation of shear cracks within mineral grains. The stress levels of the four models mentioned above are 161, 148, 172 and 181 MPa, respectively. Shear cracking is an indicator of large-scale failure or grain crushing that results from crack interaction and coalescence [50]. This conclusion is subjective for particle-based models because shear-induced sliding in PBMs is easily transformed to tensile fracturing. Shown in Fig. 7c and d are failure states and micro crack distributions for GB-FDEM, UDEC-GBM, cluster and clump models. Models with unbreakable grains have more significant dilatation, especially for UDEC-GBM (the second figure in Fig. 7c). Although this macro fragment state is similar to laboratory observations, the micro cracks shown in Fig. 7f indicate that the sample is almost pervasively fractured. However, shear bands and damage localization can be reproduced by GB-FDEM cluster models. Shown in Fig. 8 are close-up figures of microcracks observed in the uniaxial compression tests of the four models. More than six types of fracturing mechanisms have been found in grain-scale fracturing, which are comprehensively discussed in Li et al. [50]. In this study, we focus on discussing brittle cracks that propagate far away but consist of single cracks and shear bands clustered in an assembly of microcracks. Fig. 8a–c present brittle cracks in clump, cluster and GB-FDEM, respectively. Because of the obstacles of unbreakable grains, the extensional fracture in the clump model is not as long as those in the cluster model or GB-FDEM. In addition, this kind of fracture is usually presented in the form of grain isolation, which means that all bonds surrounding one mineral are ruptured. Conversely, this fracture in the cluster model is a result of the interaction of an IG:tensile crack and TG:tensile crack. This failure mode well approximates the fractures observed in the laboratory, for example, tensile fracturing and compressive testing [50] of granitic rocks. A similar result is also obtained in GB-FDEM, as shown in Fig. 8c. Many experiments have indicated that under axisymmetric loading, rock samples fail by axial splitting or dilation bands, as shown in Fig. 7e [8]. This type of failure transforms to shear bands in the form of shear localization when the confining pressure is increased [91]. As shown in Fig. 8d–g, shear bands appear in the form of grain pulverization in either grain-breakable or grain-unbreakable models. However, the shear banding in the cluster model and GB-FDEM shows more localized behaviour within several grains and

Fig. 8. Typical grain-scale cracks including single crack and shear banding cracks observed in numerical results. Single crack in (a) clump model; (b) cluster model; and (c) GB-FDEM. Clustered micro cracks forming shear bands in (d) clump model; (e) UDEC-GBM; (f) cluster model; and (g) GB-FDEM.

study of permeability in rocks as well as disturbances caused by deep coring [25]. The crack initiation stress is an important point corresponding to the onset of dilatancy, which implies the beginning of nonlinear behaviours of rocks. This stress is approximately 0.3–0.5 of the peak strength for some low-porosity rocks, such as basalt, diabase and granite [14]. From the overview of laboratory results, the average CI stress ratios (σci/σp) are 0.49, 0.48 and 0.50 for igneous, metamorphic and sedimentary rocks, respectively [51]. Although these stress levels are much lower than the peak strength, the initiation of microcracks is widely used in predicting rock bursts [2], stabilizing underground excavations [72] and controlling blasting damage [104]. Fig. 7a–d show the stress-strain curves and the propagation of microcracks as functions of axial strain for the four grain models. Using definitions of CI and CD in the study by Li et al. [50], the stress levels of crack initiation are 109, 64, 142 and 138 MPa for the GB-FDEM, UDECGBM, cluster and clump models, respectively. The grain-unbreakable capacity is accepted to increase the dilatation of samples and result in a lower crack initiation stress level [22]. Especially for the UDEC-GBM model, as shown in Fig. 7b, the generation of tensile cracks already increases to a certain level (∼10%) before the appearance of shear cracks. Most tensile cracks, which are induced by local heterogeneous tensile stresses, are parallel to the maximum principal stress [50]. These arrangements rarely decrease the stiffness of rocks and cause less nonlinearity in rock deformation. However, the proliferation of micro tensile cracks perpendicular to the extensional direction leads to

173

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

(a). 0.6

Normalized crack number

(b).

Tensile crack Shear crack

0.6

0 MPa 10 MPa 20 MPa 30 MPa 40 MPa

(c).

IG:tensile crack TG:tensile crack

0.6

0 MPa 10 MPa 20 MPa 30 MPa 40 MPa

0.4

0.4

0.4

0.2

0.2

0.2

0.0 0

100

200

300

0.0

0

Normalized crack number

(e).

Tensile crack Shear crack

0.8

160

240

320

0.0

0 MPa 10 MPa 20 MPa 30 MPa 40 MPa

0

80

Axial stress, (MPa)

Axial stress, (MPa)

(d).

80

IG:shear crack TG:shear crack

(f).

0 MPa 10 MPa 20 MPa 30 MPa 40 MPa

0.6

0.6

240

320

Axial stress, (MPa)

IG: tensile crack TG: tensile crack

0.8

0 MPa 10 MPa 30 MPa 40 MPa

160

IG: shear crack TG: shear crack

0.10

0 MPa 10 MPa 20 MPa 30 MPa 40 MPa

0.08

0.06 0.4

0.4

0.04 0.2

0.2

0.0

0.02

0

250

500

Axial stress, (MPa)

750

0.0

0

100

200

300

Axial stress, (MPa)

400

0.00

0

100

200

300

400

Axial stress, (MPa)

Fig. 9. (a) Normalized intergranular crack number vs. axial stress under different confining states for UDEC-GBM; (b) Normalized IG:tensile crack and TG:tensile crack numbers; and (c) normalized IG:shear crack and TG:shear crack numbers vs. axial stress for GB-FDEM. (d) Normalized intergranular crack number vs. axial stress under different confining states for clump; (e) normalized IG:tensile crack and TG:tensile crack numbers; and (f) normalized IG:shear crack and TG:shear crack numbers vs. axial stress for cluster.

appears to agree better with the experimental results that have been reported for sandstones [9]. The propagation of micro cracks, as well as the transition between crack types, is an interesting research topic. Fig. 9 plots the normalized crack number change as a function of axial stress, in which the normalized crack number is defined as i

=

Ni n i=1

Ni

whereas n = 2 for clump and UDEC-GBM and n = 4 for cluster and GBFDEM. In the case of UDEC-GBM shown in Fig. 9a, the following results can be obtained: (a) the proportion of IG:tensile cracks decreases from 48% to 24% when the confining pressure increases from 0 to 40 MPa; (b) confinement has an inhibitory effect on the generation of tensile cracks; and (c) the damage level at peak failure decreases as the confining stress increases. Similar conclusions can be obtained for the GBFDEM model, as shown in Fig. 9b and c. The percentage of transgranular cracks is less than 15%, and this proportion is reduced to 5% when the confining pressure is 40 MPa. Because the crack initiation stress is regarded as the stress level corresponding to tensile crack initiation, this stress increases from 109 to 275 MPa when the confining

(9)

where χi is the normalized crack number of type i cracks and Ni is the total number of type i cracks. n is the number of all possible crack types,

174

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

R1

R2

R3 R4

100

3DEC-Trigon 3PDEM PBM Clump Cluster EPBM FJM GB-FDEM Multi GB-DEM PFC-GBM RBSN UDEC-Trigon(FPZ) UDEC-GBM

Tensile strength (MPa)

Sandstone Coal

10

from that in other models. These differences are mainly embodied in three aspects: (a) tensile cracks, either IG:tensile cracks or TG:tensile cracks, account for the dominant proportion, and in general, the percentage of shear cracks is less than 10%; (b) no evident transition of micro fracturing is observed, although a slight increase in TG:shear cracks from 3% to 9% occurs when the confining pressure is increased; and (c) the total crack number is greatly increased as the confining stress increases.

R5 R6

1

5.2. Effect of transgranular fracturing on strength ratio σucs/σt

0.1

0.01 0.1

1

10

100

Considering the difficulty of measuring the tensile strength of rocks in the laboratory, the strength ratio, which is denoted as the ratio of UCS to tensile strength, is an important parameter for estimating the tensile strength in the field [16]. The strength ratio is considered an intrinsic index of rock materials that is dependent on rock type. Sheorey [92] indicated that this ratio ranges from 2.7 to 39 for different rocks, and the average value is approximately 14.7. Brook [15] also concluded that the mean ratios for mudstone, limestone, sandstone and granite are 10, 10, 15 and 20, respectively. As discussed in Potyondy and Cundall [86], the conventional PBM has the shortcomings of overpredicting the tensile strength and underpredicting the internal friction angle for cemented rocks. A grain-based model is a valid approach to overcome this defect and mimic a more realistic strength ratio [86,22]. Shown in Fig. 10a is an overview of possible simulated strengths from different grain models. PBM fails to simulate a strength ratio of more than 5.0 without modifying contact models. Extended PBM (EPBM), which modifies the contribution coefficient of the moment on tension fractures, enlarges the strength ratio range [97]. Generally, non-spherical models such as 3DEC-Trigon, 3PDEM, UDEC and RBSN can successfully capture the unilateral effect of rocks, and the possible simulated strength ratio is discretionary. As shown in the figure, clump, cluster and FJM overcome this shortcoming in PBM, especially for the clump model, and the possible strength ratio can be more than 100.0. The cluster model also has a flexible range from 3.0 to 50.0 for the strength ratio [53]. Using a Mohr-Coulomb failure envelope, the relationship between tensile strength and UCS is constant, which is determined by the friction angle φ as

1000

UCS (MPa)

(b). 250 Stress, (MPa)

200

Clump Cluster UDEC-GBM GB-FDEM

150 100 Compression test

50 Extension test

0 0.0

0.2

0.4

0.6

0.8

Axial strain [%] Fig. 10. (a). A summary of the simulated strength ratios of UCS to tensile strength using grain-based methods: 3DEC-Trigon [28]; 3PDEM [55]; PBM [86,100,93,53]; clump [22,6,53]; cluster [53]; EPBM [97]; FJM [97,103]; GBFDEM [51]; Multi GB-DEM [58]; PFC-GBM [84–85,11,34–35]; RBSN [90]; UDEC-Trigon [41–42] and UDEC-GBM [47,81,20,29]; (b) the simulated results of compressive and tension test of for different grain-based models.

stress increases from 0 to 40 MPa. From the results of the clump model shown in Fig. 9d, the transition from tensile crack to shear crack is more pronounced when the confining pressure is increased. Due to rotationinduced tensile failure in PBM, the fracturing result is rather different

t

(b).

50 40 30

20

20

15 10

10

5

0

3

0

10

20

30

40

50

60

=1(

1

3)

(10)

where σ1 and σ3 are the major and minor principal stresses and σucs and

Proability (%)

Macro tensile strength, σ t (MPa)

(a).

1 ucs

Macro UCS, σucs (MPa)

(a).

250 200 150 100 50 0

0

UDEC-GBM GB-FDEM Clump Cluster PBM

0

10

20

30

40

50

60

Micro tensile strength, σIG,t (MPa)

Micro tensile strength, σ IG,t (MPa)

Fig. 11. Numerical results for (a) macro tensile strength and (b) macro uniaxial compressive strength as functions of micro tensile strength for different grain-based models. (a) Is the probability contour of the relationship between macro tensile strength with micro intergranular tensile strength. The dots are numerical results of different samples.

175

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

(a). Stress, (MPa)

600

40

400

30 200

20 30 40

10 0 0.5

0 0.0

10

1.0

1.5

2.0

2.5

Clump Cluster UDEC-GBM GB-FDEM PBM Laboratory result

Axial strain [%] GB-FDEM UDEC-GBM

280

Stress, (MPa)

Table 4 Computed results of clump, cluster, UDEC-GBM, GB-FDEM, PBM and laboratory results.

50

0

(b). 350

20 10

210

40

0 20

140

30

ucs

10

450

1.0

1.5

150

Internal friction angle, (°)

(d).

Sedimentary rock -20

0

20

40

σ 3, (MPa) 80

Clump Cluster UDEC-GBM GB-FDEM PBM Labratory test

60

Friction decrease 40

29.3 35.3 29.9 39.8 45 28.2

57.2 43.5 49.8 35.5 29.8 51

t

=

2c cos 1 + sin

(11)

1 + sin 1 sin

(12)

The Hoek-Brown criterion is a kind of nonlinear strength envelope and is widely used to predict rock strength under different stress states [33]. The equation chosen for the failure of intact rock is

Compaction

20

0

202 183 185.6 129.4 168.4 169

5.3. Hoek-Brown parameter mi

Extension

0 -10

40 10.8 12.9 11.1 3.2 9.5

The friction angle should be more than 55° to satisfy the strength ratio range σucs/σt > 10 from Eq. (12). This large value fails to agree with experimental results for the friction angle in triaxial compression tests. Fig. 10b shows the simulated stress-strain curves of the compression test and extension test with the four models. The strength ratios are 14.2, 11.25, 12.5 and 17.8 for clump, cluster, UDEC-GBM and GB-FDEM models, respectively. Lacking the interlocking force is the main cause of the low strength ratio in PBM, and the roughness of the grain morphology increases the resistance to particle rotation [53]. The fundamental reason is the use of the beam contact model, which adds the rotation term to the tensile stress [86]. Therefore, not only the macro tensile strength but also the macro UCS is significantly dependent on the micro tensile strengths. Fig. 11 shows the dependency of macro strengths on the micro tensile strength for the four models. In the case of macro tensile strength (as shown in Fig. 11a), σt exhibits a linear relationship with σIG,t. Most calibrated tensile strengths of contacts on grain boundaries vary in the range of 15–27 MPa in grain-based models. The change of UCS to σIG,t is plotted in Fig. 11b. The difference between PBM and grain models is the high dependency of UCS on intergranular tensile strength. The slope of the linear relationship against σIG,t is approximately 3.75 for PBM, and this value is reduced to ∼1.5 for cluster and ∼0.2 for UDEC-GBM.

Igneous rock

-40

=

t

300

0

ϕ (degree)

2c cos , 1 sin

=

ucs

0.5

Clump Cluster UDEC-GBM GB-FDEM PBM Labratory test

600

c (MPa)

where c is cohesion, and φ is the friction angle. Therefore, the strength ratio σucs/σt can be expressed as

Axial strain [%]

(c). 750

σucs (MPa)

σt are the uniaxial compressive strength and tensile strength, respectively. σucs and σt can be written as

40

0

0 0.0

mi

30

70

σ 1, (MPa)

Fig. 12. Stress-strain curves of (a) clump and cluster models; (b) GB-FDEM and UDEC-GBM models. (c) is the Hoek-Brown failure envelope for simulated experimental results. The labels ①, ②, ③, ④ and ⑤ are Hoek-Brown envelopes with mi parameters of 40, 12.9, 10.8, 11.1 and 3.2, respectively. The corresponding parameters σucs are 202, 185.6, 183, 129.4 and 168.4 Mpa, respectively. (d) Is the degradation of internal friction angle computed from the Hoek-Brown envelopes as listed in Table 4.

Cluster model Clump model

10

20

30

40

50

1

σ 3, (MPa)

=

3

+

ucs

mi

3

+1

where mi is the material constant, which equals mi = tan

176

(13)

ucs

ucs /

t

in the

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

derivation of Eq. (13) considering the growth of microcracks [110], where κ is the coefficient for mixed mode fractures. Another form of mi can be obtained when σ1 = 0.0 in Eq. (13) as

mi = ((

t ucs

)2

1)

ucs

pressures. This feature is the possible reason why the actual peak strength does not follow the linear increase described by the MohrCoulomb criterion. Fig. 12c plots the numerical and experimental results of granitic rocks under a confining stress range from σ3 = 0 to 40 MPa. HoekBrown and Mohr-Coulomb strength envelopes are fitted for different models as listed in Table 4. The material parameter mi is very high (∼40) for the clump model because the strength enhancement in confined states is dramatic. Unfortunately, this model results in a rapid drop of strength in the negative stress zone (under the extension stress state) and underestimates the tensile strength, especially for the tension test (∼8.0 MPa). The linear model of the Mohr-Coulomb envelope coincides well with these abnormal results, but the macro friction angle is as high as 57.2˚. Other grain models present appropriate predictions on the failure strengths excluding the PBM. The material parameter mi computed from laboratory results is 9.5, while this value of PBM is only 3.2. From the results shown in Fig. 12c, PBM has larger tensile strength and a lower friction coefficient in comparison with other models. These

(14)

t

Eq. (14) can be simplified as mi = σucs/σt when the strength ratio is larger than 8.0 [16]. The Hoek-Brown parameter mi is important to describe the nonlinear shape of the strength envelope, especially for the nonlinear increase under high confining pressure. Fig. 12 shows the results of stress-strain curves for different models under different confining states. The curves illustrate that confinement has a strengthening effect on failure strengths as well as fracture strains. Because of the incapacity of grain crushing, the clump model exhibits the strongest confinement dependency, as shown in Fig. 12a. UDEC-GBM has lower peak strength than the clump model due to the deformability of mineral grains. Grain-breakable models allow for fracturing transition in highly confined samples and therefore show less dependency on confining

Table 5 Computed crack density results of clump, cluster, UDEC-GBM, GB-FDEM, and PBM and laboratory results. Each column has the identical contour legend as shown in row #5. σ3

GB-FDEM

UDEC-GBM

Cluster

Clump

Laboratory test

0 MPa

10 MPa

20 MPa

(continued on next page)

177

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Table 5 (continued) σ3

GB-FDEM

UDEC-GBM

Cluster

Clump

Laboratory test

30 MPa

40 MPa

Crack density per cm2

Crack density per cm2

Crack density per cm2

0.00 0.05 0.10 0.15 0.20 0.25

0.00

0.02

0.04

0.0

0.06

results agree well with linear strength envelopes such as soft rocks and sedimentary rocks, etc. On the other aspect, grain-breakable capacity has the ability to decrease the value of mi to an actual range. This conclusion reveals the underlying mechanism of the important role of grain crushing on the confinement effect. At the same time, the deformability of mineral grains is important to characterize the failure of rocks. Although grain breakage is not allowed in UDEC-GBM, the large deformation of grains can also reduce mi to the normal range. In the Mohr-Coulomb model, as presented in Eq. (10), the strength ratio is the slope of the strength envelope in principal stress space (σ1–σ3). This value is nonlinear in the Hoek-Brown model, and here, we use the secant line to represent the friction angle in a stress increment Δσ3. Therefore, the friction angle can be expressed as

d d

1 3

=

1 + sin 1 sin

=1+

3 ucs

+ 1)

1/2

0.3

0.4

0.5

0.0

0.1

0.2

5.4. Shear banding and grain pulverization As discussed above, grain-breakable models have the ability to simulate the high strength ratio of rocks as well as the nonlinear enhancement of peak strengths under confined states. Previous studies have indicated that natural cohesive materials subject to laboratory loading exhibit complex pressure-dependent stress-strain behaviours and failure modes [82]. Axial splitting and dilation bands are commonly observed for samples compressed under low confinement [111]. At higher confining stress, a macroscopic fault develops in the form of shear bands. From a continuum perspective, shear banding is attributed to plasticity localization, which is accompanied by dilatancy and plastic flow. In compression testing on sandstones, a compacting zone with pervasive grain crushing is observed when the sample is subject to high confining pressure (∼28 MPa) [26]. Natural faults also reveal that the distribution of crushed grains follows a power law, and the fractal value D ranges between approximately 1.8 and 5.5 [71]. Table 5 illustrates contours of the normalized crack density for the four models under confining pressures from 0 to 40 MPa. The ultimate fracture states of granitic rocks in the laboratory are presented in the last column. In comparison, GB-FDEM and cluster models successfully mimic the damage localization and shear bands in the specimens, while grain-

(15)

1 mi (mi 2

0.2

friction angle is limited to an unrealistically high level because of the overestimated interlocking force. Conversely, due to the lack of sufficient interlocking force, PBM exhibits a slight degradation of friction angle; moreover, the mean friction angle is under a low level.

where the friction angle can be implicitly expressed by substituting Eq. (13) into Eq. (15) as

1 + sin 1 sin

0.1

Crack density per cm2

(16)

Fig. 12d plots the change in friction angle as a function of confining pressure for different models. The friction coefficients of rocks under different confinement states are not identical. The loss of friction associated with an increase in confining pressure leads to a nonlinear strength envelope, as shown in the Hoek-Brown criterion. Although the clump model can also simulate a nonlinear decrease in friction, the

178

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

Fig. 13. Examples of ultimate fragment states under two conditions when σ3 = 0.0 MPa and σ3 = 30.0 MPa for the four grain-based models including cluster, GBFDEM, clump and UDEC-GBM. The colour in cluster and clump models indicates the ID of debris.

unbreakable models demonstrate pervasively distributed micro fractures with increasing confining pressure. To compare the macro failure states between different models, the fragment states are plotted in Fig. 13 for cluster, GB-FDEM, clump and UDEC-GBM models under uniaxial compression and confining conditions (σ3 = 30 MPa). When σ3 = 0 MPa, sliding-induced macro faults in the cluster are evident. In the cluster and clump models, the different colours represent the IDs of the fragment. In contrast to other models, the UDEC-GBM model fails in the form of evident dilatancy, which behaves as axial splitting along the grain boundaries. When the confining pressure increases to 30 MPa, clear changes can be observed: (a) dilatancy is no more pronounced; (b) grain-unbreakable models do not show shear localizations as we expected; and (c) more debris is generated in the faults of the cluster model (the black colour of fragments means this debris contains only one single particle, which is similar to powders in the laboratory). An additional effect of confinement is that the shear band angle decreases with increasing confining pressure [82,111]. In this study, the shear band inclination is defined as the angle between the shear band

normal and the principal stress. Because grain-unbreakable models fail to reveal significant shear bands, only the cluster model and GB-FDEM model are discussed in this section. The theoretical shear band inclination according to the Mohr-Coulomb criterion is

=

4

+

2

(17)

Fig. 14a plots the shear band angle as a function of confining stress. The solid lines are theoretical results derived from Eq. (17) when mi equals 20, 15, 10, 5 and 2.0. As expected, increasing confinement results in a reduction in the shear band angle both for the cluster and GBFDEM models. Grain crushing is an important factor affecting the confinement dependency on dilatancy and shear banding. In the cluster model, the decrease in θ is slight because the particles are treated as rigid entities. This assumption decreases the deformability of mineral grains and fails to simulate compaction bands in samples when the confinement is very high. This shortcoming induces higher θ magnitude and wider shear bands than laboratory results (shown in Fig. 15a and b). In contrast, a compaction band is found in GB-FDEM when the

179

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

(a).

Fig. 14. Shear-induced bands of grainbreakage models and the results for grain crushing number. (a) Plot of shear band angle vs. confining stress and (b) grain crushing number as a function of confining stress in cluster and GBFDEM. The labels ①, ②, ③, ④ and ⑤ are theoretical results when the mi parameters are equal to 20, 15, 10, 5 and 2, respectively.

(b).

confining pressure exceeds 30 MPa, and a negative dilatancy angle (θ < 45°) is considered a result of large deformations of mineral grains. At the same time, double shear bands are commonly observed in GB-FDEM, as shown in Fig. 15d. Plotted in Fig. 14b is the change in grain crushing number with confining pressure. Increasing confinement induces more crushed mineral grains during the formation of shear bands. Similar results were comprehensively discussed in previous studies from the aspect of laboratory testing [82,71]. Unfortunately, grain-unbreakable models such as clump and UDEC-GBM appear to have insufficient capability to simulate this mechanism. The simulated results in Fig. 14b suggest that the rigidity assumption of particles overestimates grain crushing and thus fails to capture the compaction phenomenon in rocks under high confining pressure. This shortcoming can be partially overcome by increasing the average particle number in each mineral grain in the cluster.

cluster (particle-based, grain is breakable) and clump (particle-based, grain is rigid), are compared in this study. The reasonability of grain models in combination with crack stresses, strength ratios, Hoek-Brown parameters and shear localization is assessed. The main conclusions are as follows: 1. The clump model treats mineral grains as rigid entities and results in evident dilatancy in the uniaxial compression test. Rough grain boundaries and unbreakable mineral grains both increase interlocking forces among particles and encourage a higher strength ratio in comparison with PBM. This character leads to clear confinement dependency of the peak strength, and the strength envelope shape follows a linear increase in the form of the Mohr-Coulomb criterion. The tensile strength simulated in the clump model is lower than those from experimental results. Due to the inability to crush grains, pervasively distributed cracks appear in samples subject to high confining pressure. This result indicates that the transgranular fracturing significantly affects the mechanical behaviour of rocks after crack initiation. 2. The advantages of the cluster model are that mineral grains are deformable and crushable compared with clumps. Transgranular contacts commonly have stronger mechanical and physical properties than intergranular cracks; therefore, the high strength ratio, nonlinear envelope and fracturing transition under high confinement can be fully revealed. The cluster model successfully reproduces most deformation behaviours, state transformations and underlying mechanisms of brittle rocks in this study. For example, the possible strength ratios by the cluster model span a range from 3.0 to 50.0; the locus of failure strength follows the Hoek-Brown criterion, and shear localization is formed as a result of grain crushing. Special factors that should be noted for the cluster model are (a) the difficulty of selecting proper micro parameters and (b) the lack of compaction bands under very high confinement. The

6. Summary and conclusions Grain-based models have been widely developed to simulate the micro heterogeneity of concretes, rocks and granular materials. These methods can be distinguished as square-based, triangle-based, polygonbased and particle-based grains according to different grain shapes and can also be categorized as grain-breakable and grain-unbreakable models in terms of the capacity for grain crushing. To systematically understand the features of different methods, Section 2 presents an overview of state-of-the-art of grain models, and the summarized contents are listed in Table 1 from aspects such as application area, methodology, grain shape modelling, grain crushing ability, and advantages/disadvantages. Subsequently, intergranular and transgranular fracturing models and some difficulties in calibration procedures are discussed. To investigate the roles of grain shape and grain crushing in rock deformation, four typical models, namely, GB-FDEM (block-based, grain is breakable), UDEC-DEM (block-based, grain is unbreakable),

180

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

(a).

(b).

σ3=0 MPa

σ3=40 MPa

8.0

Band width

0

80.

0

40.

4.8

0.0

1.6

.0

-40

.0

-1.6

-80

0.0

-12

-4.8

0.0

-16

0.0

-8.0

-20

(c).

(d). 0

0

10.

40.

5.0

20.

0

Band width

0.0

0.0 .0

-5.0

-20 0

. -10

Double shear

.0

-40

.0

-60

0

-80

.0

-15

.0

. -20

Fig. 15. Typical shear banding for the cluster and GB-FDEM models under uniaxial compression and confined state σ3 = 40.0 MPa. (a) Cluster model when σ3 = 0.0 MPa; (b) cluster model when σ3 = 40.0 MPa; (c) GB-FDEM model when σ3 = 0.0 MPa; and (d) GB-FDEM model when σ3 = 40.0 MPa.

latter point is very important in confirming the average particle number in each mineral grain because the grain deformability is governed by transgranular contacts. 3. Ignoring the high efficiency of contact detection for spherical/circular particles, polygonal grains are more appropriate for realistic rock modelling. UDEC-GBM has the advantages of mimicking most mechanical behaviours of rocks measured in the laboratory except for shear banding. Unbreakable grains lead to abnormal dilatancy under uniaxial compression. 4. Dual-layer fracturing models are implemented in GB-FDEM to

distinguish the reduction in physical properties on grain boundaries. Hybrid methods inherit the advantages of both FEM and DEM. A comparison of shear bands with those for the cluster model suggests that compaction bands occur when the confining pressure is very high. This phenomenon can be attributed to the deformability of mineral grains in the direction perpendicular to the major principal stress. At the same time, the shear band inclination rapidly decreases, and the dilatancy angle becomes negative. Double shear bands are accompanied by the appearance of compaction, as well as grain crushing. 181

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

From the crack propagation in association with increased deviatoric stress, the definition of crack initiation can be explained by the generation of intergranular cracks. For grain-unbreakable models, intergranular tensile cracks are initiated prior to shear cracks. Shear cracks, especially transgranular shear cracks, are indicators of crack interaction, propagation and coalescence. Unfortunately, this conclusion is not appropriate for particle models because tensile fracturing is a result of the interaction of rotation and opening in the beam contact model.

sandstone. Int J Rock Mech Min Sci 2002;39(7):917–32. [27] Fujii Y, Takemura T, Takahashi M, et al. Surface features of uniaxial tensile fractures and their relation to rock anisotropy in Inada granite. Int J Rock Mech Min Sci 2007;44(1):98–107. [28] Gao FQ, 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. [29] 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–17. [30] Ghazvinian E, Diederichs MS, 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–21. [31] Gui YL, Zhao ZY, Ji J, et al. The grain effect of intact rock modelling using discrete element method with Voronoi grains. Géotech Lett 2016;6(2):136–43. [32] Guo N, Zhao J. 3D multiscale modeling of strain localization in granular media. Comput Geotech 2016;80:360–72. [33] Hoek E, Martin CD. Fracture initiation and propagation in intact rock–a review. J Rock Mech Geotech Eng 2014;6(4):287–300. [34] Hofmann H, Babadagli T, Yoon JS, et al. A grain based modeling study of mineralogical factors affecting strength, elastic behavior and micro fracture development during compression tests in granites. Eng Fract Mech 2015;147:261–75. [35] Hofmann H, Babadagli T, Zimmermann G. A grain based modeling study of fracture branching during compression tests in granites. Int J Rock Mech Min Sci 2015;77:152–62. [36] Holt RM, Kjølaas J, Larsen I, et al. Comparison between controlled laboratory experiments and discrete particle simulations of the mechanical behaviour of rock. Int J Rock Mech Min Sci 2005;42(7–8):985–95. [37] Horii H, Nemat-Nasser S. Compression-induced microcrack growth in brittle solids: axial splitting and shear failure. J Geophys Res Solid Earth 1985;90(B4):3105–25. [38] Huang JY, Hu SS, Xu SL, et al. Fractal crushing of granular materials under confined compression at different strain rates. Int J Impact Eng 2017;106:259–65. [39] Hwang YK, Lim YM. Validation of three-dimensional irregular lattice model for concrete failure mode simulations under impact loads. Eng Fract Mech 2017;169:109–27. [40] Kazerani T, Zhao J. Micromechanical parameters in bonded particle method for modelling of brittle material failure. Int J Numer Anal Meth Geomech 2010;34(18):1877–95. [41] Kazerani T. A discontinuum-based model to simulate compressive and tensile failure in sedimentary rock. J Rock Mech Geotech Eng 2013;5(5):378–88. [42] Kazerani T. Effect of micromechanical parameters of microstructure on compressive and tensile failure process of rock. Int J Rock Mech Min Sci 2013;64:44–55. [43] Keulen N, Heilbronner R, Stünitz H, et al. Grain size distributions of fault rocks: A comparison between experimentally and naturally deformed granitoids. J Struct Geol 2007;29(8):1282–300. [44] Kim K, Bolander JE, Lim YM. Failure simulation of RC structures under highly dynamic conditions using random lattice models. Comput Struct 2013;125:127–36. [45] Kim K, Lim YM. Simulation of rate dependent fracture in concrete using an irregular lattice model. Cem Concr Compos 2011;33(9):949–55. [46] Kranz RL. Microcracks in rocks: a review. Tectonophysics 1983;100(1–3):449–80. [47] 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). [48] Lemrich L, Carmeliet J, Johnson PA, et al. Dynamic induced softening in frictional granular materials investigated by discrete-element-method simulation. Phys Rev E 2017;96(6):062901. [49] Li L, Lee PKK, Tsui Y, et al. Failure process of granite. Int J Geomech 2003;3(1):84–98. [50] Li XF, Li HB, Zhao J. Inter/transgranular fracturing of micro-heterogeneous rocks and its effect on rock strength. Int J Numer Anal Met 2018. (in preparation). [51] Li XF, Li HB, Zhao J. Numerical investigation of crack initiation and propagation based on grain breakable FDEM for brittle rocks. Int J Rock Mech Min Sci 2018. (in preparation). [52] Li XF, Li HB, Li JC, et al. Crack initiation and propagation simulation for polycrystalline-based brittle rock utilizing three dimensional distinct element method. Rock Dyn: Res Eng 2016:341–8. [53] Li XF, Li HB, Liu YQ, et al. Numerical simulation of rock fragmentation mechanisms subject to wedge penetration for TBMs. Tunn Undergr Space Technol 2016;53:96–108. [54] Li XF, Li HB, Zhang QB, et al. Dynamic fragmentation of rock material: Characteristic size, fragment distribution and pulverization law. Eng Fract Mech 2018;199:739–59. [55] Li XF, Li HB, Zhao J. 3D polycrystalline discrete element method (3PDEM) for simulation of crack initiation and propagation in granular rock. Comput Geotech 2017;90:96–112. [56] Li XF, Li X, Li HB, et al. Dynamic tensile behaviours of heterogeneous rocks: The grain scale fracturing characteristics on strength and fragmentation. Int J Impact Eng 2018;118:98–118. [57] Li XF, Li HB, Zhang QB, Zhao J. Research on the dynamic behaviour of rock material: rate dependency and pulverization. Rock dynamics and applications, 3. CRC Press; 2018. p. 265–70. [58] Li XF, Zhang QB, Li HB, et al. Grain-based discrete element method (GB-DEM) modelling of multi-scale fracturing in rocks under dynamic loading. Rock Mech Rock Eng 2018:1–33. [59] Liakas S, O'Sullivan C, Saroglou C. Influence of heterogeneity on rock strength and stiffness using discrete element method and parallel bond model. J Rock Mech

Acknowledgements The editorial advice and suggestions from anonymous reviewers are acknowledged. The authors would like to thank the open source code YADE (https://yade-dem.org/doc/). We particularly acknowledge the finical supports from National Natural Science Foundation of China (No., 41572307) and China Scholarship Council (to the first author, No. 201604910678). References [1] Abdelaziz A, Zhao Q, Grasselli G. Grain based modelling of rocks using the combined finite-discrete element method. Comput Geotech 2018;103:73–81. [2] Alcott Jane M, Kaiser Peter K, Simser Brad P. Use of microseismic source parameters for rockburst hazard assessment. Seismicity caused by mines, fluid injections, reservoirs, and oil extraction. Basel: Birkhäuser; 1998. p. 41–65. [3] Anders MH, Laubach SE, Scholz CH. Microfractures: a review. J Struct Geol 2014;69:377–94. [4] 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–30. [5] Bahrani N, Kaiser PK. Estimation of confined peak strength of crack-damaged rocks. Rock Mech Rock Eng 2017;50(2):309–26. [6] Bahrani N, Valley B, Kaiser PK. Numerical simulation of drilling-induced core damage and its influence on mechanical properties of rocks under unconfined condition. Int J Rock Mech Min Sci 2015;80:40–50. [7] Bai QS, Tu SH, Zhang C. DEM investigation of the fracture mechanism of rock disc containing hole (s) and its influence on tensile strength. Theor Appl Fract Mech 2016;86:197–216. [8] Bésuelle P, Desrues J, Raynaud S. Experimental characterisation of the localisation phenomenon inside a Vosges sandstone in a triaxial cell. Int J Rock Mech Min Sci 2000;37(8):1223–37. [9] Bésuelle P. Compacting and dilating shear bands in porous rock: theoretical and experimental conditions. J Geophys Res Solid Earth 2001;106(B7):13435–42. [10] Bewick RP, Kaiser PK, Bawden WF, et al. DEM simulation of direct shear: 1. Rupture under constant normal stress boundary conditions. Rock Mech Rock Eng 2014;47(5):1647–71. [11] Bewick RP, Kaiser PK, Bawden WF. DEM simulation of direct shear: 2. Grain boundary and mineral grain strength component influence on shear rupture. Rock Mech Rock Eng 2014;47(5):1673–92. [12] Bieniawski ZT. Mechanism of brittle fracture of rock: part II—experimental studies. Int J Rock Mech Min Sci Geomech Abstr Pergamon 1967;4(4):407–23. [13] Bobet A, Einstein HH. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci 1998;35(7):863–88. [14] Brace WF, Byerlee JD. Recent experimental studies of brittle fracture of rocks. The 8th US symposium on rock mechanics (USRMS). American Rock Mechanics Association 1966. [15] Brook N. The measurement and estimation of basic rock strength. In: Rock testing and site characterization; 1995. p. 41–66. [16] Cai M. Practical estimates of tensile strength and Hoek-Brown strength parameter m i of brittle rocks. Rock Mech Rock Eng 2010;43(2):167–84. [17] Chen S, Yue ZQ, Tham LG. Digital image based approach for three-dimensional mechanical analysis of heterogeneous rocks. Rock Mech Rock Eng 2007;40(2):145. [18] Chen S, Yue ZQ, Tham LG. Digital image-based numerical modeling method for prediction of inhomogeneous rock failure. Int J Rock Mech Min Sci 2004;41(6):939–57. [19] Chen W, Konietzky H, Tan X, et al. Pre-failure damage analysis for brittle rocks under triaxial compression. Comput Geotech 2016;74:45–55. [20] Chen W, Konietzky H. Simulation of heterogeneity, creep, damage and lifetime for loaded brittle rocks. Tectonophysics 2014;633:164–75. [21] Cheng YP, Nakata Y, Bolton MD. Discrete element simulation of crushable soil. Geotechnique 2003;53(7):633–41. [22] Cho N, Martin CD, Sego DC. A clumped particle model for rock. Int J Rock Mech Min Sci 2007;44(7):997–1010. [23] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979;29(1):47–65. [24] Dinç Ö, Scholtès L. Discrete analysis of damage and shear banding in argillaceous rocks. Rock Mech Rock Eng 2018;51(5):1521–38. [25] Eberhardt EB. Brittle rock fracture and progressive damage in uniaxial compression. University of Saskatchewan; 1998. [26] El Bied A, Sulem J, Martineau F. Microstructure of shear zones in Fontainebleau

182

Computers and Geotechnics 110 (2019) 161–183

X.F. Li, et al.

[86] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [87] Potyondy DO. The bonded-particle model as a tool for rock mechanics research and application: current trends and future directions. Geosyst Eng 2015;18(1):1–28. [88] Powrie W, Ni Q, Harkness RM, et al. Numerical modelling of plane strain tests on sands using a particulate approach. Géotechnique 2005;55(4):297–306. [89] Quey R, Dawson PR, Barbe F. Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing. Comput Methods Appl Mech Eng 2011;200(17–20):1729–45. [90] Rasmussen LL, de Farias MM, de Assis AP. Extended Rigid Body Spring Network method for the simulation of brittle rocks. Comput Geotech 2018;99:31–41. [91] Schöpfer MPJ, 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. [92] Sheorey PR. Empirical rock failure criteria. AA Balkema; 1997. [93] Shimizu H, Koyama T, Ishida T, et al. Distinct element analysis for Class II behavior of rocks under uniaxial compression. Int J Rock Mech Min Sci 2010;47(2):323–33. [94] Tan X, Konietzky H, Chen W. Numerical simulation of heterogeneous rock using discrete element model based on digital image processing. Rock Mech Rock Eng 2016;49(12):4957–64. [95] Tang CA, Tham LG, Lee PKK, et al. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—part II: constraint, slenderness and size effect. Int J Rock Mech Min Sci 2000;37(4):571–83. [96] Tapponnier P, Brace WF. Development of stress-induced microcracks in Westerly granite. Int J Rock Mech Min Sci Geomech Abstr 1976;13(4):103–12. [97] Vallejos JA, Salinas JM, Delonca A, et al. Calibration and verification of two bonded-particle models for simulation of intact rock behavior. Int J Geomech 2016;17(4):06016030. [98] Villeneuve MC, Diederichs MS, Kaiser PK. Effects of grain scale heterogeneity on rock strength and the chipping process. Int J Geomech 2011;12(6):632–47. [99] Wang X, Cai M. Modeling of brittle rock failure considering inter-and intra-grain contact failures. Comput Geotech 2018;101:224–44. [100] Wang Y, Tonon F. Modeling Lac du Bonnet granite using a discrete element model. Int J Rock Mech Min Sci 2009;46(7):1124–35. [101] Wong TF, Biegel R. Effects of pressure on the micromechanics of faulting in San Marcos gabbro. J Struct Geol 1985;7(6):737–49. [102] Wong TF. Micromechanics of faulting in Westerly granite. Int J Rock Mech Min Sci Geomech Abstr. Pergamon 1982;19(2):49–64. [103] Wu S, Xu X. A study of three intrinsic problems of the classic discrete element method using flat-joint model. Rock Mech Rock Eng 2016;49(5):1813–30. [104] Xia X, Li HB, Li JC, et al. A case study on rock damage prediction and control method for underground tunnels subjected to adjacent excavation blasting. Tunn Undergr Space Technol 2013;35:1–7. [105] Yang ZX, Jardine RJ, Zhu BT, et al. Sand grain crushing and interface shearing during displacement pile installation in sand. Géotechnique 2010;60(6):469. [106] Yu QL, Ranjith PG, Liu HY, et al. A mesostructure-based damage model for thermal cracking analysis and application in granite at elevated temperatures. Rock Mech Rock Eng 2015;48(6):2263–82. [107] Zhang Y, Wong LNY. A review of numerical techniques approaching microstructures of crystalline rocks. Comput Geosci 2018. [108] Zhao T, Crosta GB, Utili S, et al. Investigation of rock fragmentation during rockfalls and rock avalanches via 3-D discrete element analyses. J Geophys Res Earth Surf 2017;122(3):678–95. [109] Zheng Z, Sun WC, Fish J. Micropolar effect on the cataclastic flow and brittleductile transition in high-porosity rocks. J Geophys Res Solid Earth 2016;121(3):1425–40. [110] Zuo J, Li H, Xie H, et al. A nonlinear strength criterion for rock-like materials based on fracture mechanics. Int J Rock Mech Min Sci 2008;45(4):594–9. [111] Bésuelle P, Baud P, Wong TF. Failure mode and spatial distribution of damage in Rothbach sandstone in the brittle-ductile transition. Thermo-Hydro-Mechanical Coupling in Fractured Rock. Basel: Birkhäuser; 2003. p. 851–68. [112] Yu QL, Ranjith PG, Liu HY, Yang TH, Tang SB, Tang CA, Yang SQ. A mesostructure-based damage model for thermal cracking analysis and application in granite at elevated temperatures. Rock Mech Rock Eng 2015;48(6):2263–82.

Geotech Eng 2017;9(4):575–84. [60] Lim SS, Martin CD, Åkesson U. In-situ stress and microcracking in granite cores with depth. Eng Geol 2012;147:1–13. [61] 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 2014;65:96–115. [62] Lisjak A, Liu Q, Zhao Q, et al. Numerical simulation of acoustic emission in brittle rocks by two-dimensional finite-discrete element analysis. Geophys J Int 2013;195(1):423–43. [63] Liu HY, Roquete M, Kou SQ, et al. Characterization of rock heterogeneity and numerical verification. Eng Geol 2004;72(1–2):89–119. [64] Liu Q, Jiang Y, Wu Z, et al. A Voronoi element based-numerical manifold method (VE-NMM) for investigating micro/macro-mechanical properties of intact rocks. Eng Fract Mech 2018;199:71–85. [65] Lockner DA, Byerlee JD, Kuksenko V, et al. Observations of quasistatic fault growth from acoustic emissions. Int Geophys Ser 1992;51:3–13. [66] Lorig LJ, Cundall PA. Modeling of reinforced concrete using the distinct element method. Fracture of concrete and rock. New York, NY: Springer; 1989. p. 276–87. [67] Ma GW, Wang QS, Yi XW, et al. A modified SPH method for dynamic failure simulation of heterogeneous material. Math Probl Eng 2014;2014. [68] Mahabadi OK, Lisjak A, Munjiza A, et al. Y-Geo: new combined finite-discrete element numerical code for geomechanical applications. Int J Geomech 2012;12(6):676–88. [69] Mahabadi OK, Randall NX, Zong Z, et al. A novel approach for micro-scale characterization and modeling of geomaterials incorporating actual material heterogeneity. Geophys Res Lett 2012;39(1). [70] Mahabadi OK, Tatone BSA, Grasselli G. Influence of microscale heterogeneity and microstructure on the tensile behavior of crystalline rocks. J Geophys Res Solid Earth 2014;119(7):5324–41. [71] Mair K, Abe S. 3D numerical simulations of fault gouge evolution during shear: grain size reduction and strain localization. Earth Planet Sci Lett 2008;274(1–2):72–81. [72] Martin CD, Chandler NA. The progressive fracture of Lac du Bonnet granite. Int J Rock Mech Min Sci Geomech Abstr. Pergamon 1994;31(6):643–59. [73] Mas Ivars D, Potyondy DO, Pierce M, et al. The smooth-joint contact model. Proceedings of WCCM8-ECCOMAS, 2008; 2008: 8th. [74] Mayer JM, Stead D. Exploration into the causes of uncertainty in UDEC grain boundary models. Comput Geotech 2017;82:110–23. [75] Mehranpour MH, Kulatilake PHSW. Improvements for the smooth joint contact model of the particle flow code and its applications. Comput Geotech 2017;87:163–77. [76] Moore DE, Lockner DA. The role of microcracking in shear-fracture propagation in granite. J Struct Geol 1995;17(1):95–114. [77] Munjiza AA. The combined finite-discrete element method. John Wiley & Sons; 2004. [78] Munjiza A, Owen DRJ, Bicanic N. A combined finite-discrete element method in transient dynamics of fracturing solids. Eng Comput 1995;12(2):145–74. [79] Nicksiar M, Martin CD. Factors affecting crack initiation in low porosity crystalline rocks. Rock Mech Rock Eng 2014;47(4):1165–81. [80] O'Sullivan C. Particulate discrete element modelling: a geomechanics perspective. CRC Press; 2014. [81] Park JW, Park C, Song JW, et al. Polygonal grain-based distinct element modeling for mechanical behavior of brittle rock. Int J Numer Anal Meth Geomech 2017;41(6):880–98. [82] Paterson MS, Wong T. Experimental rock deformation-the brittle field. Springer Science & Business Media; 2005. [83] Peng J, Wong LNY, Teh CI, et al. Modeling micro-cracking behavior of Bukit Timah granite using grain-based model. Rock Mech Rock Eng 2018;51(1):135–54. [84] Peng J, Wong LNY, Teh CI. Effects of grain size-to-particle size ratio on microcracking behavior using a bonded-particle grain-based model. Int J Rock Mech Min Sci 2017;100:207–17. [85] 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–73.

183