DEM simulation of bonded granular material. Part I: Contact model and application to cemented sand

Computers and Geotechnics 75 (2016) 192–209

Contents lists available at ScienceDirect

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

Research Paper

DEM simulation of bonded granular material. Part I: Contact model and application to cemented sand Zhifu Shen a,b, Mingjing Jiang a,b,⇑, Colin Thornton c a

Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China c School of Chemical Engineering, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK b

a r t i c l e

i n f o

Article history: Received 1 July 2015 Received in revised form 15 February 2016 Accepted 16 February 2016 Available online 27 February 2016 Keywords: Discrete Element Method Cemented sand Bonded contact model Bond geometry Bond failure criterion

a b s t r a c t A three-dimensional bonded contact model is developed in this study to simulate the mechanical behavior of bonded granular material using the Discrete Element Method. Two types of bonded contact are identified, i.e. parallel bond contact where the bond material deposits at a physical contact (two particles being in contact) and serial bond contact where the bond material deposits at a virtual contact (two particles being close but not in contact). The strength and stiffness of the bond material are incorporated into a previous unbonded contact model where the complete interactions in the normal, tangential, rolling and torsional directions are included. A grain is simplified as a sphere and the bond material is idealized in shape as a cylinder with concave ends complementary with the surfaces of the two bonded spheres. The effect of such bond geometry is considered in terms of bond strength and stiffness. The bond failure criterion can capture the coupled effect of normal force, shear force, rolling moment and torque on bond breakage. Parametric studies in DEM simulations show that the bond material properties, the bond content and the bond distribution all influence the behavior of cemented sand. With an increase in bond content or decrease in confining pressure, the mechanical behavior of cemented sand evolves from dilation-dominant where the peak state coincides with the maximum dilation rate, to bond-dominant where the peak state coincides with the maximum bond breakage rate. Microscopic statistics shows that the contributions of shear force and rolling moment to bond failure are statistically equal, both significantly greater than the contribution of torque. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Bonded granular materials are common in natural deposits [1–4] as well as in engineered ground [5,6]. For example, the presence of a small amount of cementation can endow sand with a rich variety of mechanical behavior [2,3,7,8]. As a result, cementation, fabric, density, and stress history together determine the macroscopic behavior and all these factors evolve in a coupled way during deformation. Because of the difficulty in directly observing and quantifying the degradation of cementation, the phenomenological method is often used in constitutive modeling with the microscopic mechanisms unclear [9–14]. With the aid of the Discrete Element Method (DEM), first proposed in [15], bonded granular material can be approximately simulated in an ⇑ Corresponding author at: Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China. Tel./fax: +86 21 6598 0238. E-mail address: [email protected] (M. Jiang). http://dx.doi.org/10.1016/j.compgeo.2016.02.007 0266-352X/Ó 2016 Elsevier Ltd. All rights reserved.

effective way [16–23] and the cementation degradation can be tracked and interpreted quantitatively [19,24,25]. Some key components in a phenomenological constitutive model can be identified, for instance, in the description of the yield surface and plastic potential functions [18] and in the bond degradation law [19]. Basically, bonded granular material is envisioned as an assembly of bonded grains in DEM, where the bond material can be modeled either by explicitly employing relatively fine particles (as bond material) around a contact to link the neighboring grains [24,25], or by employing a bonded contact model. The latter scheme is more efficient than the former one, but a reliable and simple bonded contact model is the key to obtain realistic and efficient DEM simulations. Many bonded contact models have been used so far [16,18–23,26–29]. Two major shortcomings of existing bonded contact models may limit the accuracy of DEM simulations. First, the bond geometry between two grains was not properly considered for simplicity. Usually, a bond is simplified as a beam connecting the centers of two neighboring grains. A more

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

reasonable approximation proposed in [26] idealizes a bond as a short beam connecting the surfaces of two grains, but the beam end is assumed flat which actually should be curved (i.e. complementary with the particle surface). Second, because of the beam assumption, a beam failure criterion and its various modified forms were often used as bond failure criteria but were not fully justified. Experimental research on bond failure criteria has been used to examine some limited loading modes [30,31]. Later, a series of complex two-dimensional (2D) experiments were carried out on bonded aluminum rods under various loading modes [32–34]. These experimental results have been incorporated into bonded contact models for 2D DEM simulations [19,30,35]. Although complete experiments on three-dimensional (3D) bond failure are underway [77], numerical studies have been used recently to try to understand the features in the 3D case [36]. It has been found that the failure criterion can be unified in the same normalized form in both 2D and 3D and that it is necessary to consider the effect of bond geometry [36], which provides the rationale to overcome the two shortcomings of existing bonded contact models. In this paper, a newly developed 3D bonded contact model is introduced to simulate the interactions of two bonded spheres. This model incorporates the complete interactions in the normal, tangential, rolling and torsional directions [37]. There are four new features in this model: (i) Two types of bonded contacts, depending on the geometry of the bond, are identified and modeled separately; (ii) Bond geometry and its effect on stiffness and strength of the bond material are properly considered; (iii) A new complete 3D bond failure criterion proposed in [36] is used to fully capture the coupled effect of normal force, shear force, rolling moment and torque on bond failure; (iv) The post-breakage unbonded contact model is modified in terms of the contact detection threshold and the ability to transmit moments/torques to properly consider the effect of attached post-failure bond material on the particle surface. These features have never been considered fully in previous models but are necessary to establish a physically complete bond contact model. The completeness does not necessarily mean complexity. Some necessary assumptions will be introduced to simplify the complicated 3D bonded contact problem itself so that the proposed model can be easily implemented. The paper is organized as follows. First, the framework of the bonded contact model is introduced in Section 2. Then, the distribution of the bond material inside a sample is discussed in Section 3. The ability of this model is demonstrated by parametric studies of numerical triaxial compression tests on cemented sands in Section 4. Finally, the mechanics of cemented sand are discussed in Section 5 and the microscopic statistics for bond failure modes is presented in Section 6. The bond contact model developed in this paper would be most useful to model porous materials that can be viewed as round/subround particles cemented by bonds, such as cemented sands and

(a) Adapted from [38]

193

some sedimentary rocks. This paper only introduces the model application to low-bond-content granular materials (e.g. cemented sands). In our accompanying paper [76], the model will be extended to a high-bond-content granular material (i.e. methane hydrate bearing sand), where (1) the bonding effect and two other aspects (particle angularity and particle enlargement) are equally important, and (2) the bonding material (methane hydrate) has special mechanical behavior compared with cement in this paper.

2. Framework of the bonded contact model 2.1. Microstructure and idealization of bonded granular materials Fig. 1 presents the typical microstructure of bonded granular materials where grains are cemented by bond material at contacts [38,39]. The amount of cementation can vary greatly at different contacts. Theoretically, the bonded granular material can be considered as an assembly of bonded grains at the particle-scale. Basically, there exist two distinct types of bonded contacts. The first type refers to the case where the bond material forms after the two grains have been in physical contact. For the second type, the bond material forms in the gap between two neighboring grains that are not in physical contact (i.e. virtual contact). Based on the mechanical response of each type of bonded contact, which will be explained below, the first type is named parallel bond contact while the second type is referred to as serial bond contact, which were named zero-thickness bond contact and thick bond contact respectively in [32–34]. Fig. 2 provides the schematic illustrations of idealized microstructure in bonded granular material. The grains are idealized as spheres. For an unbonded contact where there is no bond material or the bond material has been broken in Fig. 2(a), the interaction between the two grains occurs over a contact area since real (sand) grains are irregular in shape. The intact bond material at a bonded contact can be regarded as an elastic continuum in comparison with the overall discontinuous nature of the bonded granular material. For the parallel bond contact in Fig. 2(b), the contact interaction can be considered as two components connected in parallel: one from the physical contact of the two grains which is the same as that illustrated in Fig. 2(a), and the other one from the addition of the bond material. For the serial bond contact in Fig. 2(c), the contact interaction can be considered as two components connected in series: one from the interaction of the two grains since a grain itself is deformable (the bond material is treated as a rigid load transmitter), and the other one from the deformable bond material (the grain boundaries are treated as rigid). The deformability of a grain is modeled through the overlap at a contact which is controlled by the stiffness. The deformation and related stiffness of the bond material will be discussed below.

(b) Adapted from [39]

Fig. 1. Microstructure of cemented sand.

194

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

(a) Unbonded contact

(b) Parallel bond contact

(c) Serial bond contact

Fig. 2. Idealized representation of microstructure in cemented sand.

(a) Parallel bond contact

(b) Serial bond contact Fig. 3. Mechanical components of a bonded contact model (in the online version of this paper, blue: mechanical elements for particle-particle interaction; red: mechanical elements for bond material). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 shows the mechanical components of a bonded contact model represented by elements including springs, divider, slider, roller, twister, and bonds. Interactions in four (normal, tangential, rolling and torsional) directions are considered in the contact model, as shown by the white arrows in the figure. The stiffness for each spring in Fig. 3 is denoted by Ki,j where i = n, s, r and t for normal, tangential, rolling and torsional direction, respectively, and j = p and b for particle and bond material, respectively. The loads (force F and moment/torque M) transmitted by a spring are also denoted in the same way in Fig. 3. For the parallel bond contact in Fig. 3(a), two sets of mechanical elements, one for the particle-particle interaction in physical contact (in blue in the online version of this paper) and the other one for the bond material (in red), are connected in parallel. The spring and bond elements for a bond material in each direction, enclosed by a rectangle in Fig. 3(a), enter into play when the bond material is formed and are disregarded after bond breakage. For the serial bond contact in Fig. 3(b), the springs considering the deformability of a grain (in blue) and a bond material (in red) are connected in series in each direction. The bond element acts as a lock (i.e. the divider, slider, roller and twister are deactivated) so that only an elastic response is allowed at a serial bond contact in the presence

of the bond material. Upon bond breakage, the spring and bond elements for a bond material are disregarded in each direction. Therefore, the bond material is modeled as a brittle elastic medium by its spring and bond elements. After bond breakage, the parallel bond contact and serial bond contact return to an unbonded contact. The springs for a particle and a bond material are simplified as uncoupled in deformation. Otherwise, the bonded contact interaction would become a complicated boundary value problem [40– 50]. The stiffness and bond failure criterion then may not be expressed analytically in a simple way. The contact model in this paper attempts to achieve a balance between capturing the major behavior of a bonded contact in a simple way for DEM simulation and preserving the complexity of the contact problem itself. In this sense, the mechanical systems in Fig. 3 are reasonable since they identify two fundamentally different types of bonded contacts and consider the bond geometry effect on bond stiffness and strength which will be explained below. Based on the framework in Fig. 3, four key aspects of the bonded contact model will be investigated: unbonded contact behavior, bond stiffness, bond failure criterion and the effect of post-failure bond material on contact behavior.

195

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

2.2. Unbonded contact behavior

Table 2 Summary of the behavior of bond material.

An unbonded contact law describes the interaction at a physical contact without bond material. Due to the irregular shapes of real (sand) grains, forces and moments/torques can be transmitted through an area at a contact subjected to relative motion in the four directions. A three-dimensional unbonded contact law has been proposed to consider the effect of contact area by incorporating rolling and torsional resistances [37]. In the normal direction, an elastic behavior is considered, and only a compressive force can be transmitted. In the tangential, rolling and torsional directions, an elastic behavior with a limited strength is assumed in each direction. Details of the contact law have been introduced in [37] and the key formulae are listed in Table 1. When the peak resistance is reached in one direction, the corresponding slider or roller or twister in Fig. 3 is activated. In Table 1, the parameter un denotes inter-particle overlap, Ep is the modulus of the particle material, R is the harmonic mean of the particle radii, n is the ratio of normal stiffness over tangential stiffness, l is the inter-particle friction coefficient, Rc = bcR is the contact radius with bc being a dimensionless parameter linking the contact size with the particle size (see Fig. 2), and fc describes the effect of local asperity crushing (fc = 4 if local crushing is not considered). The parameter Dt denotes the time step used in DEM computations. In Fig. 3 and Table 1, n is a unit vector in the _ r and x _ t are the relcontact normal direction, u_ n (in Table 2), u_ s , x ative normal velocity, tangential velocity, rolling rate and torsional rate at a contact, respectively. 2.3. Elastic behavior of bond material

2.3.1. Normal stiffness Assume that a bond cylinder is subjected to a normal compression of d. The bond cylinder is dissected into hollow cylinder elements for integration operation. Fig. 5(a) shows a typical element with an inner radius of r, a thickness of dr, and a length of h. The bond geometry is shown in Fig. 5(d). The normal force developed in the element is

2prEb ddr ¼ h

2prEb ddr qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hb þ Ds  2 D2s =4  r 2

ð1Þ

F n;b þ K n;b u_ n Dt F s;b þ K s;b u_ s Dt _ r Dt M r;b þ K r;b x _ t Dt M t;b þ K t;b x

Stiffness

Peak strength

Eq. Eq. Eq. Eq.

Eqs. (9a and 9b) Eqs. (11–13)

(3) (4) (7) (8)

Particle boundary

(a) Parallel bond contact

(b) Serial bond contact

Fig. 4. Geometry of a bonded contact.

where Eb is the Yong’s modulus of the bond material. It is assumed that there is no interaction between neighboring elements. Integrating dFn,b in terms of r from 0 to Db/2 leads to

F n;b ¼ Eb d

Db =2

0

2prdr qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hb þ Ds  2 D2s =4  r2

ð2Þ

and the normal stiffness Kn,b can be written as

" # pffiffiffiffiffiffiffiffiffiffiffiffiffiffi! pffiffiffiffiffiffiffiffiffiffiffiffiffiffi F n;b pEb Db n n2  1 n2  1 l ln l  þ  l lnðl  1Þ  1 K n;b ¼ ¼ 2 n n d ð3Þ where m = Hb/Db (bond slenderness ratio) and n = Ds/Db are two dimensionless quantities characterizing the geometry of a bond cylinder, and l = m/n + 1. 2.3.2. Tangential stiffness Analogous to the elastic moduli relationship Gb = 0.5Eb/(1 + mb), the tangential stiffness Ks,b is simply related to the normal stiffness Kn,b by

K s;b ¼

K n;b 2ð1 þ mb Þ

ð4Þ

where Gb is the shear modulus and vb is the Poisson’s ratio of the bond material. Eq. (4) is the same as the analytical expression given in [26]. 2.3.3. Rolling stiffness The rolling stiffness (or bending stiffness in the terminology of beam theory) is derived by applying a rolling rotation hr to the element in Fig. 5(b). The rolling moment (or bending moment) developed in the element can be written as

Table 1 Summary of the formulae for unbonded contacts. Contact response

Stiffness

Peak resistance

Normal Tangential Rolling

F n;p ¼ K n;p un n F s;p F s;p þ K s;p u_ s Dt _ r Dt M r;p M r;p þ K r;p x

K n;p ¼ 2Ep R K s;p ¼ K n;p =n

–

K r;p ¼ 0:25K n;p R2c

0:25fc F n Rc

Torsional

M t;p

K t;p ¼ 0:5K s;p R2c

0:65lF n;p Rc

_ t Dt M t;p þ K t;p x

F n;b F s;b M r;b M t;b

Normal Tangential Rolling Torsional

Z

Since grains are simplified as spheres, the bond material between two spheres can be idealized in shape as a short cylinder (referred to as bond cylinder hereafter) with spherical concave ends, as shown in Fig. 4, where Ds, Db and Hb denote the particle diameter, bond cylinder diameter and the minimum thickness of the bond cylinder, respectively. The Young’s modulus of a bond material is usually lower than that of the particle material, for instance, 10 GPa for Portland cement (bond material) compared to 70 GPa for quartz (particle material). Therefore, it is acceptable to assume that the grains are rigid and a bond cylinder is deformable when dealing with the elastic deformation of bond material. The stiffness of a bond cylinder in each direction (normal, tangential, rolling and torsional) is approximately derived as follows.

dF n;b ¼

Elastic response

lF n;p

dMr;b ¼

pEb hr h 4h

4

ðr þ drÞ  r4

i

ð5Þ

Neglecting terms of higher order than dr leads to

dMr;b ¼

pEb hr h

r 3 dr

ð6Þ

196

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

(a) Normal mode

(b) Rolling mode

(c) Torsional mode

(d) Bond geometry

Fig. 5. Integration elements for deriving stiffness.

and integration gives the rolling stiffness Kr,b ( " # 1:5 Mr;b pEb D3b n3 1 ðn2  1Þ K r;b ¼ ¼ 1 3 hr 16 n3 pffiffiffiffiffiffiffiffiffiffiffiffiffi! m2 þ 2mn n2  1 þ 1 n2 n " pffiffiffiffiffiffiffiffiffiffiffiffiffi!#) 3 2 n2  1 m þ n m þ 3m n þ 2mn2 þ þ lnðl  1Þ  ln l  ð7Þ 2n3 n3 n

2.3.4. Torsional stiffness Similarly, the torsional stiffness Kt,b is derived by applying a torsion ht to the element in Fig. 5(c). The derivation is omitted here but it is found that the following relationship holds,

K t;b ¼

K r;b 1 þ vb

ð8Þ

Rt ¼ St ðRnc þ Rnt ÞDb =2

The parameter Si is a quantity capturing the effect of normal force on strength, 1

mi3

Si ¼ mi1 f n þ mi2 f n ½lnðf n Þ

ð12Þ

where fn = (Fn,b + Rnt)/(Rnc + Rnt) is the dimensionless normal force (i.e. Fn,b normalized by the compression strength Rnc and the tension strength Rnt), and mi1, mi2, mi3 (i = s, r, t) are fitting parameters. The strength envelopes are shown in Fig. 6. Note that Eq. (11) gives the bond strength if only one of the shear, rolling and torsional loads is acting on the bond cylinder in addition to the normal force. In the general case, the shear force Fs,b, rolling moment Mr,b and torque Mt,b together bring a bond to breakage in a coupled way. A coupled failure criterion can be practically described by

 2

2

2

fs þ fr þ ft

In summary, the stiffness of a bond cylinder in each direction is a function of the elastic constants of the bond material (Eb and vb) and the geometry of a bond cylinder (Hb, Db and Ds).

ð11cÞ

<1

intact

P 1 broken

ð13Þ

where fs = Fs,b/Rs, fr = Mr,b/Rr, and ft = Mt,b/Rt. This coupled failure criterion is visualized in Fig. 6. Table 2 summarizes the formulae to describe the bond behavior.

2.4. Bond failure criterion The numerically obtained bond failure criterion in [36] can consider the bond geometry effect on bond strength and can capture the coupled effect of normal force, shear force, rolling moment, and torque on bond failure. This bond failure criterion is used in this paper and is introduced briefly below. The criterion includes five basic failure modes (compression, tension, shear, rolling and torsion) and a coupled failure mode which includes the five basic modes as special cases. In compression failure, a bond is crushed by a compressive force that exceeds the ultimate compressive force Rnc. In tension failure, a bond is pulled apart by a tensile force that exceeds the ultimate tensile force Rnt. The two ultimate forces can be expressed as

Rnc ¼ vc rc pðDb =2Þ2 ;

Rnt ¼ vt rt pðDb =2Þ2

2.5. Effects of post-failure bond material on contact behavior The contact behavior after bond breakage is considered to be that of an unbonded contact. Experiments on bonded particles [32–34,77] and numerical simulations [36] show that the bond material at a serial bond contact is still effective to transmit normal force at the instant of shear, rolling (and torsional in 3D) failures. Only under pure compression will the bond be crushed into

ð9a; bÞ

where rc and rt are the compressive and tensile strengths of the bond material. The two factors vc and vt account for the bond geometry effect, which, according to [36], can be correlated to the two dimensionless geometry factors, m = Hb/Db and n = Ds/Db, as

vc ¼ C 1 eC2 m nC3 ; vt ¼ 1

ð10a; bÞ

where C1, C2 and C3 are three fitting parameters. A bond may fail in shear, rolling or torsion mode under a normal force Fn,b with Rnt < Fn,b < Rnc. The corresponding strength Ri (i = s in shear failure, i = r in rolling failure, and i = t in torsion failure) depends on the compression strength Rnc, the tension strength Rnt and the bond diameter Db,

Rs ¼ Ss ðRnc þ Rnt Þ

ð11aÞ

Rr ¼ Sr ðRnc þ Rnt ÞDb =2

ð11bÞ

Fig. 6. Failure envelopes of bond material.

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

fragments and detached from the contact. Fig. 7 presents a conceptual diagram of the physical transition process of a bonded contact to a post-breakage unbonded contact. To ensure the physical correctness, two modifications to the unbonded contact model are introduced to describe this transition as follows. 2.5.1. Modification 1 For a serial bond contact in Fig. 7(a), the bond material forms at a negative overlap (un = Hb, also see Fig. 4) and the bonded contact law in the normal direction can be described by the solid Line 1 in Fig. 7(a). After bond breakage, the bond material can bridge the gap between the two previously bonded particles so that they can interact as if they were enlarged in size, as illustrated by the enlarged particle boundary in Fig. 7(a). A new threshold value un0 (=0 by default) has to be used in the post-breakage unbonded contact model so that the normal contact law, represented by the dashed Line 2 in Fig. 7(a), becomes:

F n;p ¼ K n;p ðun  un0 Þ

ð14Þ

The value of un0 is defined so that the normal force does not change due to bond breakage. That is, the normal force transmitted through a serial bond contact is totally carried now by the post-breakage unbonded contact. If the normal force Fn,b at the instant of bond breakage is tensile, un0 is set equal to Hb and the post-breakage unbonded contact law can be described by the dashed Line 3 in Fig. 7(a). In this case, a jump in the normal force upon bond breakage is observed because the post-breakage unbonded contact cannot resist tensile force. The shear force, rolling moment and torque previously transmitted by the bond material are now carried by the post-breakage unbonded contact. If the modification of un0 is not applied (as in previous bonded contact models in the literature), bond breakage could reduce a serial bond contact to a virtual

(a) Normal direction in a serial bond contact

197

contact that cannot transmit any load, leading to unrealistic instability at local neighboring contacts. The modification of un0 is not used for a serial bond contact if fn is greater than a specific value (e.g. 0.9) because then the bond material is crushed by the compressive force. The modification of un0 is not used in a parallel bond contact in our model because the physical contact of the two grains keeps transmitting load after bond breakage as long as un > 0. That is, the loads (normal force, shear force, rolling moment and torque) transmitted by the bond element are removed immediately upon bond breakage in a parallel bond contact. As a result, the failure of both types of bonded contacts may be accompanied by plasticlike sliding/rolling/twisting or contact splitting (tensile failure), local adjustment of structure, and seismic emission if the adjustment is strong enough. These failure modes then correspond well with the shear, tensile and mixed microscopic failure types revealed by moment tensor analysis of acoustic emission [78].

2.5.2. Modification 2 The ability of an unbonded contact model to transmit shear force, rolling moment and torque depends on the contact size Rc defined as bcR (R: harmonic mean of the particle radii, see Fig. 2 and Table 1). For a post-breakage unbonded contact, if the size of the bond cylinder Rb(= Db/2) is greater than bcR, it is reasonable to assume that Rc equals Rb at the instant of bond breakage, i.e. state ‘‘A” in Fig. 7(c). With the increase of tangential displacement S on the contact plane, the post-failure bond material can gradually get removed due to relative motion at a contact, rendering Rc monotonically reducing from Rb to bcR, i.e. state ‘‘B” in Fig. 7(c). If S reverses (decreases), Rc remains constant. If Rb is less than bcR, there is no need to consider this modification. This modification applies to both serial and parallel bond contacts.

(b) Normal direction in a parallel bond contact

(c) Evolution of contact size after bond breakage Fig. 7. Illustration of effects of post-failure bond material on contact behavior.

198

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

Our simulations (to be introduced later in the paper) indicate that the two modifications are trivial for cemented sand with very low bond content (e.g. bond material volume being several-percent of sand particle volume). However, they greatly improve model performance for sand with high bond content, such as the methane hydrate bearing sand discussed in our accompanying paper [76], which can be treated as cemented sand and where the bond material (methane hydrate) can fill up to half of the pore space [51,52]. As the macroscopic loading continues, a sample gradually loses these post-breakage unbonded contacts with the above mentioned modifications and gains newly-generated unbonded contacts without modification. This is reasonable for relatively low bond content samples since the bond material is considered to be completely detached from the grain surfaces and exists only in the pores, having no influence on the mechanical behavior. For those geomaterials with a high bond content, such as methane hydrate bearing sediment, post-failure bond material may not be completely removed from particle surface, thus influencing particle shape and effective particle size. Then, the threshold value un0 to detect new contact may need to be further modified and the shape parameter bc may need to be increased in a newly-generated unbonded contact. These issues are discussed in our accompanying paper [76]. In summary of Section 2, a combination of the formulae in each direction from Table 1 and Table 2 in either a ‘‘parallel” or a ‘‘serial” way as outlined in Fig. 3 leads to the bonded contact model proposed in this paper. 3. Distribution of bond material In Fig. 4, two quantities (Hb and Db) are used to characterize the geometry of a bond cylinder. Due to the variation of particle size and the non-homogeneous deposition of bond material, the two quantities may vary greatly at different contacts. To describe the bond material distribution and its variability in a sample, two microscopic properties of low-bond-content cemented sand are discussed preliminarily as follows. The special way of bond material (methane hydrate) distribution in methane hydrate bearing sand is discussed in our accompanying paper [76]. 3.1. Connectivity Bond material can be considered as an agent to connect neighboring grains for load transmission. Connectivity is here defined as a microscopic property of cemented sand to reflect the ability of

the bond material to establish an effective cohesive mechanical connection at both physical and virtual contacts. Fig. 1 indicates that the bond slenderness ratio (m = Hb/Db) at a contact is limited to a certain range. A slim bond with large slenderness ratio is rarely observed. This is possibly related to the depositing behavior of the bond material. Besides, a slim bond, even if it does exist occasionally, is mechanically fragile, thus contributing little to the macroscopic strength of cemented sand. Therefore, slenderness ratio m should be limited in a range of [0, m0] with m0 being the maximum possible slenderness ratio, which can be chosen as a measure of mechanical connectivity at the microscopic scale.

3.2. Variability A new dimensionless parameter bb is introduced to link the diameter of the bond cylinder Db and the diameter of the particle size Ds (=2R, see Fig. 2),

Db ¼ bb Ds

ð15Þ

The amount of cementation exhibits great variation at the different contacts, as shown in Fig. 1. It is assumed here that bb follows a simple normal distribution to capture such internal variability, 9ðb b Þ2 3  b 0 f ðbb ; b0 ; wÞ ¼ pffiffiffiffiffiffiffi e 2w2 2pw

ð16Þ

where b0 is the mean value of bb and w (equal to three times the standard deviation) is the half-width of the distribution, as shown in Fig. 8. The values of bb outside the range [b0  w, b0 + w] are not used. The above discussed connectivity and variability are considered when assigning the bonded contact model to a candidate contact (including virtual contacts) with known Hb. A random bb is first generated from the probability distribution in Fig. 8. The slenderness ratio m = Hb/Db then can be calculated and checked. If m is not in the range [0, m0], this contact is skipped. Otherwise, the candidate contact is assigned with a parallel bond contact model if Hb = 0 or a serial bond contact model if Hb > 0. Repeating the same procedure at all candidate contacts leads to a cemented sample and the total bond volume can be calculated. By adjusting the bb distribution (b0 and w) or m0, samples with different bond volume contents (Vb/Vs, Vb: volume of all bond material, Vs: volume of all particles) can be obtained.

Probability density f (βb ,β0 ,w)

4. Parametric study using DEM simulations

w

β0

βb Fig. 8. Normal distribution of bb.

The new bonded contact model was implemented using PFC3D [53] and applied to simulate the mechanical behavior of cemented sand. Fig. 9 shows the numerical assembly and its particle size distribution. Each numerical assembly consists of 20,000 spheres with dimensions of 4.18 mm  4.18 mm  8.36 mm, confined by six frictionless boundary walls. The assembly was colored in layers to observe deformation. Two assemblies with initial void ratios of 0.82 and 0.91 (under zero stress) were used, representing medium-dense and loose samples. The sample was first isotropically compressed to 50 kPa using the unbonded contact model. The bonded contact model was then implemented at physical contacts and some virtual contacts following the procedure outlined in Section 3. The sample was then isotropically compressed to different target confining pressures for triaxial compression tests. The motion of the boundary walls and the forces acting on them were measured to calculate the strains and stresses.

199

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

100

Percentage smaller than by weight (%)

80

60

40

20

0 0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

Particle diameter (mm) Fig. 9. Numerical assembly and particle size distribution.

4.1. Model parameters

4.2. Effect of bond content

The parameters for the unbonded contact model were calibrated to reproduce the common mechanical behavior of clean sand. These parameters are Ep = 0.7 GPa, n = 5, l = 0.5, bc = 0.15, fc = 4, and are the same for all simulations. The parameters related to bond have clear physical meanings so that they are possible to be determined directly as follows. First, the 12 fitting parameters are phenomenological and can be simply obtained from experiments on 3D bonded spheres [77] as in 2D [32–34]. However, no complete experiment has yet been carried out on 3D bond failure. These parameters have to come from the numerical simulations [36], and they are believed to be representative since the 3D numerical results are qualitatively consistent with the 2D experiments with proper normalization in comparable aspects. The 12 fitting parameters are C1 = 1.444, C2 = 0.702, C3 = 0.168, ms1 = 0.055, ms2 = 0.439, ms3 = 0.696, mr1 = 0, mr2 = 0.324, mr3 = 0.757, mt1 = 0.031, mt2 = 0.327, and mt3 = 0.692. The strength envelopes described by these fitting parameters are shown in Fig. 6. Second, a common Poisson’s ratio of the bond material vb = 0.2 was chosen. The common strengths of bond materials, i.e. rc = 50 MPa, 80 MPa and rt = 5 MPa, 7 MPa, was used. Eb/Ep (=0.05 and 0.2 in this paper) can be evaluated by the Young’s modulus of the bond material and sand particle material. Third, the parameter b0 was chosen according to bond volume content. The parameters related to bond distributions (w, m0) can be estimated by observing the microscopic structure of cemented granulates. However, this is difficult due to the lack of detailed microscopic structure observation at present in the literature. A calibration process is then indeed needed if experimental data are to be matched. However, because of the simplifications in particle shape and contact model, it is not expected that the DEM could give precisely matching simulation results even when these parameters can be reliably evaluated. Due to the shared bond effect and granular nature between real cemented granulate and the ‘‘DEM material”, it is more useful to apply DEM that is enriched by physically complete bond contact model to investigate the general mechanical behavior of cemented granulate, than to purely match experimental data. Following this idea, the parameters in this paper were chosen from common ranges to show the ability of the model. The basic set of parameters are Eb/Ep = 0.05, rc = 80 MPa, rt = 7 MPa, m0 = 0.1, w = 0.2. They will be subjected to change for parametric study. One parameter was changed when its effect was studied while the others remained the same as in the basic set of parameters.

The sample with an initial void ratio of 0.82 (medium-dense) was used to study the effect of bond content on the behavior of cemented sand at a confining pressure of 100 kPa. The value of b0 was chosen so that samples with bond volume contents of 0.5% (b0 = 0.235), 1% (b0 = 0.29), 2% (b0 = 0.348) and 3% (b0 = 0.387) were formed. Fig. 10 presents the numerical simulation results. The results obtained from an uncemented sample are also shown for reference. The bond breakage is the ratio of broken bond number over initial bond number. Fig. 10(a) shows that the presence of cementation can greatly increase the initial stiffness (slope) and peak strength of a sample. After the peak state, the cemented sand experiences softening with the residual strength tending to converge with the strength of uncemented sand. In Fig. 10(b), the percentage of broken bonds increases with axial strain first at an increasing rate and then at a decreasing rate. Not all the bonded contacts are broken at an axial strain beyond 30% and a higher percentage of bonds can survive if the bond content is increased. Note that the compression and tension strengths (Rnc, Rnt) and therefore the shear, rolling and torsional strengths (Rs, Rr, Rt) vary at different contacts due to their different slenderness ratio m and bb. Some contacts with higher strengths or in more favorable configurations (i.e. position and contact normal direction) have a higher chance to avoid breakage. Particle clusters are formed from these surviving bonded contacts. Fig. 11 shows the random spatial distributions of clusters which are distinguished by colors (see the online version of this paper) and grouped by N, the number of member particles in a cluster. These particle clusters are characterized by highly irregular shapes compared with an individual spherical particle. The irregular particle clusters enable a sample to reach a void ratio that is impossible for an uncemented sample, which is clearly shown by the increased dilation in Fig. 10(c). Particle clusters also slightly increase the steady-state strength, which has been observed in experiments [54–57]. The survival of particle clusters at large strain has been numerically captured in [24,25] and experimentally confirmed in [55]. In [55], cemented sands were immersed in deaired water after triaxial compression tests and were found to form a ‘‘cohesionless heap”with grain clusters and blocks. The effect of particle clusters will be further explored in Section 5. In Fig. 10(a), the increase in bond volume content leads to an increase in the initial stiffness and the peak strength. This is because both the number of bonded contacts and the strength of each bonded contact increase with the bond content. Since less bonded contacts are broken with a higher bond content in Fig. 10 (b), the residual strength, dilation rate and steady-state volumetric

200

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

1200

100

σ1−σ3 (kPa)

1000 800

Bond breakage (%)

Uncemented 0.5% 1% 2% 3%

600 400

80 60 40 30

40

20 10

20

200

0 0

0

0

5

10

15

20

25

30

0

35

0

5

10

15

2

20

3

4

25

5

30

35

ε1 (%)

ε1 (%)

(a) Stress

1

(b) Bond breakage evolution

strain response -12 -10

εv (%)

-8 -6 -4 -2 0 2

0

5

10

15

20

25

30

35

ε1 (%)

(c) Volumetric response Fig. 10. Effect of bond content on the behavior of cemented sand (e0 = 0.82, r3 = 100 kPa).

N = 45

291

N = 15

N = 3, 4

44

N=9

N=2

14

N=5

8

N=1

Fig. 11. Spatial distributions of particle clusters (e0 = 0.82, r3 = 100 kPa, bond content = 1%, e1 = 15%, N: number of member particles in a cluster.).

201

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

underestimate the peak strength of a cemented sand. The new model in this paper can implement the bond effect at a virtual contact (i.e. between two separate particles) as well. This is reasonable since in reality bond materials can bridge the gap between two neighboring grains, as shown in Fig. 1. The parameter m0 (bond slenderness ratio limit) in the new model controls the possible maximum bond slenderness ratio. Fig. 14 presents the numerical simulation results for samples with the same bond volume content (1%) but different values of m0 (0, 0.05, 0.1). The value of b0 increases with the decrease of m0 because of the constraint of equal bond volume content, i.e. b0 = 0.29, 0.307, 0.368 for m0 = 0.1, 0.05, 0.0, respectively. A sample with greater m0 contains more bonded (virtual) contacts but each bond is generally weaker (due to smaller b0) than that in a sample with lower m0. Fig. 14 shows that an increase in m0 leads to an increase in the initial stiffness, peak strength and dilation. The simulation results indicate that ignoring bond effects at virtual contacts in previous models in the literature underestimates the effects of bond material on the behavior of a cemented sand.

strain all increase with the bond content. These observations are consistent with the general trend in experiments [56,58–61]. 4.3. Effect of bond material modulus The sample with an initial void ratio of 0.82 (medium-dense) and a bond volume content of 1% was used to study the effect of bond material modulus on the behavior of cemented sand at a confining pressure of 100 kPa. Fig. 12 presents the numerical simulation results obtained for two values of bond material modulus, i.e. Eb/Ep = 0.05 and 0.2. Fig. 12(a) shows that the sample stiffness and post-peak softening rate increase with Eb but the peak strength is the same. With a higher Eb, less inter-particle relative motion is required to reach bond failure at a contact. Therefore, in Fig. 12(b), an increase in Eb leads to an increase in the bond breakage rate and thus an increase in the post-peak softening rate in the stress-strain curve. The steady-state stress and volumetric dilation are lower in the sample with a higher Eb because more bonds are broken. This parametric study indicates that bond material with a relatively low modulus can improve the ductility of cemented sand because bonds are broken at a relatively low rate.

4.6. Effect of bond distribution variability

4.4. Effect of bond material strength

One major internal variable of a cemented sand comes from the non-uniform deposition of bond material at different contacts. In DEM simulations, since bb is randomly chosen from [b0  w, b0 + w] for each candidate contact (see Fig. 8), w can therefore reflect such internal variability of a cemented sand. Fig. 15 presents the effect of w on the mechanical behavior of a cemented sand with 1% bond volume content. Due to the constraint of equal bond volume content, the value of b0 increases slightly with the decrease of w, i.e. b0 = 0.29, 0.305, 0.31 for w = 0.2, 0.1, 0.0, respectively. Fig. 15 (a) shows that an increase in w leads to a decrease in the peak strength but an increase in the post-peak ductility. Fig. 15(b) shows that the effect of w on volumetric dilation and bond breakage is minimal. Fig. 15(c) presents the probability distributions of bb for intact bonded contacts at different levels of axial strain in a sample with w = 0.2. It is clear that the minimum and mean values of bb increase with axial strain. That is, relatively weaker bonds (with lower bb) tend to break earlier, which can explain the effect of w on peak strength and ductility as follows. The mobilization of peak strength in a cemented sand is accompanied by bond breakage. For the sample with w = 0.2 and b0 = 0.29, the wide distribution of bb (from 0.09 to 0.49) makes the peak strength depend largely on the relatively weak bonded contacts (with low bb). On the other hand, for

The effect of bond material strength (rt and rc) was investigated using a medium-dense sample (e0 = 0.82) with 1% bond volume content and r3 = 100 kPa. Two tensile strengths (5 MPa and 7 MPa) and two compressive strengths (50 MPa and 80 MPa) were used. These strength values are common for typical bond materials. Fig. 13 presents the numerical simulation results. It is expected that an increase in bond strength makes it difficult to break bonded contacts. As a result, in Fig. 13(a), the peak strength and volumetric dilation increase while the bond breakage decreases slightly with the increase in tensile strength. The same trend is observed when the compressive strength is increased in Fig. 13(b). Note that the effect of the tensile strength is more significant than the compressive strength, which will be explained in Section 6. This parametric study may provide theoretical support for the use of fiber in engineered cemented sand because fiber is believed to be effective in increasing the tensile strength of bond material [62–65]. 4.5. Effect of bond connectivity In most existing bonded contact models in the literature, bonds were considered only at physical contacts, which may

-8

600

1%, Eb/Ep = 0.2

400

-6 -4 -2 0

300

Bond breakage (%)

σ1−σ3 (kPa)

500

εv (%)

Uncemented 1%, Eb/Ep = 0.05

200 100 0

0

5

10

15

20

25

30

35

2 80 60 40 20 0

0

5

10

15

20

ε1 (%)

ε1 (%)

(a)

(b)

25

30

35

Fig. 12. Effect of bond material modulus (Eb/Ep) on the behavior of cemented sand (e0 = 0.82, r3 = 100 kPa, bond content = 1%): (a) stress–strain response, and (b) volumetric response and bond breakage evolution.

202

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

-8 Uncemented 1%, σt = 7MPa

400

1%, σt = 5MPa

-6

εv (%)

500

-4 -2 0

300 200 100 0

0

5

10

15

20

25

30

35

Bond breakage (%)

σ1−σ3 (kPa)

600

2 80 60 40 20 0 0

5

10

ε1 (%)

20

25

30

35

25

30

35

ε1 (%)

(a) Effect of tensile strength

t

-8

600

-6

εv (%)

Uncemented 1%, σc = 80MPa

500

1%, σc = 50MPa

400

-4 -2 0

300 200 100 0

0

5

10

15

20

25

30

35

Bond breakage (%)

σ1−σ3 (kPa)

15

2 80 60 40 20 0 0

5

10

ε1 (%)

15

20

ε1 (%)

(b) Effect of compressive strength

c

Fig. 13. Effect of bond material strength (rt and rc) on the behavior of cemented sand (e0 = 0.82, r3 = 100 kPa, bond content = 1%).

-8

600

-6

εv (%)

1%, m0= 0.05

400

1%, m0= 0.0

200

0

0

5

10

15

20

25

30

-4 -2 0

35

Bond breakage (%)

σ1−σ3 (kPa)

Uncemented 1%, m0= 0.1

2 80 60 40 20 0 0

5

10

15

20

ε1 (%)

ε1 (%)

(a)

(b)

25

30

35

Fig. 14. Effect of bond slenderness ratio limit (m0) on the behavior of cemented sand (e0 = 0.82, r3 = 100 kPa, bond content = 1%): (a) stressstrain response, and (b) volumetric response and bond breakage evolution.

the sample with w = 0.0 and b0 = 0.31, bb is the same for all bonded contacts and the peak strength is mobilized with the breakage of bonded contacts mainly in unfavorable configurations (i.e. position and contact normal direction). Obviously, the bond strength associated with the peak strength in the former case is generally smaller than in the latter one, and therefore the peak strength decreases with the increase of w. After the peak state, on the contrary, the remaining bonded contacts in the sample with a larger w are generally stronger that those in the sample with a smaller w, thus presumably the former case being slightly more effective in maintaining particle clusters which lead to relatively lower softening rate. This parametric study implies that a certain level of internal variability can increase the post-peak ductility but it is accompanied by a reduction in the peak strength.

In existing research on cemented sands, the type of bond material and its content are the major concerns while the effect of bond distribution is not fully understood. However, in both naturally occurring and engineered cemented sands, the property, content and distribution of the bond material are subjected to spatial and temporal variations. The microstructure is also found to be dependent on the type of bond material [66]. It is therefore fundamentally important to obtain a complete understanding of the effects of bond material property (strength and modulus), bond content and microscopic distribution (connectivity and variability) on the behavior of cemented sand. As shown in the above parametric studies, the new model in this paper can consider these factors in a consistent way, and allow to quantitatively evaluate the effect of each factor, which is difficult to investigate experimentally.

203

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

800

-8 -6 -4 -2 0

400

200

0

0

5

10

15

20

25

30

35

Bond breakage (%)

σ1−σ3 (kPa)

600

εv (%)

Uncemented 1%, w = 0.2 1%, w = 0.1 1%, w = 0.0

2 80 60 40 20 0 0

5

10

15

ε1 (%)

20

25

30

35

ε1 (%)

(a)

(b) 7

ε1 =3.3%

Probability density

6 5

ε1 =0%

4

Eq. (16)

}

ε1 =9.7%

ε1 =32%

3 2 1

0 0.0

0.49

0.29

0.09 0.1

0.2

0.3

0.4

0.5

0.6

βb

(c) Fig. 15. Effect of bond distribution variability (w) on the behavior of cemented sand (e0 = 0.82, r3 = 100 kPa, bond content = 1%): (a) stress–strain response, (b) volumetric response and bond breakage evolution, and (c) probability densities of bb at different axial strains (b0 = 0.29, w = 0.2).

5. Strength, bond breakage and dilation of cemented sand 5.1. Strength envelopes Triaxial compression tests under confining pressures from 50 kPa to 400 kPa and unconfined compression tests were simulated on samples with initial void ratios of 0.91 and 0.82. Four bond volume contents are considered, i.e. 0.5%, 1%, 2% and 3%. The corresponding values of b0 are 0.235, 0.29, 0.348, 0.387 for e0 = 0.82 and 0.244, 0.3, 0.36, 0.4 for e0 = 0.91. The basic set of parameters was used. Fig. 16 shows the effect of confining pressure on the mechanical behavior of cemented sand in triaxial compression tests. Fig. 16 shows that the post-peak softening rate and volumetric dilation are suppressed by increasing confining pressure. In Fig. 16(a), softening and dilation are always observed under different confining pressures since the host sand is medium-dense with e0 = 0.82. For the loose cemented samples with e0 = 0.91, Fig. 16(b) shows that if the confining pressure does not exceed 100 kPa, volumetric dilation and post-peak softening are observed; if the confining pressure reaches 400 kPa, the sample experiences strainhardening and contraction, similar to the loose host sand; at a confining pressure of 200 kPa, the mechanical behavior exhibits a transition characterized by strain-softening and contraction. Similar effects of confining pressure have been observed in experiments [55,56,58–61,67–70]. It is interesting to note that, for the sample with e0 = 0.91 and r3 = 100 kPa, the incremental volumetric strain becomes contractive beyond an axial strain of 20%, following the previous dilation resulted from irregular particle clusters. Similar behavior was observed in [71]. This is because

the loose structure produced by dilation becomes unstable and has to collapse toward a denser state when the clusters are largely degraded with the breakage of bonded contacts. In unconfined compression tests, the cemented samples were formed at a confining pressure of 50 kPa, and then unloaded to zero stress state before axial loading. Fig. 17 shows that the compressive stress first increases linearly with the axial strain. After yielding, the stress increases at a decaying rate until a peak is reached. Both the initial stiffness and peak strength increase with bond content and sample density. The strain required to reach the peak state also increases with bond content. Fig. 18 presents the peak strength envelopes in p–q plane where p = (r1 + 2r3)/3, q = r1  r3. The critical state strength of the host sand is also plotted with gcs = (q/p)cs = 1.15, i.e. ucs = 28.8°. For the sample with a bond content of 0.5% and e0 = 0.91, the peak strength lies above the host sand strength if r3 < 200 kPa and almost merges with the host sand strength when r3 increases up to 200 kPa. The same merging trend is observed for other samples, indicating that the bond effect is less significant under higher confining pressure. Fig. 19(a) shows that the experimentally obtained peak strength envelopes of cemented sands are curved [7,54,55,60,69,72–75]; that is, the strength envelope of a cemented sand first diverges from and then tends to merge with the envelope of the host sand. Our DEM simulation results in Fig. 18 only capture the merging part of a strength envelope. This is presumably because the cemented sample has no ‘‘defect”. A certain level of defect can be created if some of the eligible candidate contacts, which are chosen randomly and not spatially correlated, are not assigned with the bonded contact model. Fig. 19(b) shows that a 30–50% material

204

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

εv (%)

1000

σ3 = 400 kPa

600 400

σ3 = 200 kPa

200 0

Bond breakage (%)

σ1−σ3 (kPa)

800

σ3 = 100 kPa σ3 = 50 kPa 0

5

10

15

20

25

30

ε1 (%)

80 60 40 20 0 0

35

10

15

20

25

30

35

25

30

35

ε1 (%)

-6 -4

εv (%)

σ3 = 400 kPa 600

σ3 = 200 kPa σ3 = 100 kPa 200

σ3 = 50 kPa 0

5

10

15

20

25

30

ε1 (%)

-2 0 2

400

35

Bond breakage (%)

σ1−σ3 (kPa)

5

(a) e0 = 0.82

800

0

-12 -10 -8 -6 -4 -2 0 2

4 80 60 40 20 0 0

5

10

15

20

ε1 (%)

(b) e0 = 0.91

Fig. 16. Mechanical responses of cemented sands with different void ratios in triaxial compression tests (bond content = 1%).

1000

800

0.5% 1% 2% 3%

600 400

400

200

200 0

0.5% 1% 2% 3%

600

σ1 (kPa)

σ1 (kPa)

800

0

2

4

6

8

10

0

0

2

4

6

ε1 (%)

ε1 (%)

(a) e0 = 0.82

(b) e0 = 0.91

8

10

Fig. 17. Mechanical responses of cemented sands with different void ratios in unconfined compression tests (bond content = 0.5%, 1%, 2%, 3%).

defect enables a curved strength envelope to be captured in DEM simulations. Such defects are ubiquitous in naturally cemented sand possibly due to non-uniform deposition and in artificially cemented sand possibly due to poor mixing. However, defects should not change the general trend of cemented sand behavior. The current work is not meant to match experimental results but to show the ability of the proposed model. To obtain a better validation, experimental works must be carried out to obtain not only macroscopic responses but also microscopic information of a bonded assembly. Such works are challenging but are under way in the authors’ group, such as 3D bond contact behavior, as reported in [77].

5.2. Correlations between strength, bond breakage and dilation Fig. 20 presents the evolutions of stress, bond breakage rate and dilation rate with axial strain for samples with different bond contents (e0 = 0.82, r3 = 100 kPa). The bond breakage rate is defined as the tangent slope of the ‘‘bond breakage-axial strain” curve in Fig. 10. The dilation rate is defined using the total strain. For the two samples with bond contents of 2% and 3%, the peak strength coincides with the maximum bond breakage rate, prior to the maximum dilation rate. This observation implies that the peak strength is controlled by bonded contacts. For the sample with a bond content of 0.5%, the stress and dilation rate maintain almost constant

205

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

1600

1200

Tension cutoff

800 600 Uncemented 0.5% 1% 2% 3%

400 200 0

0

200

Tension cutoff

1400

q = σ1−σ3 (kPa)

q = σ1−σ3 (kPa)

1000

400

600

1200 1000 800 600

Uncemented 0.5% 1% 2% 3%

400 200 800

p = (σ1+2σ3)/3 (kPa)

0

0

200

400

600

800

1000

p = (σ1+2σ3)/3 (kPa)

(a) e0 = 0.91

(b) e0 = 0.82

Fig. 18. Simulated strength envelopes of cemented sands with different void ratios (bond content = 0.5%, 1%, 2%, 3%).

4000

Tension cutoff

q = σ1−σ3 (kPa)

3000

2000

Uncemented 1.5% 3.0% 4.5% 6.0% 9.0%

1000

0

0

400

800

1200

1600

2000

p = (σ1+2σ3)/3 (kPa)

(a) 1000

Tension cutoff

q = σ1−σ3 (kPa)

800

600

400

Uncemented no defect 30% defect 50% defect

200

0

0

200

400

600

800

1000

p = (σ1+2σ3)/3 (kPa)

(b) Fig. 19. Strength envelopes of cemented sands: (a) experiment data [54] (bond content = 1.5%, 3%, 4.5%, 6%, 9%), and (b) DEM simulations with different degrees of defect (e0 = 0.82, bond content = 1%).

where the post-peak stress–dilatancy curves are inaccurate due to the formation of shear bands and difficulty in measuring in-band deformations. Fig. 21 presents the stress–dilatancy relationship in our triaxial compression simulations with macroscopic uniform deformation, where des = 2(e1  e3)/3. Fig. 21 shows that the stress ratio g increases with D nonlinearly at a decreasing rate before the peak state. The dilation curves present a hook at the maximum dilation rate and then D decreases gradually to zero. This pattern of stress–dilatancy behavior is qualitatively consistent with and should be more reliable (due to the macroscopic uniform deformation) than that in previous work. Fig. 21 shows that under confining pressures of 50 kPa and 100 kPa, the peak strength occurs prior to the maximum dilation rate. When the confining pressure increases up to 200 kPa, the peak strength coincides with the maximum dilation rate. Therefore, an increase in confining pressure drives the behavior transition from bond-dominant to dilation-dominant. Though dilation appears to be an independent factor in the above correlations, bond breakage plays a critical role behind volumetric change, which is explored in Fig. 22. When the remaining bonds in a sample were set unbreakable at a certain axial strain (marked with a circle in Fig. 22), further loading causes the deviator stress and volume dilation to increase to a higher value than the original one. If bond breakage is allowed, less effort is generally required to overcome interlocking of particle clusters by breaking bonded contacts than by overriding a cluster (which causes dilation). Therefore, bond breakage is the core factor behind the mechanical behavior of cemented sand, either bond-dominant or dilation-dominant. A correlation between the bond breakage rate and the plastic strain has been established based on 2D DEM simulations and is found promising as a bond degradation law in a phenomenological constitutive model [19]. The same correlation will be investigated based on 3D DEM simulations using this new bonded contact model in future work.

6. Statistics for bond failure conditions in the range of axial strain from 1.4% to 3.3%, and the peak strength coincides with the maximum dilation rate, later than the maximum bond breakage rate at an axial strain of 0.77%. This observation implies that the peak strength is controlled by dilation. The behavior of the sample with an intermediate bond content of 1% represents a transition from bond-dominant to dilationdominant. Similar bond-dominant behavior has been observed in experiments [25,61] and numerical simulations [24,25]. Stress–dilatancy behavior has been previously reported in experiments [7,55,68,72,75] and numerical simulations [25],

Although bond breakage is found important in the mechanics of cemented sand, little is known about how a bond is broken. The bond failure criterion in this paper provides full access to the conditions when a bond fails, i.e. the combination of normal force, shear force, rolling moment and torque. The first quantity examined is the bond normal force Fn,b at the instant of breakage, which can be normalized as fn = (Fn,b + Rnt)/(Rnc + Rnt). Fig. 23 presents the frequency distribution of fn at three different axial strains and confining pressures. The statistics are carried out over all broken bonds. Note that fn = 0 represents pure tension failure and fn = 1

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

σ1−σ3 (kPa)

1200 1000 800 600 400 200 0 Uncemented 0.5% 1% 2% 3%

20 15 10

1200 1000 800 600 400 200 0

-1.0 -0.5

5 0

Uncemented 0.5% 1% 2% 3%

-1.5

− dεv /dε1

Bond breakage rate

σ1−σ3 (kPa)

206

0.0 0

1

2

3

4

0.5

5

0

1

2

3

4

5

ε1 (%)

ε1 (%)

(a) Stress and bond breakage rate

(b) Stress and dilation rate

2.5

2.5

2.0

2.0

1.5

1.5

η = q/p

η = q/p

Fig. 20. Evolutions of stress, bond breakage rate and dilation rate with axial strain (e0 = 0.82, r3 = 100 kPa, bond content = 0.5%, 1%, 2%, 3%). Arrows locate the peak deviator stresses and circles locate the maximum bond breakage rates or the maximum dilation rates.

1.0 50kPa 100kPa 200kPa 400kPa

0.5 0.0 -0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0 50kPa 100kPa 200kPa 400kPa

0.5 0.0 -0.6

1.0

-0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

D = − d ε v /d ε s

D = − d ε v /d ε s

(b) Bond content=2%, e0 = 0.91

(a) Bond content=1%, e0 = 0.82

Fig. 21. Stress–dilatancy relationship in triaxial compression tests on samples with different bond contents and void ratios.

800

σ1−σ3 (kPa)

σ1−σ3 (kPa)

800 600 400 200

600 400 200

0

0

-5

-8 -6

ε v (%)

ε v (%)

-4 -3 -2 -1

-2 0

0 1

-4

0

5

10

15

ε1 (%)

(a) e0 = 0.91

20

25

2

0

5

10

15

20

25

30

35

ε1 (%)

(b) e0 = 0.82

Fig. 22. Behavior with unbreakable bonded contacts in samples with different void ratios (r3 = 100 kPa, bond content = 1%).

207

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

15 10 5 0.05

0.10

0.15

0.20

Borken bond (%)

20

0 0.00

25

25

Borken bond (%)

Borken bond (%)

25

20 15 10 5 0 0.00

0.25

0.05

0.10

0.15

0.20

20 15 10 5 0 0.00

0.25

0.05

0.10

0.15

0.20

fn

fn

fn

Axial strain = 3.3%

Axial strain = 9.7%

Axial strain = 32%

0.25

25

25

20

20

20

15 10 5 0 0.00

0.05

0.10

0.15

0.20

0.25

Borken bond (%)

25

Borken bond (%)

Borken bond (%)

(a) σ 3 = 100kPa

15 10 5 0 0.00

0.05

0.10

fn

0.15

0.20

15 10 5 0 0.00

0.25

0.05

0.10

fn σ 3 = 200kPa

σ 3 = 50kPa

0.15

0.20

0.25

fn σ 3 = 400kPa

(b) Axial strain = 32% Fig. 23. Frequency distribution of dimensionless normal force at the instant of bond breakage in triaxial compression tests (bond content = 1%, e0 = 0.82).

Fig. 24. Triangular representation of bond failure condition for coupled failure.

pure compression failure. However, pure tension and compression failures are two extreme cases which rarely occur. Fig. 23 shows that, regardless of the axial strain and confining pressure, fn is smaller than 0.25 and mostly concentrates in the range (0, 0.15), implying that the normal force at the instant of bond breakage is far from causing compression failure. Therefore, the macroscopic shear strength is less sensitive to compressive strength rc than to tensile strength rt, as observed in Fig. 13. 2

2

2

In the coupled failure criterion given by Eq. (13), f s , f r , and f t represent the respective contribution of shear, rolling and torsion modes to the overall failure of a bond. Note that, for an equilateral triangle, the summation of distances from any interior point to the three sides must equal to the altitude of this triangle. By constructing an equilateral triangle with an altitude of unit length, 2

2

2

the failure condition ðf s þ f r þ f t ¼ 1Þ can be represented by a point inside the triangle as shown in Fig. 24. For any interior

2

point P in Fig. 24, its distance to the side \f s " represents the 2

value of f s and the same applies to the other two sides. Each one-sixth sector of this equilateral triangle represents one possible 2

2

2

sequence of contribution to bond failure from f s , f r , and f t . A vertex of the triangle represents the shear, rolling or torsion failure mode. Fig. 25 presents the mapping of bond failure conditions at three different axial strains and confining pressures. Each dot represents 2

a broken bond. These dots concentrate along the \f t " side and 2 \f t "

around the two corners on the side, implying that the torsion mode contributes little whereas the shear and rolling modes are the two major contributions to bond failure. The dot density is almost symmetric about the perpendicular bisector (dash line in 2

Fig. 25) of the \f t " side, indicating that the effect of shear and rolling modes are statistically equal.

208

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209

Axial strain = 3.3%

Axial strain = 9.7%

Axial strain = 32%

(a) σ 3 = 100kPa

σ 3 = 50kPa

σ 3 = 200kPa

σ 3 = 400kPa

(b) Axial strain = 32% Fig. 25. Bond failure condition for coupled failure obtained from triaxial compression tests (bond content = 1%, e0 = 0.82).

7. Conclusions A 3D bonded contact model has been developed by modifying a previous unbonded contact model to incorporate the effect of bond material on contact behavior. The new model has been incorporated into a DEM code to simulate the mechanical behavior of cemented sand. Based on observations of the micro-structure of cemented sand, two types of bonded contact, parallel bond contact and serial bond contact, are identified and modeled. The bond geometry effect is captured in terms of bond strength and stiffness. The complete interactions in the normal, tangential, rolling and torsional directions are considered, and the coupled effect of loads in the four directions is taken into account in the bond failure criterion. By considering the effect of attached post-failure bond material on the particle surface, the bonded contact model can reduce gradually to an unbonded contact model. These new features have never been considered fully in previous models. The readers would have the flexibility to choose some of the features in this model if not all features are necessary in their specific research. Also, the model in this paper could serve as a reference to compare different contact models in future work. Several conclusions can be reached. (1) A parametric study of DEM simulations shows that the bonded contact model can capture the major behavior of cemented sands. The presence of a small amount of attached bond material can greatly increase the modulus, strength and dilation of sand. (2) With increase in bond content or decrease in confining pressure, the bond effect becomes significant and the mechanical behavior evolves from dilation-dominant, where the peak state coincides with the maximum dilation rate, to bonddominant where the peak state coincides with the maximum bond breakage rate. Bond breakage plays a crucial role in both mechanisms.

(3) Bonds fail in a coupled way due to shear, rolling and torsional loads, but the shear and rolling modes are the two statistically equal dominant modes. The normal force at the instant of bond breakage is far from causing compression failure.

Acknowledgements The research is funded by China National Funds for Distinguished Young Scientists with Grant No. 51025932, which is sincerely appreciated. The authors would also thank Prof. Richard Wan for his help during the first author’s visit to University of Calgary. References [1] Collins BD, Sitar N. Stability of steep slopes in cemented sands. J Geotech Geoenviron Eng 2011;137(1):43–51. [2] Coop MR, Atkinson JH. The mechanics of cemented carbonate sands. Géotechnique 1993;43(1):53–67. [3] Leroueil S, Vaughan PR. The general and congruent effects of structure in natural soils and weak rocks. Géotechnique 1990;40:467–88. [4] Airey DW. Triaxial testing of naturally cemented carbonate soil. J Geotech Eng (ASCE) 1993;119(9):1379–98. [5] Consoli NC, Vendrusco MRA, Prietto PDM. Behavior of plate load tests on soil layers improved with cement and fiber. J Geotech Geoenviron Eng 2003;129 (1):96–101. [6] Rattley M, Lehane B, Consoli N, Richards D. Uplift of shallow foundations with cement-stabilised backfill. Proc ICE-Ground Improve 2008;161(2):103–10. [7] Cuccovillo T, Coop M. On the mechanics of structured sands. Géotechnique 1999;49(6):741–60. [8] Huang JT, Airey DW. Properties of artificially cemented carbonate sand. J Geotech Geoenviron Eng 1998;124(6):492–9. [9] Haeri SM, Hamidi A. Constitutive modelling of cemented gravelly sands. Geomech Geoeng 2009;4(2):123–39. [10] Pekau OA, Gocevski V. Elasto-plastic model for cemented and pure sand deposits. Comput Geotech 1989;7(3):155–87. [11] Abdulla AA, Kiousis PD. Behavior of cemented sands-II. Modelling. Int J Numer Anal Met Geomech 1997;21(8):549–68.

Z. Shen et al. / Computers and Geotechnics 75 (2016) 192–209 [12] Lagioia R, Nova R. An experimental and theoretical study of the behavior of a calcarenite in triaxial compression. Géotechnique 1995;45(4):633–48. [13] Nova R, Castellanza R, Tamagnini C. A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation. Int J Numer Anal Met Geomech 2003;27(9):705–32. [14] Desai CS, Toth J. Disturbed state constitutive modeling based on stress-strain and nondestructive behavior. Int J Solids Struct 1996;33(11):1619–50. [15] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979;29(1):47–65. [16] Brown NJ, Chen J-F, Ooi JY. A bond model for DEM simulation of cementitious materials and deformable structures. Granul Matter 2014;16(3):299–311. [17] de Bono J, McDowell G, Wanatowski D. Investigating the micro mechanics of cemented sand using DEM. Int J Numer Anal Met Geomech 2014. http://dx.doi. org/10.1002/nag.2340. [18] Estrada N, Taboada A. Yield surfaces and plastic potentials of cemented granular materials from discrete element simulations. Comput Geotech 2013;49:62–9. [19] Jiang MJ, Zhang FG, Sun YG. An evaluation on the degradation evolutions in three constitutive models for bonded geomaterials by DEM analyses. Comput Geotech 2014;57:1–16. [20] Jiang MJ, Zhang WC, Sun YG, Utili S. An investigation on loose cemented granular materials via DEM analyses. Granul Matter 2012;15(1):65–84. [21] Utili S, Nova R. DEM analysis of bonded granular geomaterials. Int J Numer Anal Met Geomech 2008;32(17):1997–2031. [22] Obermayr M, Dressler K, Vrettos C, Eberhard P. A bonded-particle model for cemented sand. Comput Geotech 2013;49:299–313. [23] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min 2004;41(8):1329–64. [24] Wang YH, Leung SC. A particulate-scale investigation of cemented sand behavior. Can Geotech J 2008;45(1):29–44. [25] Wang YH, Leung SC. Characterization of cemented sand by experimental and numerical investigations. J Geotech Geoenviron Eng 2008;134(7):992–1004. [26] Brendel L, Torok J, Kirsch R, Brockel U. A contact model for the yielding of caked granular materials. Granul Matter 2011;13(6):777–86. [27] Carmona HA, Wittel FK, Kun F, Herrmann HJ. Fragmentation processes in impact of spheres. Phys Rev E 2008;77(5):051302. [28] Kun F, Herrmann HJ. A study of fragmentation processes using a discrete element method. Comput Method Appl Mech Eng 1996;138(1–4):3–18. [29] Wang YC, Alonso-Marroquin F. A finite deformation method for discrete modeling: particle rotation and parameter calibration. Granul Matter 2009;11 (5):331–43. [30] Delenne JY, El Youssoufi MS, Cherblanc F, Benet JC. Mechanical behaviour and failure of cohesive granular materials. Int J Numer Anal Met Geomech 2004;28 (15):1577–94. [31] Kirsch R, Bröckel U, Brendel L, Török J. Measuring tensile, shear and torsional strength of solid bridges between particles in the millimeter regime. Granul Matter 2011;13(5):517–23. [32] Jiang MJ, Sun YG, Xiao Y. An experimental investigation on the mechanical behavior between cemented granules. Geotech Test J 2012;35(5):678–90. [33] Jiang MJ, Sun YG, Li LQ, Zhu HH. Contact behavior of idealized granules bonded in two different interparticle distances: an experimental investigation. Mech Mater 2012;55:1–15. [34] Jiang MJ, Liu F, Zhou YP. A bond failure criterion for DEM simulations of cemented geomaterials considering variable bond thickness. Int J Numer Anal Met Geomech 2014;38(18):1871–97. [35] Jiang MJ, Sun YG, Yang QJ. A simple distinct element modeling of the mechanical behavior of methane hydrate-bearing sediments in deep seabed. Granul Matter 2013;15(2):209–20. [36] Shen ZF, Jiang MJ, Wan R. Numerical study of inter-particle bond failure by 3D discrete element mothod. Int J Numer Anal Met Geomech 2016;40(4):523–45. [37] Jiang MJ, Shen ZF, Wang JF. A novel three-dimensional contact model for granulates incorporating rolling and twisting resistances. Comput Geotech 2015;65:147–63. [38] Ismail MA, Joer HA, Randolph MF, Meritt A. Cementation of porous materials using calcite. Géotechnique 2002;52(5):313–24. [39] Cuccovillo T, Coop M. Yielding and pre-failure deformation of structured sands. Géotechnique 1997;47(3):491–508. [40] Dvorkin J, Mavko G, Nur A. The effect of cementation on the elastic properties of granular material. Mech Mater 1991;12(3–4):207–17. [41] Dvorkin J, Nur A, Yin H. Effective properties of cemented granular materials. Mech Mater 1994;18(4):351–66. [42] Dvorkin J, Yin HZ. Contact laws for cemented grains: implications for grain and cement failure. Int J Solids Struct 1995;32(17–18):2497–510. [43] Yin HZ, Dvorkin J. Strength of cemented grains. Geophys Res Lett 1994;21 (10):903–6. [44] Zhu H, Chang CS, Rish JW. Normal and tangential compliance for conforming binder contact I: elastic binder. Int J Solids Struct 1996;33(29):4337–49. [45] Zhu H, Chang CS, Rish JW. Normal and tangential compliance for conforming binder contact II: visco-elastic binder. Int J Solids Struct 1996;33(29):4351–63.

209

[46] Zhu H, Chang CS, Rish JW. Rolling compliance for elastic and visco-elastic conforming binder contact. Int J Solids Struct 1997;34(31–32):4073–86. [47] Pilkavicˇius S, Kacˇianauskas R, Norkus A. Investigation of normal contact interaction between two bonded spherical particles with interface layer. Mechanika 2012;18(6):632–9. [48] Wong TF, Wu LC. Tensile-stress concentration and compressive failure in cemented granular material. Geophys Res Lett 1995;22(13):1649–52. [49] Zang A, Wong TF. Elastic stiffness and stress-concentration in cemented granular material. Int J Rock Mech Min 1995;32(6):563–74. [50] Zhong XX, Chang CS. Micromechanical modeling for behavior of cementitious granular materials. J Eng Mech (ASCE) 1999;125(11):1280–5. [51] Miyazaki K, Masui A, Sakamoto Y, Aoki K, Tenma N, Yamaguchi T. Triaxial compressive properties of artificial methane-hydrate-bearing sediment. J Geophys Res 2011;116(B6). [52] Masayuki H, Jun Y, Norimasa Y, Yukio N. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed. Soils Found 2013;53(2):299–314. [53] Itasca Consulting Group. PFC3D (particle flow code in three dimensions) user’s guide. Minneapolis; 2008. [54] Haeri SM, Hosseini SM, Toll DG, Yasrebi SS. The behaviour of an artificially cemented sandy gravel. Geotech Geol Eng 2005;23(5):537–60. [55] Lo SR, Wardani SP. Strength and dilatancy of a silt stabilized by a cement and fly ash mixture. Can Geotech J 2002;39(1):77–89. [56] Schnaid F, Prietto PDM, Consoli NC. Characterization of cemented sand in triaxial compression. J Geotech Geoenviron Eng 2001;127(10):857–68. [57] Lambe TW. A mechanistic picture of shear strength in clay. In: Research conference on shear strength of cohesive soils (ASCE); 1960. p. 555–80. [58] Abdulla AA, Kiousis PD. Behavior of cemented sands-I. Testing. Int J Numer Anal Met Geomech 1997;21(8):533–47. [59] Clough GW, Sitar N, Bachus RC, Rad NS. Cemented sands under static loading. J Geotech Eng Div (ASCE) 1981;107(6):799–817. [60] Marri A, Wanatowski D, Yu HS. Drained behaviour of cemented sand in high pressure triaxial compression tests. Geomech Geoeng 2012;7(3):159–74. [61] Asghari E, Toll D, Haeri S. Triaxial behaviour of a cemented gravely sand, tehran alluvium. Geotech Geol Eng 2003;21(1):1–28. [62] Consoli N, Montardo J, Prietto P, Pasa G. Engineering behavior of a sand reinforced with plastic waste. J Geotech Geoenviron Eng 2002;128(6):462–72. [63] Consoli NC, Vendruscolo MA, Fonini A, Dalla RosaF. Fiber reinforcement effects on sand considering a wide cementation range. Geotext Geomembr 2009;27 (3):196–203. [64] Kaniraj SR, Havanagi VG. Behaviour of cement-stabilized fiber-reinforced fly ash-soil mixtures. J Geotech Geoenviron Eng 2001;127(7):574–84. [65] Tang CS, Shi B, Gao W, Chen FJ, Cai Y. Strength and mechanical behavior of short polypropylene fiber reinforced and cement stabilized clayey soil. Geotext Geomembr 2007;25(3):194–202. [66] Ismail MA, Joer HA, Sim WH, Randolph MF. Effect of cement type on shear behavior of cemented calcareous soil. J Geotech Geoenviron Eng 2002;128 (6):520–9. [67] Dano C, Hicher PY, Tailliez S. Engineering properties of grouted sands. J Geotech Geoenviron Eng 2004;130(3):328–38. [68] Porcino D, Marcianò V, Granata R. Static and dynamic properties of a lightly cemented silicate-grouted sand. Can Geotech J 2012;49(10):1117–33. [69] Rahman ZA, Toll DG, Gallipoli D, Taha MR. Micro-structure and engineering behaviour of weakly bonded soil. Sains Malays 2010;39(6):989–97. [70] O’Rourke T, Crespo E. Geotechnical properties of cemented volcanic soil. J Geotech Eng 1988;114(10):1126–47. [71] Consoli NC, da Fonseca AV, Cruz RC, Heineck KS. Fundamental parameters for the stiffness and strength control of artificially cemented sand. J Geotech Geoenviron Eng 2009;135(9):1347–53. [72] Coop M, Willson S. Behavior of hydrocarbon reservoir sands and sandstones. J Geotech Geoenviron Eng 2003;129(11):1010–9. [73] Malandraki V, Toll D. Drained probing triaxial tests on a weakly bonded artificial soil. Géotechnique 2000;50(2):141–51. [74] Malandraki V, Toll DG. Triaxial tests on weakly bonded soil with changes in stress path. J Geotech Geoenviron Eng 2001;127(3):282–91. [75] Rios S, Viana da Fonseca A, Baudet BA. On the shearing behaviour of an artificially cemented soil. Acta Geotech 2014;9(2):215–26. [76] Shen ZF, Jiang MJ. DEM simulation of bonded granular material. Part II: Extension to grain-coating type methane hydrate bearing sand. Comput Geotech. 2016 [accepted]. [77] Jiang MJ, Jin SL, Shen ZF, Coop MR, Liu W. Experimental study on threedimensional contact behavior of bonded granules. In: 2015 International Symposium on Geohazards and Geomechanics, IOP Conference Series: Earth and Environmental Science, vol. 26; 2015. doi: http://dx.doi.org/10.1088/ 1755-1315/26/1/012007. [78] Chang SH, Lee CI. Estimation of cracking and damage mechanisms in rock under triaxial compression by moment tensor analysis of acoustic emission. Int J Rock Mech Min 2014;41(7):1069–86.