The 3D numerical simulation of damage localization of rocks using General Particle Dynamics

Engineering Geology 224 (2017) 29–42

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

The 3D numerical simulation of damage localization of rocks using General Particle Dynamics

MARK

X.P. Zhoua,b,c,⁎, X.Q. Zhanga,b,c a b c

State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400045, PR China School of Civil Engineering, Chongqing University, Chongqing 400045, PR China Key Laboratory of New Technology for Construction of Cities in Mountain Area, Chongqing University, Chongqing 400045, PR China

A R T I C L E I N F O

A B S T R A C T

Keywords: General Particle Dynamics (GPD) The failure of rocks Damage localization 3D numerical simulation The intermediate principal stress

The damage localization of rocks is an open issue in rock engineering. Moreover, it is difficult to observe directly the micro-scale process of damage localization and damage evolution of rocks. Therefore, in this paper, the General Particle Dynamic method (GPD) is developed to simulate the damage localization and failure process of 3D rock samples. The effect of the intermediate principal stress on the strength, the stress-strain curve and the damage localization modes of a cubical rock specimen subjected to the true triaxial compressive loads is investigated. The different failure modes including the single shear damage localization mode, axial splitting damage mode, X-shaped conjugate shear damage localization mode and Cone-shaped damage localization are found. It is found that the numerical results obtained from GPD are in good agreement with experimental observations. It is implied that the formation of damage localization band induced by the initiation and propagation and coalescence of micro-cracks in rocks can be well simulated using GPD.

1. Introduction It is well known that rock materials undergo damage localization and failure under the compressive loads. The better understanding of the damage localization of rocks is of a significant importance due to its widespread applications to many types of geotechnical and geological engineering. To study the local damage and failure of rocks, many experimental investigations were performed. For instance, Bésuelle et al. (2000) conducted triaxial compression experiments on sandstone and analyzed the development of damage localization phenomenon. Klein et al. (2001) carried out an experimental study on sandstone under conventional triaxial compression, and studied the effect of confining pressure on the damage localization process. Yang and Jing (2013) conducted triaxial compression experiments on red sandstone to investigate its strength and deformation failure behaviors. In parallel with the experimental studies, various theoretical models were also proposed to describe the damage localization and failure in rock materials. A micromechanics-based model was proposed to analyze the relationship between the localization damage and deformation of brittle rocks subjected to unloading (Zhou, 2005). Hu et al. (2010) studied an anisotropic plastic damage model for semi-brittle materials, and revealed that the evolution of damage is related to growth of weakness

⁎

Corresponding author. E-mail address: [email protected] (X.P. Zhou).

http://dx.doi.org/10.1016/j.enggeo.2017.04.021 Received 7 July 2016; Received in revised form 2 March 2017; Accepted 25 April 2017 Available online 11 May 2017 0013-7952/ © 2017 Elsevier B.V. All rights reserved.

planes. Zhu and Shao (2015) remarked that inelastic deformation and damage evolution at microdefects can govern macroscopic behaviors of brittle solids, and raised a refined micromechanical anisotropic unilateral damage model for quasi-brittle geomaterials under compression. Along this line, a micromechanical model considering the damagefriction coupling for brittle materials was firstly proposed by Zhu et al. (2016) and Qi et al. (2016). It is found that the coupled model can be more accurate description of the damage localization evolution law. These physical and empirical models provide quite a direct interpretation of damage localization and evolution. However, the formulation of such models is complex, and the experimental identification of parameters is not easy. Therefore, in order to describe more conveniently and more properly the damage localization of materials, a number of numerical tools were proposed to study the process from the damage localization to failure of rock-like materials. The localized degradation model was adopted to study the failure of the rock using FLAC3D (Fang and Harrison, 2002). In order to investigate the failure mechanism and mode of rocks, Wong et al. (2006) carried out a model which contains pre-existing flaws under uniaxial compressive loads using Rock Failure Process Analysis code (RFPA2D). Pan et al. (2012) studied the effect of the intermediate principal stress on rock failure process using ElastoPlastic Cellular Automaton (EPCA3D), and found that local failure is

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 3. Geometry model. Fig. 1. Failure loci in the deviatoric plane for the nonlinear unified strength criterion. (a) Undamaged configuration. (b) Damaged configuration.

formed when the intermediate principal stress reaches a certain value. Wu et al. (2015) studied the damage localization for solid using the extended embedded finite elements (XFEM). Shao et al. (2006) and Jiang et al. (2011) investigated the deformation localization and crack formation in elastic-brittle materials using the discrete numerical method (DEM). Although the above numerical methods succeeded in one or another aspects, most of these methods were based on the finite element (FEM) or discrete element methods (DEM), which have many limitations. For instance, the numerical modeling of fracture with traditional mesh-based technique (e.g. FEM) requires a very fine mesh to model the singular stress field in the neighborhood of damaged material (e.g. crack tip) (Zhou et al., 2015a,b). As the damage evolves, the structure needs to be re-meshed to consider the localized change in geometry and certain inaccurate results (Fernández-Méndez et al., 2005). In addition, it is very difficult to investigate the 3D localization phenomenon of geomaterials using FEM. Moreover, it is found that the particle size has a strong influence on numerical results in the discrete element-based numerical methods (DEM) (Yao et al., 2016). Therefore, there is a need to develop other efficient numerical methods for minimizing effect of mesh size and particle size. It is attractive to adopt the mesh-less approach to simulate damage localization. The advantage of this approach is that no remeshing is needed. A number of mesh-less methods including Smoothed Particle Hydrodynamics (SPH) (Leblanc et al., 2014), the continuous/discontinuous deformation analysis (CDDA) method (Cai et al., 2014), the meshless Shepard and least squares (MSLS) (Zhu et al., 2011), a numerical

Fig. 4. The particle model.

procedure based on the edge-based smoothed finite element (Vu-Bac et al., 2013) and the three dimensional discontinuous deformation analysis (3D DDA) (Wu et al., 2014) were developed. In the above mesh-less methods, SPH has the advantage of robustly computing the material point history even at severely deformed configurations while avoiding the large computational costs of re-meshing in a Lagrangian

(a) Undamaged configuration.

(b) Damaged configuration.

Fig. 2. The distribution of the 3D particles.

30

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Experiment data

1200

The peak point 1000

σ1-σ3/ MPa

800 600

Bifurcation point 400 200 0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

ε1/ %

(a) Experiment results (Lee and Haimson, 2011). 1200

The peak point

B

800

Rock

Density (Kg/ m3)

Uniaxial compressive strength (MPa)

Young's modulus (GPa)

Poisson's ratio

2707.8

161.8

49

0.24

σ1-σ3/ MPa

Table 1 The main material parameters of granodiorite. Intermediate stress parameter

E

C

1000

Fig. 5. The uniaxial strength distribution of 3D particles in the numerical model. (a) Experiment results (Lee and Haimson, 2011). (b) The numerical results.

GPD damage model

D

F 600

Bifurcation point 400

A

200 Granodiorite

1.0

0 0.0

framework or an Eulerian framework. However, for SPH method, the interaction between any two particles is only controlled by the kernel function, the interaction is automatically terminated if one is out of the influence domain of the other. Therefore, SPH can not well simulate the damage localization and failure process of materials. In this paper, General Particle Dynamics (GPD), which overcomes the shortcomings of SPH, is developed to simulate the damage localization of rocks. In General Particle Dynamics method, a new parameter, named as the “interaction factor”, is introduced to define the level of interaction between the ith and jth particles. It is assumed that any particle is damaged when its stresses satisfy the critical value, and the new boundaries are formed between damaged particles. Therefore, GPD has the advantage of robustly simulating the damage localization process of rocks subjected to true triaxial compressive loads. This paper is organized as follows: The damage model of General Particle Dynamics (GPD) is given in Section 2. The numerical model of three-dimensional particle distribution is described in Section 3. The damage localization process of rocks subjected to true triaxial compressive load is studied in Section 4. The effect of the intermediate principal stress on the damage localization process of rock is investigated in Section 5. The conclusions are then drawn in Section 6.

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

ε1/ %

(b) The numerical results. Fig. 6. The complete stress-strain curves obtained from: (a) Experiment data (Lee and Haimson, 2011); (b)The numerical results.

the damage initiation and growth in rocks are determined by the damage particles. In this paper, the non-linear unified strength (failure) criterion is applied to determine the condition of damage initiation. Damage is believed to be initiation from one particle when its stresses satisfy the non-linear unified strength (failure) criterion (Yu et al., 2002). It is clear that there are three principal shear stresses τ13, τ12 andτ23 in the three-dimensional principal stress state σ1, σ2 and σ3. However, only two principal shear stresses are independent variables among τ13, τ12, τ23 because the maximum principal shear stress equals the sum of the other two, that is to say, τ13 = τ12 + τ23. Since there are only two independent principal shear stresses, the shear stress state can also be converted into the twin-shear stress state (τ13, τ12, σ13, σ12) or (τ13, τ23, σ13, σ23). Therefore, the twin-shear strength theory was established (Yu, 1983). However, the twin-shear strength theory (Yu, 1983) can not reflect the fundamental characteristics of rocks, i.e., different tensile and compressive strengths, hydrostatic pressure effects and nonlinear behaviors. Thus, the non-linear unified strength (failure) criterion is proposed as follows (Yu et al., 2002):

2. The damage model of General Particle Dynamics (GPD) 2.1. Damage strength criterion

⎡ ⎤2 1 mσc F = ⎢σ1 − (bσ2 + σ3 ) ⎥ − × (bσ2 + σ3) − sσc2 = 0 when F ⎣ ⎦ 1+b 1+b

In this paper, only the damage model implemented in the GPD is described in detail because the general formulations of this method such as constitutive model, governing equations, correction for consistency were stated in the previous works (Zhou et al., 2015a,b). In most cases, rock-like materials fail in a brittle manner. Therefore,

(1)

≥ F′

⎡ 1 ⎤2 F′ = ⎢ (bσ2 + σ1) − σ3 ⎥ − mσc σ3 − sσc2 = 0 ⎣1 + b ⎦ 31

when F < F′

(2)

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 7. The distribution of the maximum principal stress σ1 at different time steps (Unit: Pa).

Fig. 8. Comparison between numerical results and experimental data: (a) A macro-failure mode represented by arrows of granodiorite specimen subjected to true triaxial compressive loads. (b) The fault plane profiles and damage localization during true triaxial tests from the SEM micrographs. (c)The single shear damage localization mode obtained from the proposed numerical method.

32

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

σ1/σ2/σ3=190/50/10MPa

C

C

200

B

σ1-σ3/ MPa

particles. This interaction factor f is determined based on the damage in particles. Initially, for undamaged particles, f = 1, which implies ‘full interaction’. With the progress of damage, f finally becomes zero for fully damaged particle.

σ1/σ2/σ3=220/80/30MPa

250

σ1/σ2/σ3=175/15/10MPa

C

B D

150

B

D 100

D

A E A

0 0.0

E 0.2

0.4

0.6

0.8

1.0

ε1/ % Fig. 9. The numerical stress-strain curve of numerical sample under the different stress condition. (a) X-shaped conjugate damage mode. (b) Axial splitting damage mode. (c) Cone-shaped damage mode on left of the rock specimen.

(if σi < σmax )

(5)

D = 1, and

f = 0,

(if particle damaged)

(6)

dρi = dt

where b is the intermediate stress parameter, which can be determined σ by (1 + α ) τ0 − σt , in which α = σt , σt , σc are respectively the uniaxial σt − τ 0

f = 1,

where D is the damage factor, and f is the interaction factor. In the proposed method, influence of damaged particles on neighbors is proposed to determine the level of interaction between the ith and the jth particles. The statistical distribution of the mechanical parameters of particles is described by the Weibull distribution function, and effects of the heterogeneity of rock materials are considered. Therefore, GPD can simulate the damage localization process of brittle rocks subjected to true triaxial compressive loads. Then, the discrete conservation equations for General Particle Dynamics (GPD) are expressed as (Zhou et al., 2015a,b):

E 50 A

D = 0, and

∑ mj vijβ Wij,β + f ⋅ ∑ j∈U

c

tensile strength and uniaxial compressive strength of rocks; σc is uniaxial compressive strength of an intact rock material;τ0is shear strength of rocks; m and s are the material parameter as the same as those in the Hoek-brown criterion, which are defined by Hoek (1990) and taken the following form:

dviα =− dt

⎛ GIS − 100 ⎞ m = mi exp ⎜ ⎟ ⎝ 28 − 14D1 ⎠

1 dei = 2 dt

⎛ GIS − 100 ⎞ s = exp ⎜ ⎟ ⎝ 9 − 3D1 ⎠

mj vijβ Wij, β (7)

j ∈ Da

⎛ σ αβ

∑ mj ⎜⎜ j∈U

⎝

i

ρi2

+

σ jαβ ρj2

⎞ + Πij⎟⎟⋅Wij, β − f ⋅ ⎠

⎞ ⎛ σ αβ σ jαβ mj ⎜⎜ i 2 + 2 + Πij⎟⎟ ρj ⎠ j ∈ Da ⎝ ρi

∑

⋅Wij, β

(3)

⋅ (4)

where D1 is a disturbance coefficient which varies from 0.0 for the undisturbed in situ rock masses to 1.0 for very disturbed rock masses (Hoek, 1990), mi is the value of m for intact rock and can be obtained from experiments, the parameter mi varies from 4 for very fine weak rock like clay-stone to 33 for coarse igneous light-colored rock like granite, and in this paper m = 10 and s = 0.5. In Eq. (1), the physical meaning of F is that the strength of rocks is controlled by the shear stress state of (τ13, τ12, σ13, σ12). That is to say, τ12 ≥ τ23, the failure of rocks is controlled by the maximum principal shear stress τ13 and the intermediate principal shear stress τ12. In Eq. (2), the physical meaning of F′ is that the strength of rocks is controlled by the shear stress state of (τ13, τ23, σ13, σ23).That is to say, τ23 > τ12, the failure of rocks is controlled by the maximum principal shear stress τ13 and the intermediate principal shear stress τ23. It is evident that the intermediate principal stress is considered in the non-linear unified strength criterion. The non-linear twin-shear failure criterion can be introduced by the non-linear unified strength criterion when the parameter b = 1.0 and the single-shear strength criterion can be reduced when b = 0. According to Eqs. (1) and (2), the limit loci in the deviatoric plane can be plotted in Fig. 1. The limit locus of the single-shear strength criterion is the lower bound in the deviatoric plane (b = 0). The upper bound is the non-linear twin-shear strength criterion (b = 1.0). Therefore, the non-linear unified strength criterion forms a full spectrum of new criteria when b takes values between 0 and 1.0 (i.e. 0 < b < 1), which in turn can suit for various different materials.

1 2

(8)

⎛ σ αβ

∑ mj vijα ⎜⎜ j∈U

⎝

i

ρi2

+

σ jαβ ρj2

⎞ + Πij ⎟⎟ Wij, β + f ⎠

⎛ σ αβ ⎞ σ jαβ mj vijα ⎜⎜ i 2 + 2 + Πij ⎟⎟ Wij, β ρj j ∈ Da ⎝ ρi ⎠

∑

(9)

where U denotes undamaged particles, Da denotes damaged particles, mj is the mass of the jth particle, ρi, viα and ei are respectively the mass density, velocity vector and specific internal energy of the ith particle, σαβ denotes the Cauchy stress tensor (σ) and it is positive if the stress is ∂W (xj − xi , h ) is the gradient of kernel function with smoothtensile, Wij, β = β ∂xi

ing length h, vijα = viα − vjα, Πij is the artificial viscosity which can be used to stabilize computation in shock problems, d/dt is the derivative versus time used in the moving Lagrangian framework, the superscripts α , β = 1 , 2 , 3 are integer indices for the three spatial directions. Once damage is initiated in one particle, its interactions will no longer be the same as those in the undamaged material. Although the damaged particles do not disappear and are still in the influence domain of others, the damaged particles are traction free and the damaged particles have no influence on the surrounding particles in the stress calculation because the interaction factor f is zero. For the damaged particles, because their stress becomes zero, the surrounding living particles have no influence on them. A new boundary is formed, as shown in Fig. 2.

3. The numerical model of three-dimensional particle distribution 3.1. The numerical model and loading conditions In this paper, damage localization features of rocks subjected to true triaxial compressive loads are modeled. The numerical model of 50 × 100 × 30 = 150,000 particles represents a sample geometry of 50 (width) × 100 (length) × 30 (thickness) mm in scale, and the loading conditions for a rock specimen subjected to the true triaxial compressive loads are plotted in Fig. 3. The rock specimen is discretized into a system composed of 3D particles, as shown in Fig. 4.

2.2. Damage model Now we introduce a parameter f, coined as the ‘interaction factor’ which defines the level of interaction between the ith and the jth 33

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

(a) X-shaped conjugate damage mode.

(b) Axial splitting damage mode.

(c) Cone-shaped damage mode on left of the rock specimen. Fig. 10. The maximum principal stress contour of samples at different time steps (Unit: Pa). (a) σ1/σ2/σ3 = 175/15/10 MPa; (b) σ1/σ2/σ3 = 190/50/10 MPa; (c) σ1/σ2/σ3 = 220/80/ 30 MPa.

3.2. Weibull distribution function

simulate the random microstructures in rock materials, rock heterogeneity can be well characterized by using statistical approaches. In GPD method, it is assumed that particles are of the same shape and size, and that there is no geometric priority in any orientation (Zhou et al.,

When a heterogeneous rock material is involved, the disorder of microstructures in rocks should be modeled in a numerical model. To 34

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 11. Sketch of the numerical macro damage localization mode under the different compressive stress. (a) X-shaped conjugate damage mode. (b) Axial splitting damage mode. (c) Cone-shaped damage mode.

(a) X-shaped conjugate damage mode

(b) Axial splitting damage mode

(c) Cone-shaped damage mode

Fig. 12. The experimental failure mode of specimens under the different compressive stress (Chen and Feng, 2006).

GPD resolution 1.0mm GPD resolution 1.25mm GPD resolution 1.5mm

200

W (x ) =

σ1-σ3/ MPa

(10)

where, m1 defines the shape of the Weibull distribution function, and it can be referred to as the homogeneity index, x is the mechanical parameter of one particle, and x0 (the uniaxial compressive strength) is the even value of the parameter of all the particles. According to the Weibull distribution (Weibull, 1951), a larger value of m1 indicates that more particles have mechanical properties that are approximated to the mean value, which describes a more homogeneous rock sample. In the present numerical models, the heterogeneity index is m1 = 10, as shown in Fig. 5.

150

100

50

0 0.0

m −1 ⎡ ⎛ x ⎞m1⎤ m1 ⎛ x ⎞ 1 exp ⎢ −⎜ ⎟ ⎥ ⎜ ⎟ x0 ⎝ x0 ⎠ ⎣ ⎝ x0 ⎠ ⎦

0.2

0.4

0.6

0.8

4. Damage localization process of rocks subjected to true triaxial compressive loads

1.0

ε1/ %

4.1. The numerical results of the damage localization process of 3D specimens

Fig. 13. Steady stress-strain variation for different particle resolutions.

In this paper, the main objective is to analyze the damage localization process of granodiorite subjected to the true triaxial compression loads using the GPD. In order to verify the proposed numerical method, the numerical condition is consistent with the test one. According to works by Lee and Haimson (2011), the mechanical properties of granodiorite are listed in Table 1. The true triaxial test data show that b = 0.5–1.0 for the dolomite and the marble, b = 1.0 for granodiorite, granite and the trachyte (Yu et al., 2002). Therefore, the value of b is taken 1.0. True triaxial compression tests were conducted by simultaneously raising all three principal stresses at a constant rate until σ3

2015a,b). Disorder can be defined by randomly distributing the mechanical properties of particles. The statistical distribution of the elemental mechanical parameters is described by the Weibull distribution function (Weibull, 1951). These mechanical parameters of particles include the uniaxial compressive strength, Young's modulus, Poisson ratio and unit weight. For the localized damage problem, only the uniaxial compressive strength of particles is described by the Weibull's distribution. The Weibull distribution function is expressed as follow (Weibull, 1951): 35

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 14. The maximum principal stress contour of samples at the different particle resolutions (Unit: Pa) (a) particle resolution is 1 mm; (b) particle resolution is 1.25 mm; (c) particle resolution is 1.5 mm

900 800

Experimental data GPD damage model

1200

GPD damage model Expermental data

1000

700

48% of σ1, peak

800

σ1-σ3/ MPa

σ1-σ3 / MPa

600 500 400 300 200

58% of σ1, πεακ

600 400 200

100 0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0 0.0

1.6

0.2

0.4

0.6

(a)

σ 2 = σ 3 = 100 MPa.

(c) GPD damage model Experimental data

1.0

1.2

1.4

1.6

σ 2 = 440 MPa, σ 3 = 100 MPa. GPD damage model Experimental data

1200

1000

1000

800

800

50% of σ1, peak

σ1-σ3/ MPa

σ1-σ3/ MPa

1200

600 400 200 0 0.0

0.8

ε1/ %

ε1/ %

60% of σ1, peak

600 400 200

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0 0.0

1.6

0.2

0.4

0.6

ε1/ %

(b)

0.8

1.0

1.2

1.4

1.6

ε 1/ %

σ 2 = 300 MPa, σ 3 = 100 MPa.

(d)

σ 2 = 600 MPa , σ 3 = 100 MPa.

Fig. 15. The overall stress (σ1 − σ3) vs strain obtained from the numerical simulation and experiment data with σ3 = 100 MPa (Haimson and Chang, 2000): (a) σ2 = 100 MPa, (b) σ2 = 300 MPa (c) σ2 = 440 MPa and (d) σ2 = 600 MPa.

reached its pre-determined value (i.e.160 MPa). Then, σ1 and σ3 were kept constant, σ2 was increased until it reached its pre-determined value (i.e. 320 MPa). Thereafter, σ1 was increased in a way that the axial loading velocity of 0.05 m/s and a time step of 5 × 10− 6 s are adopted. The complete stress-strain curve of granodiorite obtained from experimental data is depicted in Fig. 6 (a). The complete stress-strain curve of granodiorite sample obtained from the proposed numerical

method is plotted in Fig. 6 (b). Compared the numerical stress-strain curves with the experimental stress-strain curves of granodiorite, the numerical stress-strain curves of granodiorite sample are in good agreement with the experimental stress-strain curve of granodiorite sample. Moreover, the numerical strength of granodiorite is equal to 1040 MPa, which is in good agreement with 1050 MPa obtained from the experiments. 36

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

900

Table 2 The true triaxial stress conditions and the numerical results. Numerical maximum principal stress σ1/MPa

0 0 0 0 0 25 25 25 25 25 65 65 65 65 65 65 105 105 105 105 105

0 40 80 150 200 30 90 200 300 400 50 100 200 300 400 500 100 200 300 400 500

163.4 195.8 220.3 216.8 216.7 296.0 338.4 399.0 450.1 400.9 435.1 465.6 524.2 575.1 618.0 622.2 682.8 782.1 859.6 877.5 850.6

σ3=105MPa

σ2=σ3 750

σ3=65MPa 600 /MPa

σ2/MPa

1

σ3/MPa

450

σ3=25MPa

300

σ3=0MPa Dunham dolomite

150 0

100

200

300

400

500

600

700

σ2/MPa

(a) The numerical results. Fig. 7 shows the distribution of the maximum principal stress of the 3D sample at different time steps as well as their damage evolution. In addition, markers A, B, C, D, E and F on the stress-strain curve in Fig. 6(b) correspond to the different damaged stages in Fig. 7. Fig. 7 (a) shows the stress distribution at t = 300 μs, which corresponds to the marker “A” on the stress-strain curve in Fig. 6 (b). At this stage, uniform damage initiates, and stress concentration occurs around the damage points. Fig. 7 (b) shows the stress distribution at t = 500 μs, which corresponds to the marker “B” on the stress-strain curve in Fig. 6 (b). It is concluded from Fig. 7 (b) that a few of damage points coalesce with each other, localized damage initiates, and the obvious stress concentration around the localized damage zone occurs. Fig. 7 (c) shows the stress distribution at t = 700 μs, which corresponds to the marker “C” on the stress-strain curve in Fig. 6 (b). It is found from Fig. 7 (c) that multi-localized-damage zones initiate, the obvious stress concentration around the multi-localized-damage zones occurs, which leads to the rapid growth of the localized damage zone. Fig. 7 (d) shows the stress distributions at t = 900 μs, which corresponds to the marker “D” on the stress-strain curve in Fig. 6 (b). In Fig. 7 (d), multi-localized-damage zones coalesce with each other, which leads to the occurrence of the damage localization band, and the load-bearing capacity of granodiorite sample reaches the peak value. Fig. 7 (e) shows the stress distributions at t = 1000 μs, which corresponds to the marker “E” on the stress-strain curve in Fig. 6 (b). In Fig. 7 (e), the damage localization band passes through the whole granite specimen, which leads to the instability of granodiorite specimen and decrease of the load-bearing capacity of granodiorite sample. Fig. 7 (f) shows the stress distributions at t = 1200 μs, which corresponds to the marker “F” on the stress-strain curve in Fig. 6 (b). It is observed from Fig. 7 (f) that the failure of granodiorite sample occurs, and the load-bearing capacity of granodiorite sample rapidly drops to the point “F”, which corresponds to the residual strength. A dry granodiorite sample subjected to the true triaxial test exhibits a through-going single shear fracture steeply, as shown Fig. 8 (a). The scatter electron micrographs (SEM) was adopted to observe the damage localization zone, as shown in Fig. 8 (b). It is found from Fig. 8 (b) that the vertical microcrack density is the highest adjacent to the main fault. We infer that the localized micro-cracks are coalesced to form damage localization and ultimately result in through-going main fault as σ1 approaches its peak value. General micromechanical features described here are analogous to what observed in the present numerical simula-

900

σ3=105MPa

σ1=σ2

750

σ3=65MPa

σ1/MPa

600

σ3=25MPa

450 300

σ3=0MPa 150

Dunham dolomite 0

0

100

200

300

400

σ2/MPa

(b) The experimental data. Fig. 16. A comparison between the numerical results and experiment data for the compressive strength of Dunham dolomite under the true triaxial compression: (a) the numerical results (b) the comprehensive true triaxial tests conducted by Mogi (1971). Table 3 The main material parameters of granite.

37

Rock

Density (kg/m3)

Uniaxial compressive strength (MPa)

Young's modulus (GPa)

Poisson's ratio

Intermediate stress parameter

Granite

3200.0

190.0

50

0.25

1.0

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

uniform damage points. In the third stage, the maximum principal stress non-uniformly distributes, and the maximum principal stress concentrates around localized damage points. In the fourth stage, the maximum principal stress concentrates around the localized damage band, and the maximum principal stress around other damage points reduce. In the fifth stage, the concentration of the maximum principal stress reduces around the localized damage as well as other damage points. Fig. 11 shows the macroscopic failure mode of rock samples. Fig. 11 (a) shows the X-shaped conjugate shear damage localization bands. Fig. 11 (b) shows the axial splitting damage localization mode. Fig. 11 (c) shows cone-shaped damage localization mode. The numerical results are in good agreement with the experimental observations (Chen and Feng, 2006), as shown in Fig. 12.

Table 4 The true triaxial stress conditions and the numerical results. σ3/MPa

σ2/MPa

Numerical maximum principal stress σ1/MPa

100 100 100 100 100 100 100

100 300 440 600 800 1000 1200

742.7 798.8 900.3 996.8 1086.7 1174.0 1131.4

4.3. Effects of different particle resolutions on the damage localization and failure of rocks The application of GPD to analyze the damage localization and failure of rocks is relatively new. Hence, it is instructive to perform accuracy analysis by different particle resolutions. The results presented in this paper contain three different particle resolution of 1 mm, 1.25 mm and 1.5 mm in the GPD model for σ1 of 190 MPa, σ2 of 50 MPa and σ3 of 10 MPa. The numerical results for three different particle resolutions are plotted in Figs. 13 and 14. It can be found from Fig. 13 that the stress-strain curves are seen to be same initially and obvious bifurcation arises in the end of curves, and the accuracy for the particle resolution of 1 mm is the highest, while the accuracy for the particle resolution of 1.5 mm is the lowest. Moreover, Fig. 13 shows that the numerical result for the particle resolution of 1.25 mm is very closed to that for the particle resolution of 1.0 mm. That is to say, numerical solutions for the particle resolution of 1.0 mm are not obviously mesh-dependent Therefore, in this paper, the particle resolution of 1.0 mm is adopted. Fig. 14 shows the damaged localization for different particle resolutions. It is observed form Fig. 14 that the failure modes slightly depend on GPD particle resolutions, and the characteristics of damaged localization is more and more obvious as the particle resolution decreases.

Fig. 17. The relationship between σ1 and the intermediate principal stress at different intermediate stress parameter b when σ3 = 100 MPa.

tion under the same stress condition, as shown in Fig. 8 (c). It is observed from Fig. 8 that the damaged profile obtained from the proposed numerical method is in good agreement with the experimental observation (Lee and Haimson, 2011). 4.2. Damage localization modes of the 3D specimens subjected to different stress conditions In this section, the proposed GPD is applied to describe the damage localization modes of brittle rocks subjected to different triaxial loads. In additional, the ultimate strength properties are also studied. For this purpose, a numerical sample with the size of 60 mm × 20 mm × 60 mm is adopted. The representative mechanical parameters of sample in this section are the same as those described in Table 1 in Section 4.1. The value of b is taken 1.0. In the numerical simulation tests, σ2 and σ3 are kept constant, which are respectively equal to 15 MPa and 10 MPa, 50 MPa and 10 MPa, 80 MPa and 30 MPa. The value of σ1 is raised up to the failure of sample at a constant displacement control rate of 0.15 × 10− 3 mm/s. The numerical axial stress-axial strain curves are plotted in Fig. 9. The distribution of the maximum principal stress at different time steps is depicted in Fig. 10. Fig. 9 shows the numerical stress-strain curves under different compressive stress. Generally, four stages are observed during the whole process of compression. In the first stage, very few damage points occur, and rock material exhibits elastic behaviors before point A. In the second stage, the uniform damage forms from point A to point B, and rock material exhibits nonlinear hardening behaviors. In the third stage, localized damage occurs from point B to point C, and the load-bearing capacity of rock material reaches the peak value at point C. In the fourth stage, damage localization band develops from point C to point D, the rock material exhibits the strain softening behaviors, and the load-bearing capacity of rock material reduces. In the fifth stage, damage localization band forms, the failure of the rock material occurs. The distribution of the maximum principal stress σ1 at different time-steps is plotted in Fig. 10. In the first stage, the maximum principal stress σ1 uniformly distributes due to occurrence of no damage. In the second stage, the maximum principal stress concentrates around the

5. Effects of the intermediate principal stress 5.1. The effect of the intermediate principal stress on the stress-strain curve of rocks The sample size and loading rate and the mechanical parameters of granite in the numerical simulation are similar to those used in the true triaxial tests (Haimson and Chang, 2000). The true triaxial test data show that b = 1.0 for granite(Yu et al., 2002). In the true triaxial tests, σ3 was equal to 100 MPa, and σ2 was varied from 100 MPa to 600 MPa. Fig. 15 shows the numerical results of relationship between the deviator stress (σ1 − σ3) and the axial strain (ε1). In Fig. 15 (a), when the maximum principal stress σ1 reaches 48% of the peak value, the first damage point occurs. In Fig. 15 (b), when the maximum principal stress σ1 reaches 50% of the peak value, the first damage point occurs. In Fig. 15 (c), when the maximum principal stress σ1 reaches 58% of the peak value, the first damage point occurs. In Fig. 15 (d), when the maximum principal stress σ1 reaches 60% of the peak value, the first damage point occurs. That is to say, the ratio of the maximum principal stress at the first damage point to the peak strength increases with increasing the intermediate principal stress. Moreover, the experimental stress-strain curves of samples were also plotted in Fig. 15 (Haimson and Chang, 2000). Comparing the numerical results with experimental data, it is found that the numerical results are in good agreement with the experiment data. In addition, a similar observation in sandstone was made by Takahashi and Koide (1989). 38

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 18. The maximum principal stress contour of samples at different intermediate stress parameter (Unit: Pa) (a) σ2 = σ3 = 100MPa; (b) σ2 = 300MPa; σ3 = 100MPa; (c) σ2 = 440MPa; σ3 = 100MPa; (d) σ2 = 600MPa; σ3 = 100MPa.

39

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 18. (continued)

Dunham dolomite. For this purpose, the specimens subjected to the true triaxial numerical test are made on Dunham dolomite, whose mechanical parameters are the same as those given by Mogi (1967a,b). The true triaxial stress conditions and the numerical results are listed in Table 2.

5.2. Effects of the intermediate principal stress on the strength of rocks In this section, the proposed numerical model is implemented to study the typical effects of the intermediate stresses on the strength of 40

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

Fig. 19. The test results of the failure angle with different σ2 (Feng et al., 2016).

Fig. 20. The numerical results of the failure angle with different σ2.

110

5.3. Effects of the intermediate stress parameter on the strength and damaged localization of rocks

Numerical results Experiment data

105

Failure angle θ/°

100

The representative mechanical parameters of sample in this section are shown in Table 3. The true triaxial stress conditions and the numerical results are listed in Table 4. Fig. 17 shows that the effect of the intermediate principal stress on the peak strength at different intermediate stress parameter b. It is found from Fig. 17 that the peak strength increases with increasing the intermediate stress parameter b, and the peak strength first increases, then decreases with an increase in the intermediate principal stress. Fig. 18 shows that the effect of the intermediate principal stress on the maximum principal stress distribution and the damaged localization of samples at different intermediate stress parameter b. It is found from Fig. 18 that the maximum principal stress concentration around the damage localization zone is more and more obvious with increasing the intermediate principal stress, the failure of rock samples varies from splitting failure to shear failure as the intermediate principal stress increases, and the magnitude of the damage localization zone decreases with increasing the intermediate stress parameter.

95 90 85 80 75 70 65

100

200

300

400

500

Intermediate principal stress σ 2 /MPa Fig. 21. The failure angle θ at the different value of σ2.

The true triaxial test data show that b = 1.0 for Dunham dolomite (Yu et al., 2002). The effect of the intermediate stress on the strength of Dunham dolomite is plotted in Fig. . 16. In Fig. . 16, the solid line represents the conventional triaxial strength criterion for σ3 = σ2, dashed lines show the true triaxial strength trend for each set of tests for which σ3 is kept constant. Fig. . 16 shows that the strength of Dunham dolomite increases with increasing the minimum principal stress, the effect of the intermediate principal stress has two zones: (1) rock strength increases with increasing the intermediate principal stress from σ2 = σ3 to a maximum value (reaches the peak strength when σ2 = σ2′); (2) rock strength decreases with the further increase of the intermediate principal stress from σ2 = σ2′ to σ2 = σ1. The numerical results are in good agreement with the experiment data (Mogi, 1967a,b). Remarkably, this similar behavior is also predicted by the theoretical model proposed by Wiebols and Cook (1968).

5.4. The effects of the intermediate principal stress on the failure angle In order to study the effect of intermediate principal stress on the failure angle, granite sample subjected to the true triaxial compression is adopted by Feng et al. (2016). The strike and dip of the failure are shown in Fig. 19. Fig. 20 shows that the numerical results of the failure angle of rock samples when σ3 is fixed. In Figs. 20 and 21, θ is the angle between the damage localization band and the minimum principal stress direction, which is called the failure angle. It is observed from Fig. 21 that the failure angle θ increases when the intermediate principal stress varies from 100 MPa to 300 MPa, then the failure angle θ decreases when the intermediate principal stress is > 300 MPa. The numerical results 41

Engineering Geology 224 (2017) 29–42

X.P. Zhou, X.Q. Zhang

localization in cemented sands by two-dimensional distinct element method analyses. Comput. Geotech. 38 (1), 14–29. Klein, E., Baud, P., Reuschl, T., Wong, T.F., 2001. Mechanical behaviour and failure mode of bentheim sandstone under triaxial compression. Phys. Chem. Earth A 26 (1–2), 21–25. Leblanc, S., Boyer, P., Joslin, C., 2014. Modelling and animation of impact and damage with smoothed particle hydrodynamics. Vis. Comput. 30 (6), 909–917. Lee, H., Haimson, B.C., 2011. True triaxial strength, deformability, and brittle failure of granodiorite fromthe San Andreas Fault Observatory at Depth. Int. J. Rock Mech. Min. Sci. 48 (7), 1199–1207. Mogi, K., 1967a. Effect of the intermediate principal stress on rock failure. J. Geophys. Res. Atmos. 72 (20), 5117–5131. Mogi, K., 1967b. Effect of the triaxial stress system on rock failure. Rock Mech. Jpn. 1, 53–55. Pan, P.Z., Feng, X.T., Hudson, J.A., 2012. The influence of the intermediate principal stress on rock failure behaviour: a numerical study. Eng. Geol. 124 (4), 109–118. Qi, M., Shao, J.F., Giraud, A., Zhu, Q.Z., Colliat, J.B., 2016. Damage and plastic friction in initially anisotropic quasi brittle materials. Int. J. Plast. 82, 260–282. Shao, J.F., Chau, K.T., Feng, X.T., 2006. Modeling of anisotropic damage and creep deformation in brittle rocks. Int. J. Rock Mech. Min. Sci. 43 (4), 582–592. Takahashi, M., Koide, H., 1989. Effect of the intermediate principal stress on strength and deformation behavior of sedimentary rocks at the depth shallower than 200m. Rock Great Depth 8, 28–31. Vu-Bac, N., Nguyen-Xuan, H., Chen, L., Lee, C.K., Zi, G., Zhuang, X., Liu, G.R., Rabczuk, T., 2013. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. J. Appl. Math. 2013 (1), 1–12. Weibull, W., 1951. A statistical distribution function of wide applicability. J. Appl. Mech 18 (2), 293–297. Wiebols, G.A., Cook, N.G.W., 1968. An energy criterion for the strength of rock in polyaxial compression. Int. J. Rock Mech. Min. Sci. Geomech. 5 (6), 529–549. Wong, R.H., Tang, C.A., Chau, K.T., Lin, P., 2006. Splitting failure in brittle rocks containing pre-existing flaws under uniaxial compression. Eng. Fract. Mech. 69 (17), 1853–1871. Wu, W., Zhu, H., Zhuang, X., Ma, G., Cai, Y., 2014. A multi-shell cover algorithm for contact detection in the three dimensional discontinuous deformation analysis. Theor. Appl. Fract. Mech. 72 (1), 136–149. Wu, J.Y., Li, F.B., Xu, S.L., 2015. Extended embedded finite elements with continuous displacement jumps for the modeling of localized failure in solids. Comput. Methods Appl. Mech. Eng. 285 (3), 346–378. Yang, S.Q., Jing, H.W., 2013. Evaluation on strength and deformation behavior of red sand-stone under simple and complex loading paths. Eng. Geol. 164 (393), 1–17. Yao, C., Jiang, Q.H., Shao, J.F., Zhou, C.B., 2016. A discrete approach for modeling damage and failure in anisotropic cohesive brittle materials. Eng. Fract. Mech. 155, 102–118. Yu, M.H., 1983. Twin shear stress yield criterion. Int. J. Mech. Sci. 25 (1), 71–74. Yu, M.H., Zan, Y.W., Zhao, J., Yoshimine, M., 2002. A unified strength criterion for rock material. Int. J. Rock Mech. Min. Sci. 39 (8), 975–989. Zhou, X.P., 2005. Localization of deformation and stress-strain relation for mesoscopic heterogeneous brittle rock materials under unloading. Theor. Appl. Fract. Mech. 44 (1), 27–43. Zhou, X.P., Zhao, Y., Qian, Q.H., 2015a. A novel meshless numerical method for modeling progressive failure processes of slopes. Eng. Geol. 192 (46), 139–153. Zhou, X.P., Bi, J., Qian, Q.H., 2015b. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech. Rock. Eng. 48 (3), 1097–1114. Zhu, Q.Z., Shao, J.F., 2015. A refined micromechanical damage-friction model with strength prediction for rock-like materials under compression. Int. J. Solids Struct. 60–61 (5), 75–83. Zhu, H.H., Zhuang, X.Y., Cai, Y.C., Ma, G.W., 2011. High rock slope stability analysis using the enriched meshless shepard and least squares method. Int. J. Comput. Methods 8 (2), 209–228. Zhu, Q.Z., Zhao, L.Y., Shao, J.F., 2016. Analytical and numerical analysis of frictional damage in quasi brittle materials. J. Mech. Phys. Solids 92, 137–163.

obtained from GPD are in good agreement with the previous experimental observations (Feng et al., 2016). 6. Conclusion In this paper, General Particle Dynamics (GPD) is developed to model the localized damage and failure of 3D rock sample under the true triaxial compression condition (TTC). The conclusions are summarized as follows: a realistic evaluation of the damage localization process of a cubical specimen subjected to the true triaxial compressive loads is proposed. The damage failure mechanisms are revealed using the proposed damage model. The effect of the intermediate principal stress on the peak strength, the stress-strain curve, the damage localization mode and the failure angle of rocks is investigated. It is found that the numerical results are in good agreement with experimental observations. In brief, the proposed GPD method can be applied to simulate the damage localization and predict the failure of brittle rock materials. Acknowledgment This work was supported by Project 973 (Grant no. 2014CB046903), the National Natural Science Foundation of China (Nos. 51325903 and 51279218), Natural Science Foundation Project of CQ CSTC (Nos CSTC, cstc2015jcyjys30001, cstc2016jcyjys0005 and cstc2015jcyjys30006), Graduate Scientific Research and Innovation foundation of Chongqing, China (Grant No. CYB16012). References Bésuelle, P., Desrues, J., Raynaud, S., 2000. Experimental characterization of the localization phenomenon inside a Vosges sandstone in a triaxial cell. Int. J. Rock Mech. Min. Sci. 37 (8), 1223–1237. Cai, Y., Zhu, H., Zhuang, X., 2014. A continuous/discontinuous deformation analysis (CDDA) method based on deformable blocks for fracture modelling. Front. Struct. Civ. Eng. 7 (4), 369–378. Chen, J., Feng, X.T., 2006. Ture triaxial experimental study on rock with high geo-stress. J. Rock Mech. Eng. 25 (8), 1538–1539. Fang, Z., Harrison, J., 2002. Development of a local degradation approach to the modelling of brittle fracture in heterogeneous rocks. Int. J. Rock Mech. Min. Sci. 39, 443–457. Feng, X.T., Zhang, X., Kong, R., Wang, G., 2016. A novel mogi type true triaxial testing apparatus and its use to obtain complete stress-strain curves of hard rocks. Rock Mech. Rock. Eng. 49 (5), 1649–1662. Fernández-Méndez, S., Bonet, J., Huerta, A., 2005. Continuous blending of SPH and finite elements. Comput. Struct. 83 (17–18), 1448–1458. Haimson, B., Chang, C., 2000. A new true triaxial cell for testing mechanical properties of rock and its use to determine rock strength and deformability of Westerly granite. Int. J. Rock Mech. Min. Sci. 37 (1), 285–296. Hoek, E., 1990. Estimating Mohr-Coulomb friction and cohesion valuesfrom the HoekBrown failure criterion. Int. J. Rock Mech. Min. Sci. 27 (3), 227–229. Hu, D.W., Zhu, Q.Z., Zhou, H., Shao, J.F., 2010. A discrete approach for anisotropic plasticity and damage in semi-brittle rocks. Comput. Geotech. 37 (3), 658–666. Jiang, M.J., Yan, H.B., Zhu, H.H., Utili, S., 2011. Modeling shear behavior and strain

42