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.

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.

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° [16]. 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.

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. [1720], 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 [2731]. 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, , can be expressed by the displacements of centroids of the two blocks aswherein which and are, respectively, the relative normal and tangential displacements between point P1 and point P2.where is the unit normal vector of the interface.where , , , and .where , , and are, respectively, the translational displacement in the x direction, y direction, and the rotational angle of the centroid of block 1, while , , and are displacement components for block 2.

Stress induced by relative displacements between point P1 and point P2 can be expressed aswherein which and are normal stress and tangential stress.where is the stiffness of the normal spring and is the stiffness of the tangential spring.

According to the Virtual Work Theorem, for the block system, the relationship is as follows:where , , and , respectively, stand for interfaces between blocks, force boundaries, and domains of blocks; , , , and 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:

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 Figure 4 and is composed of two parts, the compressive and the tensile components. In this chapter, tension is defined as positive.

In compression, σn is given aswhere is the relative displacement between interacting blocks and 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

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 aswhere is the updated shear stress, is the tangential stiffness, proportional to , and 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 , according to Yao et al. [27]:where and are intermediate variables and and 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 Figure 5. 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

Shear rupture takes place when , and then the interaction becomes purely frictional, with a maximum shear force defined by

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 distance . 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 vertex V 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 . Then, draw auxiliary lines which have an angle of with respect to each segment. After that, draw a circle with a radius of 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.

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 than . Shown in Figure 7(c) is the step of segment definition.

3.1.3. Insertion of Points to the Whole Domain

With the constraint of , 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.

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 Tables 2 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.

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 Figure 10. 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:where indicates all the microparameters, that is, , , , , , , and ; the subscript index , , and , 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 . When , 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, 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.

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 , 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.

4.4. Anisotropy in Elasticity

The macroelastic modulus obtained for specimens with different layer orientations is illustrated in Figure 15. 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 Figure 10 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, , and axial loading is controlled by displacement.

Shown in Figure 16 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.

Illustrated in Figure 18 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.

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 Figures 17 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 of bmin, 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.

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.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The work reported in this paper has received financial support from the National Natural Science Foundation of China (nos. U1765207, 51509104, and 41762020), Key Laboratory of New Technology for Construction of Cities in Mountain Area, Administration of Education, China (0902071812102/016), and State Key Laboratory of Geomechanics and Geotechinical Engineering, China (Z017021). These supports are gratefully acknowledged.