# Failure Process Simulation of Interlayered Rocks under Compression

Authors: Chi Yao; Sizhi Zeng; Jianhua Yang
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Anisotropy in strength and deformation of rock mass induced by bedding planes and interlayered structures is a vital problem in rock mechanics and rock engineering. The modified rigid block spring method (RBSM), initially proposed for modeling of isotropic rock, is extended to study the failure process of interlayered rocks under compression with different confining pressures. The modified rigid block spring method is used to simulate the initiation and propagation of microcracks. The Mohr–Coulomb criterion is employed to determine shear failure events and the tensile strength criterion for tensile failure events. Rock materials are replaced by an assembly of Voronoi-based polygonal blocks. To explicitly simulate structural planes and for automatic mesh generation, a multistep point insertion procedure is proposed. A typical experiment on interlayered rocks in literature is simulated using the proposed model. Effects of the orientation of bedding planes with regard to the loading direction on the failure mechanism and strength anisotropy are emphasized. Results indicate that the modified RBSM model succeeds in capturing main failure mechanisms and strength anisotropy induced by interlayered structures and different confining pressures.


---

## Body

## 1. Introduction

Foliated metamorphic rocks and sedimentary rocks are ubiquitous in nature. Due to the presence of weak structural planes or interlayered structures, deformation and strength of this kind of rocks are often viewed as transversely isotropic. It is a major concern that the compressive strength and failure mode of such rocks vary as the change of confining pressures and loading orientations with regard to the structural planes. Here, for simplicity, let the angle between the first principal stress and structural planes beα, as is shown in Figure 1.Figure 1
Definition of structural plane orientation.After years of research on anisotropic rocks, researchers found that the maximum strength took place whenα = 0° or 90° for most rocks of this type, while the minimum strength took place when α was between 30° and 45° [1–6]. Curves of strength versus α can be categorized into three types [7] (Figure 2), namely, U type, undulatory type, and shoulder type, among which, the third one is the most common type. Scholars have proposed a number of failure criteria to describe such anisotropy on strength. Jeager [8] viewed rock matrix as isotropic, assumed that rock failure was due to the frictional sliding along structural planes, and proposed a failure criterion based on the Mohr–Coulomb criterion. But Jeager’s criteria could not take into consideration the case that strength for the orientation α = 0° is not equal to which for the orientation α = 90°, which is often encountered in laboratory experiments. Bearing this in mind, other researchers extended these criteria, and McLamore et al. [9] added two parameters to account for this strength discrepancy. Duveau et al. [10] provided another modification by replacing the Mohr–Coulomb criterion with a nonlinear Barton model. Tien and Kuo [11] introduced the maximum axial strain theory into Jaeger’s criterion and succeeded in modeling various types of transversely isotropic rocks.Figure 2
Classification of anisotropy for transversely isotropic rocks (after Ramamurthy [7]): (a) U type; (b) undulatory type; (c) shoulder type.
(a)
(b)
(c)These criteria play quite well in predicting strength anisotropy quantitatively, but they fail to describe propagation processes of microcracks and thus cannot produce failure modes in an explicit manner. Moreover, quite a lot of parameters are often involved in such criteria. As a result, large amount of tests are needed for parameter calibration.In recent years, great progress has been made in numerical modeling of rock fracturing, and some groups tried to simulate the failure process of transversely isotropic rocks at mesoscale in an explicit manner. Debecker and Vervoort [12] conducted a series of numerical tests on behaviors of slate under uniaxial compression and Brazilian disc test using UDEC. By combining the parallel bonded model and the smooth joint model, Chiu et al. [13] modeled the failure process of a rock containing persistent joints under compression with different confining pressures with PFC3D. Wang et al. [14] analyzed the Brazilian tensile splitting tests of stratified biotite granulite rocks with PFC2D. The computed stress-strain curves were compared well with experimental results. Shen et al. [15] simulated a layered rock under triaxial compression using a finite element method. Lisjak et al. [16] simulated the anisotropic behaviors of Opalinus argillite by the combined FEM/DEM model. Zhao et al. [17–20], Wang et al. [21], and Wei et al. [22] have also made great contributions to this field.The rigid block spring method (RBSM), which was initially proposed by Kawai [23], provides an alternative discrete approach for modeling crack growth and fracturing in cohesive brittle materials. In this method, the cohesive material is replaced by an assemblage of polygonal discrete elements (rigid blocks) by using the Voronoi diagram. The blocks are bonded by contact interfaces. Three springs are defined on each interface to describe the relative normal and tangential displacements and rotations between two neighboring blocks. The initial RBSM method was so far successfully applied to modeling macroscopic behaviors of cement-based materials [24, 25]. This method has been successively improved by Qian and Zhang [26]. In the modified RBSM model, a continuous distribution of stress and relative displacements is defined on each interface through suitable interpolation functions. This allows to modeling the progressive failure of interfaces. The modified method was recently successfully applied to modeling the damage and failure of brittle rocks [27–31]. The advantages of this method have been discussed. This method is well applicable for small-deformation and quasi-static problems. Since the RBSM is an implicit discrete method without updating contact distribution, its implementation is easy and it provides a fast convergence. And it is not necessary to introduce an artificial damping coefficient.However, the modified RBSM model has been so far applied to initially isotropic materials. The objective of the present work was to extend this method to rock-like materials exhibiting initial inherent anisotropy. Both tensile and shear failure of interfaces will be considered. The proposed anisotropic RBSM model will be applied to modeling the mechanical behavior of typical artificial interlayered rocks, which are carried out by Tien et al. [6]. In particular, influences of the initial structural anisotropy on the induced damage and failure process will be discussed.


## 2. Principles of the Modified Rigid Block Method

The governing equations for the modified RBSM are derived from the Virtual Work Theorem, which is simply given below. Suppose that two arbitrary blocks, block 1 and block 2, are next to each other. The centroid of the two blocks are (x1, y1) and (x2, y2), respectively. Point P1 on block 1 and point P2 on block 2 coincide with the same point (x, y) on the interface between the two blocks, as can be seen in Figure 3. Assuming rotations are small, the relative displacements between point P1 and point P2, Δu, can be expressed by the displacements of centroids of the two blocks U12 as(1)Δu=BNU12,where(2)Δu=ΔunΔusT,in which Δun and Δus are, respectively, the relative normal and tangential displacements between point P1 and point P2.(3)B=lm−ml,where l,m is the unit normal vector of the interface.(4)N=−10y1010−y200−1−x1001x20,where x10=x−x1, y10=y−y1, x20=x−x2, and y20=y−y2.(5)U12=U1xU1yU1θU2xU2yU2θT,where U1x, U1y, and U1θ are, respectively, the translational displacement in the x direction, y direction, and the rotational angle of the centroid of block 1, while U2x, U2y, and U2θ are displacement components for block 2.Figure 3
Local deformation of an interface.Stress induced by relative displacements between pointP1 and point P2 can be expressed as(6)σ=DΔu,where(7)σ=σnσsT,in which σn and σs are normal stress and tangential stress.(8)D=kn00ks,where kn is the stiffness of the normal spring and ks is the stiffness of the tangential spring.According to the Virtual Work Theorem, for the block system, the relationship is as follows:(9)∑e∫l0eδΔuTσdl+∫lσeδuTpdl+∫∫seδuTfds=0,where l0e, lσe, and se, respectively, stand for interfaces between blocks, force boundaries, and domains of blocks; p, f, δu, and δΔu are, respectively, external pressure, body force, virtual displacement, and virtual relative displacement.Together with above relations, the global equilibrium equation for the whole block system can be obtained:(10)KU=F.A linear elastic relation is employed to describe the stress and displacement relation between two neighboring blocks. The normal stress is calculated through local constitutive law shown in Figure4 and is composed of two parts, the compressive and the tensile components. In this chapter, tension is defined as positive.Figure 4
Normal stress displacement relationship on interfaces.In compression,σn is given as(11)σn=knΔun,where Δun is the relative displacement between interacting blocks and kn is the normal stiffness.In tension, the normal stress is still computed with the same stiffness in compression. The maximum tensile stressσn max is equal to the tensile strength T, such that(12)σnmax=T.After the maximum tensile stress is reached, the normal stress is set to be zero.Due to the possible change in orientation during iteration, the shear stress is computed incrementally and is defined as(13)σs=σsupdated+ksΔus,where σsupdated is the updated shear stress, ks is the tangential stiffness, proportional to kn, and Δus is increment relative tangential displacement.Uniformly distributed random Voronoi diagram is used as mesh to discretize the concerning domain of targeted rocks. With this kind of mesh,kn and ks can be determined by macroelastic parameters, that is, elastic modulus E and Poisson’s ratio v, according to Yao et al. [27]:(14)r=ks/kn=4.025ν4−6.087ν3+6.022ν2−3.966v+1,E/E0=−0.6291r4+1.617r3−1.678r2+1.174r+0.5162,kn=E0h1+h2,ks=r⋅kn,where r and E0 are intermediate variables and h1 and h2 denote the distances from the centroids of two neighboring blocks to their connecting interface.To simulate the fracturing and damage process of geo materials, a modified Mohr–Coulomb model is adopted, as shown in Figure5. The maximum allowable shear stress is computed by the normal stress σn, the cohesion C, the critical normal stress σn critical, the local frictional angle φ1, and the local residual frictional angle φ2. The critical normal stress σn critical is defined to limit the frictional strengthening effects. Before rupture, the maximum shear stress is computed by(15)σsmax⁡=σntan⁡φ1+C,σn>σn critical,σsmax⁡=σntan⁡φ1+C,σn≤σn critical.Figure 5
Failure criterion of interfaces.Shear rupture takes place whenσs>σsmax, and then the interaction becomes purely frictional, with a maximum shear force defined by(16)σsmax⁡=σntan⁡φ2,σn>σn critical,σsmax⁡=σntan⁡φ2,σn≤σn critical.In the modified rigid block spring method, structural planes and virtual cracks are treated with the same failure criterion but different failure parameters.


## 3. Mesh Generation

Two main steps are taken to generate Voronoi diagram-based mesh: point insertion and tessellation. Each step is described in detail as follows.


### 3.1. Point Insertion

The process of point insertion is based on the concept of point saturation. Point saturation is achieved by maintaining distance between neighboring points under a minimum admissible distancelmin. To explicitly model bedding planes, a multistep insertion procedure is adopted. Here, we define a segment of boundaries or bedding planes as a segment and an intersection between boundaries or bedding planes as a vertex.


#### 3.1.1. Insertion of Points to Define Vertexes

Suppose there is a vertexV connected by four segments, as shown in Figure 6. To define this vertex, firstly, find out the minimum angle between two neighboring segments, the value of which is assumed to be 2β. Then, draw auxiliary lines which have an angle of β with respect to each segment. After that, draw a circle with a radius of 0.5lmin around vertex V. Insert points at intersections between the circle and all auxiliary lines around vertex V. The vertex V is defined by all the inserted points. An example is also given for illustration of the process of mesh generation, which is shown in Figure 7. Shown in Figure 7(b) is the step of definition of all vertexes.Figure 6
Illustration of vertex definition.Figure 7
Illustration of mesh generation process.
(a)
(b)
(c)
(d)
(e)
(f)


#### 3.1.2. Insertion of Points to Define Segments

Pairs of points are symmetrically inserted to define each segment. Distances between each neighboring points on the same side of preexisting line are larger thanlmin. Shown in Figure 7(c) is the step of segment definition.


#### 3.1.3. Insertion of Points to the Whole Domain

With the constraint oflmin, points are sequentially inserted to saturate the whole domain, as shown in Figure 7(d).


### 3.2. Tessellation

There are various methods available for Voronoi tessellation with a set of points. In the present study, a sweep line algorithm proposed by Fortune [32] is used. With points inserted in the previous step, a Voronoi diagram is generated, as shown in Figure 7(e). After removing line segments outside the boundaries, the mesh we need for simulation is obtained (Figure 7(f)).


## 3.1. Point Insertion

The process of point insertion is based on the concept of point saturation. Point saturation is achieved by maintaining distance between neighboring points under a minimum admissible distancelmin. To explicitly model bedding planes, a multistep insertion procedure is adopted. Here, we define a segment of boundaries or bedding planes as a segment and an intersection between boundaries or bedding planes as a vertex.


### 3.1.1. Insertion of Points to Define Vertexes

Suppose there is a vertexV connected by four segments, as shown in Figure 6. To define this vertex, firstly, find out the minimum angle between two neighboring segments, the value of which is assumed to be 2β. Then, draw auxiliary lines which have an angle of β with respect to each segment. After that, draw a circle with a radius of 0.5lmin around vertex V. Insert points at intersections between the circle and all auxiliary lines around vertex V. The vertex V is defined by all the inserted points. An example is also given for illustration of the process of mesh generation, which is shown in Figure 7. Shown in Figure 7(b) is the step of definition of all vertexes.Figure 6
Illustration of vertex definition.Figure 7
Illustration of mesh generation process.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.1.2. Insertion of Points to Define Segments

Pairs of points are symmetrically inserted to define each segment. Distances between each neighboring points on the same side of preexisting line are larger thanlmin. Shown in Figure 7(c) is the step of segment definition.


### 3.1.3. Insertion of Points to the Whole Domain

With the constraint oflmin, points are sequentially inserted to saturate the whole domain, as shown in Figure 7(d).


## 3.1.1. Insertion of Points to Define Vertexes

Suppose there is a vertexV connected by four segments, as shown in Figure 6. To define this vertex, firstly, find out the minimum angle between two neighboring segments, the value of which is assumed to be 2β. Then, draw auxiliary lines which have an angle of β with respect to each segment. After that, draw a circle with a radius of 0.5lmin around vertex V. Insert points at intersections between the circle and all auxiliary lines around vertex V. The vertex V is defined by all the inserted points. An example is also given for illustration of the process of mesh generation, which is shown in Figure 7. Shown in Figure 7(b) is the step of definition of all vertexes.Figure 6
Illustration of vertex definition.Figure 7
Illustration of mesh generation process.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1.2. Insertion of Points to Define Segments

Pairs of points are symmetrically inserted to define each segment. Distances between each neighboring points on the same side of preexisting line are larger thanlmin. Shown in Figure 7(c) is the step of segment definition.


## 3.1.3. Insertion of Points to the Whole Domain

With the constraint oflmin, points are sequentially inserted to saturate the whole domain, as shown in Figure 7(d).


## 3.2. Tessellation

There are various methods available for Voronoi tessellation with a set of points. In the present study, a sweep line algorithm proposed by Fortune [32] is used. With points inserted in the previous step, a Voronoi diagram is generated, as shown in Figure 7(e). After removing line segments outside the boundaries, the mesh we need for simulation is obtained (Figure 7(f)).


## 4. Numerical Example

Using artificial rock like materials, Tien et al. [6] conducted a series of tests in order to study failure modes and strength anisotropy of interlayered rocks under compression with different confining pressures. The artificial rock is formed from stratified layers made of material A and material B. Thicknesses of layers for the two materials are the same. But mechanical properties of them are quite different. Model material A was composed of kaolinite, cement, and water in ratio of 4 : 1 : 1.2, while B in ratios of 1 : 1 : 0.6, and mechanical parameters of these two materials are listed in Table 1.Table 1
Macromechanical parameters of materials in Tien et al.’s experiment [6].
E (GPa)
v
UCS (MPa)
Frictional angle (°)
Cohesion (MPa)
Material A
21.7
0.23
104.2
29
30.7
Material B
11.9
0.21
43.3
25
13.8
UCS: uniaxial compressive strength.


### 4.1. Calibration of Input Parameters for Material A and Material B

Input parameters of the two distinct materials are calibrated independently. They are identified through a trial and error optimization algorithm against the strength envelops of conventional triaxial compression test. In each step of calibration, simulation results are compared with experimental data at all confining pressures. Calibration progresses until simulation results match well with experimental data on the whole. Shown in Tables2 and 3 are, respectively, trial parameters for calibration of material A and material B. The corresponding numerical strength envelops comparing with experimental data are shown in Figures 8 and 9. The calibration process is through adjusting input parameters to gradually approach the experimental data. The finally obtained input parameters are listed in Table 4.Table 2
Trial input parameters for calibration of material A.
E (GPa)
v
tan
φ
1
tan
φ
2
C (MPa)
T (MPa)
σ
n
critical (MPa)
Trial 1
21.7
0.23
0.5
0.5
40
12
90
Trial 2
21.7
0.23
0.5
0.45
45
13.5
90
Trial 3
21.7
0.23
0.45
0.4
45
13.5
90
Trial 4
21.7
0.23
0.35
0.3
45
13.5
90Table 3
Trial input parameters for calibration of material B.
E (GPa)
v
tan
φ
1
tan
φ
2
C (MPa)
T (MPa)
σ
n
critical (MPa)
Trial 1
11.9
0.21
0.6
0.6
50
15
90
Trial 2
11.9
0.21
0.4
0.3
20
6
45
Trial 3
11.9
0.21
0.35
0.25
20
6
40
Trial 4
11.9
0.21
0.25
0.1
20
6
60
Trial 5
11.9
0.21
0.25
0.1
20
6
80
Trial 6
11.9
0.21
0.3
0.1
18
5.4
60Figure 8
Calibration process of interlayered materials.Figure 9
Calibration process of interlayered materials.Table 4
Calibration results of microparameters.
E (GPa)
v
tan
σ
n
critical
tan
σ
n
critical
C (MPa)
T (MPa)
σ
n
critical (MPa)
Material A
21.7
0.23
0.35
0.3
45
13.5
90
Material B
11.9
0.21
0.25
0.1
20
6
80


### 4.2. Microparameters of Layer Interfaces

To determine microparameters of layer interfaces, sensitivity analysis on macroresponse was conducted. Seven specimens with different bedding plane orientations are generated as shown in Figure10. The size of these specimens is 100 cm × 50 cm. The layer thickness is 10 cm. The following expression is used to define the microparameters of layer interfaces:(17)Parai=γ⋅ParaA+1−γParaB,where Para indicates all the microparameters, that is, E, v, tanφ1, tanφ2, C, T, and σncr; the subscript index i, A, and B, respectively, represent layer interface, material A and material B; and γ is a ratio coefficient. Letting γ be 0, 0.25, 0.5, 0.75, and 1 and using microparameters listed in Table 2, compressive strengths of the seven specimens are computed with the confining pressure of 0 MPa and 14 MPa. Results of this numerical experiment are shown in Figures 11 and 12. From Figure 11, it can be observed that, with the confining pressure of 0 MPa, increase of interface parameters can only lead to strength increase at α = 15°,30°, and 45°, when γ<0.5. When γ≥0.5, varying interface parameters has no effects on compressive strength at any orientation. From Figure 12, it can be observed that, with the confining pressure of 14 MPa, compressive strength generally increases as the increase of γ at all orientations apart from α = 0° and α = 90°. Finally, γ=0.5 is chosen as a coefficient for determination of layer interface parameters, since with which the computed compressive strengths match much better with experimental results than any other value of γ considering both the confining pressure of 0 MPa and 14 MPa.Figure 10
Numerical specimens with structural planes at different orientation: (a)α = 0°; (b) α = 15°; (c) α = 30°; (d) α = 45°; (e) α = 60°; (f) α = 75°; (g) α = 90°.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 11
Uniaxial compressive strength with different microparameters on material interface.Figure 12
Compressive strength under the confining pressure of 14 MPa with different microparameters on material interface.


### 4.3. Determination of Layer Thickness

In order to evaluate effects of layer thickness on strength anisotropy and failure modes, three groups of numerical tests are carried out. Each group has 7 specimens with same layer thickness, but different bedding plane orientation, respectively, 0°, 15°, 30°, 45°, 60°, 75°, and 90°. Layer thicknesses of different groups are not the same, respectively, 5 cm, 10 cm, and 16.67 cm. Here,d is used to represent layer thickness. With microparameters listed in Table 2 and γ=0.5, failure process of these 21 specimens are computed under uniaxial compression. Uniaxial compressive strengths of these specimens are illustrated in Figure 13. From this figure, it can be observed that, at α = 0°, 30°, 45°, and 60°, strengths of the group of d = 5 cm almost coincide with the group of d = 10 cm. At other orientations, strengths of d = 10 cm are a little lower than those of d = 5 cm. At α = 0° and 45°, strengths of the group of d = 16.67 cm are the same as those of the other two groups. But at other orientations, strengths of this group are a little lower than those of the other two groups. However, generally speaking, difference of strength among the three groups is not quite great. Illustrated in Figure 14 are failure modes when α = 30° from the three different groups. Shown in Figure 15 are microcrack distributions after failure for α = 0°, α = 45°, and α = 90° under triaxial compression with different confining pressures. Tensile and shear crack are, respectively, represented by a blue and red line. We can see that failure modes for different layer thicknesses are almost the same: tensile splitting across the weak layer. Failure modes at other orientations also share this point. Consequently, it can be said that layer thickness has limited influence on both strength anisotropy and failure modes in some extent. But from a theoretical viewpoint, to make simulation more representative, there should be enough gaps between block size, layer thickness, and specimen size. But for numerical efficiency, block size should not be too large. Finally, to make a compromise between representativeness and efficiency, d = 10 cm is chosen as the layer thickness.Figure 13
Uniaxial strength of numerical specimens at different orientations with different layer thickness.Figure 14
Failure mode of specimens atα = 30° with different layer thickness under uniaxial compression: (a) d = 5 cm; (b) d = 16.67 cm; (c) d = 16.67 cm.
(a)
(b)
(c)Figure 15
Macroelastic modulus obtained by numerical modeling.


### 4.4. Anisotropy in Elasticity

The macroelastic modulus obtained for specimens with different layer orientations is illustrated in Figure15. It can be noted that the peak value of elastic modulus takes place when α = 0°, and the lowest value occurs when α = 45°. Elastic property shows an obvious anisotropic behavior.


### 4.5. Failure Process Modeling

Failure processes of the 7 specimens illustrated in Figure10 are simulated under compression with different confining pressures, respectively, 0 MPa, 6 MPa, 14 MPa, and 35 MPa in accordance with Tien et al.’s experiments [6]. In the simulation, microparameters listed in Table 2 are employed, γ=0.5, and axial loading is controlled by displacement.Shown in Figure16 are stress-strain curves of different specimens under uniaxial compression. It can be observed that strength and deformation behaviors both show obvious anisotropy. Illustrated in Figure 17 is the comparison of strength between numerical tests and experimental results at different orientations and different confining pressures. As can be seen, strength anisotropy is well produced with the numerical model, and numerical results generally match well with experimental results.Figure 16
Stress-strain curves for each specimen under uniaxial compression.Figure 17
Comparison of strength between simulation and experiments under different confining pressures (continuous lines represent simulation results, and scattered points represent experiment result).Illustrated in Figure18 are macrofailure modes of the seven specimens under compression with different confining pressures. From these figures, failure modes under uniaxial compression can be categorized into three groups: (1) tensile splitting in the weak material, at α = 0° and 15°; (2) tensile splitting and sliding along weak planes, at α = 30°, 45°, and 60°; and (3) tensile splitting across weak planes, at α = 75° and 90°. These failure modes are basically in line with experimental observations shown in Figure 19.Figure 18
Failure modes of specimens with different orientations under compression with different confining pressures.Figure 19
Failure modes for specimens with different layer orientations under uniaxial compression observed in laboratory (from Tien et al. [6]): (a) α = 45°; (b) α = 30°; (c) α = 15°; (d) α = 90°; (e) α = 75°; (f) α = 60°.
(a)
(b)
(c)
(d)
(e)
(f)With confining pressures, tensile splitting effects are constrained. As confining pressure increases, failure modes tend to be dominated by shearing and can be roughly categorized into two groups: (1) shearing along the layer interface, atα = 30° and 45° with all confining pressures and at α = 15° with the confining pressure of 6 MPa and 14 MPa and (2) shearing across layer interfaces, at α = 0°, 60°, 75°, and 90° with all confining pressures and at α = 15° with the confining pressure of 35 MPa. Specifically, at α = 15°, failure modes are in the first group with relatively low confining pressure, but as confining pressure increases to 35 MPa, failure mode becomes shearing across weak planes, which indicates that effects of weak planes on the failure mode are weakened by confining pressure. These failure modes are also in line with experimental observations.From Figures17 and 18, it can be observed that the failure mode is directly linked to strength. At the same confining pressure, compressive strength for specimens with a failure mode of shearing sliding along weak planes is much lower than that of shearing across weak planes. The specimen of α = 30° has the lowest compressive strength at all confining pressures.


### 4.6. Mesh Dependency Analysis

In order to investigate the mesh dependency of mechanical strength obtained by the proposed RBSM model, additional simulations are conducted on two other groups of specimens generated with different values ofbmin, say 0.2 mm and 0.3 mm, respectively. Comparisons of peak deviatoric stress between three different groups of specimens are illustrated in Figure 20 for different loading orientations and confining pressures. As one can see, the differences of strength between three groups of specimens with different mesh sizes are quite small, in particular for low confining pressures (0 MPa). For the confining pressure of 20 MPa, some larger scatters are observed for the loading orientations between 45° and 75° but remain in an acceptable range.Figure 20
Comparison of failure strength for different mesh sizes with different layer orientations and confining pressures.


## 4.1. Calibration of Input Parameters for Material A and Material B

Input parameters of the two distinct materials are calibrated independently. They are identified through a trial and error optimization algorithm against the strength envelops of conventional triaxial compression test. In each step of calibration, simulation results are compared with experimental data at all confining pressures. Calibration progresses until simulation results match well with experimental data on the whole. Shown in Tables2 and 3 are, respectively, trial parameters for calibration of material A and material B. The corresponding numerical strength envelops comparing with experimental data are shown in Figures 8 and 9. The calibration process is through adjusting input parameters to gradually approach the experimental data. The finally obtained input parameters are listed in Table 4.Table 2
Trial input parameters for calibration of material A.
E (GPa)
v
tan
φ
1
tan
φ
2
C (MPa)
T (MPa)
σ
n
critical (MPa)
Trial 1
21.7
0.23
0.5
0.5
40
12
90
Trial 2
21.7
0.23
0.5
0.45
45
13.5
90
Trial 3
21.7
0.23
0.45
0.4
45
13.5
90
Trial 4
21.7
0.23
0.35
0.3
45
13.5
90Table 3
Trial input parameters for calibration of material B.
E (GPa)
v
tan
φ
1
tan
φ
2
C (MPa)
T (MPa)
σ
n
critical (MPa)
Trial 1
11.9
0.21
0.6
0.6
50
15
90
Trial 2
11.9
0.21
0.4
0.3
20
6
45
Trial 3
11.9
0.21
0.35
0.25
20
6
40
Trial 4
11.9
0.21
0.25
0.1
20
6
60
Trial 5
11.9
0.21
0.25
0.1
20
6
80
Trial 6
11.9
0.21
0.3
0.1
18
5.4
60Figure 8
Calibration process of interlayered materials.Figure 9
Calibration process of interlayered materials.Table 4
Calibration results of microparameters.
E (GPa)
v
tan
σ
n
critical
tan
σ
n
critical
C (MPa)
T (MPa)
σ
n
critical (MPa)
Material A
21.7
0.23
0.35
0.3
45
13.5
90
Material B
11.9
0.21
0.25
0.1
20
6
80


## 4.2. Microparameters of Layer Interfaces

To determine microparameters of layer interfaces, sensitivity analysis on macroresponse was conducted. Seven specimens with different bedding plane orientations are generated as shown in Figure10. The size of these specimens is 100 cm × 50 cm. The layer thickness is 10 cm. The following expression is used to define the microparameters of layer interfaces:(17)Parai=γ⋅ParaA+1−γParaB,where Para indicates all the microparameters, that is, E, v, tanφ1, tanφ2, C, T, and σncr; the subscript index i, A, and B, respectively, represent layer interface, material A and material B; and γ is a ratio coefficient. Letting γ be 0, 0.25, 0.5, 0.75, and 1 and using microparameters listed in Table 2, compressive strengths of the seven specimens are computed with the confining pressure of 0 MPa and 14 MPa. Results of this numerical experiment are shown in Figures 11 and 12. From Figure 11, it can be observed that, with the confining pressure of 0 MPa, increase of interface parameters can only lead to strength increase at α = 15°,30°, and 45°, when γ<0.5. When γ≥0.5, varying interface parameters has no effects on compressive strength at any orientation. From Figure 12, it can be observed that, with the confining pressure of 14 MPa, compressive strength generally increases as the increase of γ at all orientations apart from α = 0° and α = 90°. Finally, γ=0.5 is chosen as a coefficient for determination of layer interface parameters, since with which the computed compressive strengths match much better with experimental results than any other value of γ considering both the confining pressure of 0 MPa and 14 MPa.Figure 10
Numerical specimens with structural planes at different orientation: (a)α = 0°; (b) α = 15°; (c) α = 30°; (d) α = 45°; (e) α = 60°; (f) α = 75°; (g) α = 90°.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 11
Uniaxial compressive strength with different microparameters on material interface.Figure 12
Compressive strength under the confining pressure of 14 MPa with different microparameters on material interface.


## 4.3. Determination of Layer Thickness

In order to evaluate effects of layer thickness on strength anisotropy and failure modes, three groups of numerical tests are carried out. Each group has 7 specimens with same layer thickness, but different bedding plane orientation, respectively, 0°, 15°, 30°, 45°, 60°, 75°, and 90°. Layer thicknesses of different groups are not the same, respectively, 5 cm, 10 cm, and 16.67 cm. Here,d is used to represent layer thickness. With microparameters listed in Table 2 and γ=0.5, failure process of these 21 specimens are computed under uniaxial compression. Uniaxial compressive strengths of these specimens are illustrated in Figure 13. From this figure, it can be observed that, at α = 0°, 30°, 45°, and 60°, strengths of the group of d = 5 cm almost coincide with the group of d = 10 cm. At other orientations, strengths of d = 10 cm are a little lower than those of d = 5 cm. At α = 0° and 45°, strengths of the group of d = 16.67 cm are the same as those of the other two groups. But at other orientations, strengths of this group are a little lower than those of the other two groups. However, generally speaking, difference of strength among the three groups is not quite great. Illustrated in Figure 14 are failure modes when α = 30° from the three different groups. Shown in Figure 15 are microcrack distributions after failure for α = 0°, α = 45°, and α = 90° under triaxial compression with different confining pressures. Tensile and shear crack are, respectively, represented by a blue and red line. We can see that failure modes for different layer thicknesses are almost the same: tensile splitting across the weak layer. Failure modes at other orientations also share this point. Consequently, it can be said that layer thickness has limited influence on both strength anisotropy and failure modes in some extent. But from a theoretical viewpoint, to make simulation more representative, there should be enough gaps between block size, layer thickness, and specimen size. But for numerical efficiency, block size should not be too large. Finally, to make a compromise between representativeness and efficiency, d = 10 cm is chosen as the layer thickness.Figure 13
Uniaxial strength of numerical specimens at different orientations with different layer thickness.Figure 14
Failure mode of specimens atα = 30° with different layer thickness under uniaxial compression: (a) d = 5 cm; (b) d = 16.67 cm; (c) d = 16.67 cm.
(a)
(b)
(c)Figure 15
Macroelastic modulus obtained by numerical modeling.


## 4.4. Anisotropy in Elasticity

The macroelastic modulus obtained for specimens with different layer orientations is illustrated in Figure15. It can be noted that the peak value of elastic modulus takes place when α = 0°, and the lowest value occurs when α = 45°. Elastic property shows an obvious anisotropic behavior.


## 4.5. Failure Process Modeling

Failure processes of the 7 specimens illustrated in Figure10 are simulated under compression with different confining pressures, respectively, 0 MPa, 6 MPa, 14 MPa, and 35 MPa in accordance with Tien et al.’s experiments [6]. In the simulation, microparameters listed in Table 2 are employed, γ=0.5, and axial loading is controlled by displacement.Shown in Figure16 are stress-strain curves of different specimens under uniaxial compression. It can be observed that strength and deformation behaviors both show obvious anisotropy. Illustrated in Figure 17 is the comparison of strength between numerical tests and experimental results at different orientations and different confining pressures. As can be seen, strength anisotropy is well produced with the numerical model, and numerical results generally match well with experimental results.Figure 16
Stress-strain curves for each specimen under uniaxial compression.Figure 17
Comparison of strength between simulation and experiments under different confining pressures (continuous lines represent simulation results, and scattered points represent experiment result).Illustrated in Figure18 are macrofailure modes of the seven specimens under compression with different confining pressures. From these figures, failure modes under uniaxial compression can be categorized into three groups: (1) tensile splitting in the weak material, at α = 0° and 15°; (2) tensile splitting and sliding along weak planes, at α = 30°, 45°, and 60°; and (3) tensile splitting across weak planes, at α = 75° and 90°. These failure modes are basically in line with experimental observations shown in Figure 19.Figure 18
Failure modes of specimens with different orientations under compression with different confining pressures.Figure 19
Failure modes for specimens with different layer orientations under uniaxial compression observed in laboratory (from Tien et al. [6]): (a) α = 45°; (b) α = 30°; (c) α = 15°; (d) α = 90°; (e) α = 75°; (f) α = 60°.
(a)
(b)
(c)
(d)
(e)
(f)With confining pressures, tensile splitting effects are constrained. As confining pressure increases, failure modes tend to be dominated by shearing and can be roughly categorized into two groups: (1) shearing along the layer interface, atα = 30° and 45° with all confining pressures and at α = 15° with the confining pressure of 6 MPa and 14 MPa and (2) shearing across layer interfaces, at α = 0°, 60°, 75°, and 90° with all confining pressures and at α = 15° with the confining pressure of 35 MPa. Specifically, at α = 15°, failure modes are in the first group with relatively low confining pressure, but as confining pressure increases to 35 MPa, failure mode becomes shearing across weak planes, which indicates that effects of weak planes on the failure mode are weakened by confining pressure. These failure modes are also in line with experimental observations.From Figures17 and 18, it can be observed that the failure mode is directly linked to strength. At the same confining pressure, compressive strength for specimens with a failure mode of shearing sliding along weak planes is much lower than that of shearing across weak planes. The specimen of α = 30° has the lowest compressive strength at all confining pressures.


## 4.6. Mesh Dependency Analysis

In order to investigate the mesh dependency of mechanical strength obtained by the proposed RBSM model, additional simulations are conducted on two other groups of specimens generated with different values ofbmin, say 0.2 mm and 0.3 mm, respectively. Comparisons of peak deviatoric stress between three different groups of specimens are illustrated in Figure 20 for different loading orientations and confining pressures. As one can see, the differences of strength between three groups of specimens with different mesh sizes are quite small, in particular for low confining pressures (0 MPa). For the confining pressure of 20 MPa, some larger scatters are observed for the loading orientations between 45° and 75° but remain in an acceptable range.Figure 20
Comparison of failure strength for different mesh sizes with different layer orientations and confining pressures.


## 5. Conclusions

The modified rigid block spring method (RBSM) is employed for simulation of failure process of an artificial interlayered rock. Rock mass is replaced by the Voronoi based block system. Microcracks are explicitly modeled and propagate along block interfaces. A combined criterion is used to take into consideration crack events: the Mohr–Coulomb criterion is used to detect shearing failure events and the tensile strength criterion for detection of tensile failure events. Comparison with Tien et al.’s experimental results [6] leads to the following conclusions:(1)
The modified RBSM model has the capacity to describe anisotropic behaviors of interlayered rocks in both deformation and compressive strength qualitatively and quantitatively.(2)
Using the modified RBSM model, failure modes are successfully captured under compression with different confining pressures. As confining pressure increases, effects of structural planes on failure modes become weakened.The present work is based on two-dimensional conditions. In the near future, this model will be extended to three-dimensional conditions.



---

*Source: 9615457-2018-08-01.xml*