Failure Process Simulation of Interlayered Rocks under Compression

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. )e 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. )e modified rigid block spring method is used to simulate the initiation and propagation of microcracks. )e 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.


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 ° [1][2][3][4][5][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.
ese 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 di erent con ning pressures with PFC3D.Wang et al. [14] analyzed the Brazilian tensile splitting tests of strati ed biotite granulite rocks with PFC2D.e computed stress-strain curves were compared well with experimental results.Shen et al. [15] simulated a layered rock under triaxial compression using a nite element method.Lisjak et al. [16] simulated the anisotropic behaviors of Opalinus argillite by the combined FEM/DEM model.Zhao et al. [17][18][19][20], Wang et al. [21], and Wei et al. [22] have also made great contributions to this eld.
e 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.e blocks are bonded by contact interfaces.ree springs are de ned on each interface to describe the relative normal and tangential displacements and rotations between two neighboring blocks.e initial RBSM method was so far successfully applied to modeling macroscopic behaviors of cement-based materials [24,25].
is method has been successively improved by Qian and Zhang [26].In the modi ed RBSM model, a continuous distribution of stress and relative displacements is de ned on each interface through suitable interpolation functions.
is allows to modeling the progressive failure of interfaces.e modi ed method was recently successfully applied to modeling the damage and failure of brittle rocks [27][28][29][30][31]. e advantages of this method have been discussed.is 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 arti cial damping coe cient.
However, the modi ed RBSM model has been so far applied to initially isotropic materials.
e 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.
e proposed anisotropic RBSM model will be applied to modeling the mechanical behavior of typical arti cial interlayered rocks, which are carried out by Tien et al. [6].In particular, in uences of the initial structural anisotropy on the induced damage and failure process will be discussed.

Principles of the Modified Rigid
Block Method e governing equations for the modi ed RBSM are derived from the Virtual Work eorem, which is simply given below.Suppose that two arbitrary blocks, block 1 and block 2, are next to each other.e centroid of the two blocks are (x 1 , y 1 ) and (x 2 , y 2 ), respectively.Point P 1 on block 1 and point P 2 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 P 1 and point P 2 , Δu { }, can be expressed by the displacements of centroids of the two blocks U { } 12 as where in which Δu n and Δu s are, respectively, the relative normal and tangential displacements between point P 1 and point P 2 .
[B] l m where l, m { } is the unit normal vector of the interface.
where U 1x , U 1y , and U 1θ are, respectively, the translational displacement in the x direction, y direction, and the rotational angle of the centroid of block 1, while U 2x , U 2y , and U 2θ are displacement components for block 2.  Advances in Civil Engineering Stress induced by relative displacements between point P 1 and point P 2 can be expressed as where in which σ n and σ s are normal stress and tangential stress. [D] where k n is the sti ness of the normal spring and k s is the sti ness of the tangential spring.
According to the Virtual Work eorem, for the block system, the relationship is as follows: where l e 0 , l e σ , and s e , 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: A linear elastic relation is employed to describe the stress and displacement relation between two neighboring blocks.e 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 de ned as positive.
In compression, σ n is given as where Δu n is the relative displacement between interacting blocks and k n is the normal sti ness.
In tension, the normal stress is still computed with the same sti ness in compression.e 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 de ned as

Advances in Civil Engineering
where σ s updated is the updated shear stress, k s is the tangential sti ness, proportional to k n , and Δu s 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, k n and k s can be determined by macroelastic parameters, that is, elastic modulus E and Poisson's ratio v, according to Yao et al. [27]: where r and E 0 are intermediate variables and h 1 and h 2 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 modi ed Mohr-Coulomb model is adopted, as shown in Figure 5. e 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 .e critical normal stress σ n critical is de ned to limit the frictional strengthening e ects.Before rupture, the maximum shear stress is computed by Shear rupture takes place when σ s > σ s max , and then the interaction becomes purely frictional, with a maximum shear force de ned by In the modi ed rigid block spring method, structural planes and virtual cracks are treated with the same failure criterion but di erent failure parameters.

Mesh Generation
Two main steps are taken to generate Voronoi diagrambased mesh: point insertion and tessellation.Each step is described in detail as follows.

Point Insertion.
e 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 l min .To explicitly model bedding planes, a multistep insertion procedure is adopted.Here, we de ne a segment of boundaries or bedding planes as a segment and an intersection between boundaries or bedding planes as a vertex.

Insertion of Points to De ne Vertexes.
Suppose there is a vertex V connected by four segments, as shown in Figure 6.To de ne this vertex, rstly, nd out the minimum angle between two neighboring segments, the value of which is assumed to be 2β.en, draw auxiliary lines which have an angle of β with respect to each segment.After that, draw a circle with a radius of 0.5l min around vertex V. Insert points at intersections between the circle and all auxiliary lines   Advances in Civil Engineering

Tessellation.
ere 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 )).

Numerical Example
Using arti cial 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 di erent con ning pressures.e arti cial rock is formed from strati ed layers made of material A and material B.
icknesses of layers for the two materials are the same.But mechanical properties of them are quite di erent.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.

Calibration of Input Parameters for Material A and Material B.
Input parameters of the two distinct materials are calibrated independently.ey are identi ed 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 con ning 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.
e corresponding numerical strength envelops comparing with experimental data are shown in Figures 8 and 9. e calibration process is through adjusting input parameters to gradually approach the experimental data.
e nally obtained input parameters are listed in Table 4.

Microparameters of Layer Interfaces.
To determine microparameters of layer interfaces, sensitivity analysis on macroresponse was conducted.Seven specimens with di erent bedding plane orientations are generated as shown in Figure 10.e size of these specimens is 100 cm × 50 cm.e layer thickness is 10 cm. e following   6 Advances in Civil Engineering expression is used to de ne the microparameters of layer interfaces: where Para indicates all the microparameters, that is, E, v, tan φ 1 , tan φ 2 , C, T, and σ n cr ; the subscript index i, A, and B, respectively, represent layer interface, material A and material B; and c is a ratio coe cient.Letting c 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 con ning 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 con ning pressure of 0 MPa, increase of interface parameters can only lead to strength increase at α 15 °,30 °, and 45 °, when c < 0.5.When c ≥ 0.5, varying interface parameters has no e ects on compressive strength at any orientation.From Figure 12, it can be observed that, with the con ning pressure of 14 MPa, compressive strength generally increases as the increase of c at all orientations apart from α 0 °and α 90 °.Finally, c 0.5 is chosen as a coe cient for determination of layer interface parameters, since with which the computed compressive strengths match much better with experimental  Advances in Civil Engineering results than any other value of c considering both the con ning pressure of 0 MPa and 14 MPa.But at other orientations, strengths of this group are a little lower than those of the other two groups.However, generally speaking, di erence of strength among the three groups is not quite great.Illustrated in Figure 14 are failure modes when α 30 °from the three di erent groups.Shown in Figure 15 are microcrack distributions after failure for α 0 °, α 45 °, and α 90 °under triaxial compression with di erent con ning pressures.Tensile and shear crack are, respectively, represented by a blue and red line.We can see that failure modes for di erent 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 in uence 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 e ciency, block size should not be too large.Finally, to make a compromise between representativeness and e ciency, d 10 cm is chosen as the layer thickness.

Anisotropy in Elasticity.
e macroelastic modulus obtained for specimens with di erent layer orientations is  Advances in Civil Engineering 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.16 are stress-strain curves of di erent 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 di erent orientations and di erent con ning 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 di erent conning pressures.From these gures, 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 °. ese failure modes are basically in line with experimental observations shown in Figure 19.
With con ning pressures, tensile splitting e ects are constrained.As con ning 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 con ning pressures and at α 15 °with the con ning pressure of 6 MPa and 14 MPa and (2) shearing across layer interfaces, at α 0 °, 60 °, 75 °, and 90 °with all con ning pressures and at α 15 °with the con ning pressure of 35 MPa.Speci cally, at α 15 °, failure modes are in the rst group with relatively low con ning pressure, but as con ning pressure increases to 35 MPa, failure mode becomes shearing across weak planes, which indicates that e ects of weak planes on the failure mode are weakened by con ning pressure.
ese 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 con ning 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.e specimen of α 30 °has the lowest compressive strength at all con ning 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 di erent values of b min , say 0.2 mm and 0.3 mm, respectively.Comparisons of peak deviatoric stress between three di erent groups of specimens are illustrated in Figure 20 for di erent loading orientations and con ning pressures.As one can see, the di erences of strength between three groups of specimens with di erent mesh sizes are quite small, in particular for low con ning  Advances in Civil Engineering pressures (0 MPa).For the con ning pressure of 20 MPa, some larger scatters are observed for the loading orientations between 45 °and 75 °but remain in an acceptable range.

Conclusions
e modi ed rigid block spring method (RBSM) is employed for simulation of failure process of an arti cial 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) e modi ed RBSM model has the capacity to describe anisotropic behaviors of interlayered rocks in both deformation and compressive strength qualitatively and quantitatively.(2) Using the modi ed RBSM model, failure modes are successfully captured under compression with different con ning pressures.As con ning pressure increases, e ects of structural planes on failure modes become weakened.
e present work is based on two-dimensional conditions.In the near future, this model will be extended to three-dimensional conditions.
Data Availability e data used to support the ndings of this study are available from the corresponding author upon request.

Figure 1 :
Figure 1: De nition of structural plane orientation.

Figure 3 :
Figure 3: Local deformation of an interface.

Figure 7 :
Figure 7: Illustration of mesh generation process.

4. 3 .
Determination of Layer ickness.In order to evaluate e ects 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 di erent 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 c 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 gure, 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.

Figure 13 :Figure 14 :Figure 15 :
Figure 13: Uniaxial strength of numerical specimens at di erent orientations with di erent layer thickness.

Table 2 :
Trial input parameters for calibration of material A.

Table 3 :
Trial input parameters for calibration of material B.
[6].Failure Process Modeling.Failure processes of the 7 specimens illustrated in Figure10are simulated under compression with di erent con ning 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 Table2are employed, c 0.5, and axial loading is controlled by displacement.Shown in Figure