Abstract

Bedding planes are the primary control on the anisotropy of mechanical characteristics and fracture patterns in rock. To analyze the influence of the geometrical properties of bedding planes on the direct shear strength characteristics and fracture patterns of transversely isotropic rocks, numerical models were established using an improved modeling method using Particle Flow Code. The results of the numerical model were in good agreement with those of the physical experiments of an artificial rock mass containing a single bedding plane. The results show that the shear fractures with a range of bedding plane geometries can be divided into two patterns. When the inclination angles of the bedding planes were larger or smaller, a thoroughgoing fracture plane was formed near the preexisting shear fracture plane. On the other hand, the intact rock was broken into many parallel sheets.

1. Introduction

Transversely isotropic rocks have undergone a variety of geological and tectonic evolution processes during the process of formation. They often contain weak planes, such as bedding, foliation, or schistosity, which contribute significantly to anisotropic deformation and strength [14]. Furthermore, this type of rock is often encountered in numerous geological applications and various engineering projects, for example, exploitation of shale gas [57] and stability analyses of large underground caverns and other rock surrounding mining or traffic tunnels [8, 9]. In addition, both the mechanical properties and damage mechanisms of transversely isotropic rocks show strong anisotropy, which poses a huge challenge to practical engineering. Therefore, it is of great practical importance to study the mechanical properties of transversely isotropic rock masses.

With the development of experimental techniques and numerical simulation methods, the study of transversely isotropic rocks has been renewed [1013]. Many physical experimental studies have been performed on different rock-like and natural materials. Tien and Tsao prepared artificial transversely isotropic rocks and carried out uniaxial and triaxial compression tests to study the compressional deformation anisotropy of layered structures [14, 15]. Cho et al. investigated the dependency of the determined elastic constants by using a number of strain measurements and mechanical properties from uniaxial compression tests, triaxial compression tests, and Brazilian splitting tests with three anisotropic rocks that showed clear evidence of transverse anisotropy and that were from Korea [3]. Liao et al. conducted a series of direct tensile tests on argillite to study the tensile properties of transversely isotropic rocks [16]. Vervoort and Tavallali [1719] studied the induced fracture patterns during Brazilian tests of various transversely isotropic rock materials. However, all the studies mentioned above focused on the anisotropy of the strength and elastic deformation parameters and the fracture characteristics under compression or tension. There are few little studies on the mechanical properties in the shear failure process of transversely isotropic rocks. In addition, studies on the shear strength of a rock mass are often limited to the shear mechanical properties of the weakest structures [2025].

Shuai et al. [26, 27] investigated the shear strength of shale at four different inclinations (0°, 30°, 60°, and 90°) and explored the relationship between the peak shear strength and the inclination. However, due to the limitations of the experimental specimens, no research has been conducted on the influence of a wider range of inclination or other geometrical properties of the bedding planes on the shear strength characteristics and fracture patterns of transversely isotropic rocks. Thus, based on previous studies [2833], improved numerical models were established by using the (Particle Flow Code) software to explore the shear strength and fracture patterns of transversely isotropic rocks under various normal pressures and geometrical parameters of the bedding planes. The numerical models were verified by the direct shear test of the artificial rocks that contained a single bedding plane. can monitor the characteristic of microcracks in the samples (position, direction, type, etc.) during the loading process to analyze the fracture patterns at the microscale. In addition, compared with experimental tests, numerical simulation has the advantage of lower costs, shorter experimental runtimes, and sufficient data collection. Therefore, a numerical method such as a Particle Flow Code can better analyze the fracture pattern and mechanism of transversely isotropic rocks during direct shear deformation.

2. Improved Numerical Model

The bond removal method in the has been widely used to model the bedding planes in a rock mass by PFC [20, 21, 3437], which sets the bond strength to zero within bedding planes (usually by command JSET). The drawbacks of this method can be clearly seen in the simulation of the direct shear tests within a planar joint. Due to the random arrangement of particles, the interlocking effect of particles on either side of the bedding plane is significant, which results in the inherent roughness of the bedding planes at the microscale, as shown in Figure 1(a). In addition, the particles within the planar joint have a “climbing” effect during the direct shear test due to the rigidity of the particles, as shown in Figure 2(b). This microscale movement of particles leads to a local stress concentration and unreasonable normal expansion. The shear stress primarily increases due to the multiple stress concentrations, and the value of the shear stress is generally much higher than that in nature. After the failure of the structure with the “climbing” particles or after the “climbing” has been completed, the shear capacity of the joint suddenly decreases, which releases the stored energy and causes energy dissipation. Therefore, in the resulting shear stress-shear displacement curve, this process is reflected by a much higher peak shear strength followed by a rapid decrease [28, 35], as shown in Figure 3.

To overcome the shortcomings of this approach, the smooth joint model (SJM) was introduced into the [38]. In numerical models, the SJM can be applied at the contacts between particles that have centers positioned on opposite sides of the joint interface, similar to the JSET command in . However, unlike the bond removal method the contacting particles can overlap and pass through each other, which can eliminate the “climbing” effect, as shown in Figure 2(c).

However, Bahaaddini et al. [29] and Mehranpour and Kulatilake [30] noted certain defects when using the SJM. For example, the stress concentrations also occurred during the shearing process in the SJM if the shear displacement exceeded a certain value. When the SJM was applied to the contacts, the original contact models were replaced. Due to the nonuniformity of the particle radii and the randomness of the particle positions, as shown in Figure 1(b), the positions of the particles changed during the shearing process, resulting in the formation of new contacts. In addition, some new contacts were not applied with the SJM, since both particles might be on the same side of the joint interface. Therefore, a slight “climbing” effect would also appear, which led to an unreasonable increase and drop in the shear stress curve, as shown in Figure 3. The details of the SJM defects were analyzed in [29].

In this paper, an improved smooth joint model was adopted based on the previous research [28, 29]. In the process of generating the specimen, temporary walls of a specific stiffness were added where the bedding planes would be later added to separate the particles in the upper and lower blocks, resulting in a substantial reduction in the connectivity between particles at the bedding planes. Then, a specific isotropic stress was applied and the particles with less than three contacts (floaters) were eliminated. Next, the temporary walls were deleted. At this stage, the particles on either side of the bedding planes were completely separated, which caused the mechanical parameters of the bedding planes to be weak. Therefore, in order to force the particles at the bedding surfaces to be in contact, the floater-elimination procedure was repeated.

The floater-elimination procedure guaranteed generating a dense particle pack, but not just by deleting the floaters. Firstly, all particles except the floaters were fixed and given velocities of zero. Floaters were then sufficiently expanded (30%) to ensure contact with all their immediate neighbors. After being cycled to local equilibrium, the floaters were contracted by a calculated amount to decrease their maximum contact force within a given tolerance. The contraction step was applied repeatedly until all floaters were inactive [39].

The final numerical model is shown in Figure 1(c), in which the particles within the bedding planes were arranged more neatly and can eliminate the “climbing” effect more effectively than using simple SJM, as shown in Figure 1(b). Furthermore, the particles on the both sides of the bedding planes are ensured to be in good contact with each other. In addition, the shear strength envelopes of a planar joint were also the most reasonable with the final numerical model, as shown in Figure 3.

3. Physical Direct Shear Tests of Rock Masses with a Single Bedding Plane

In this paper, the physical direct shear tests of rock masses with single joints were conducted and compared with the numerical simulations to verify the numerical method of simulating the direct shear tests of jointed rock masses.

3.1. Preparation of the Samples

According to the previous physical experimental research, the gypsum is a good testing material because of its simple structure, low price, and short curing time [40]. To obtain more reasonable strength characteristics of the rock mass, the high-strength gypsum was used to produce the test specimens. In addition, the water : gypsum ratio was set as 0.35 : 1.0, which created a high-strength material while ensuring a certain fluidity to facilitate pouring during specimen preparation. To delay the curing time of the gypsum to allow enough time for stirring and pouring, a retarder was added at a weight ratio of 0.1%. The effect of the additional retarder on the overall mechanical properties of the specimen was negligible. The size of the specimens used in the direct shear tests was 100 m00 m00 mm. The prepared cube specimens were cut along the horizontal plane at inclination angles of 0°, 30°, 60°, and 90°, as shown in Figure 4. In addition, the separate specimens were then plastered together; the plaster included a custom cement material mixed with putty powder and construction glue in a 2 : 1 weight ratio. After the glue dried, the specimens containing single bedding planes with inclination angles of 0°, 30°, 60°, and 90° were ready for testing, as shown in Figure 5.

3.2. Direct Shear Tests

The direct shear tests were carried out using an RMT-301 servohydraulic mechanical testing system for rock and concrete (Figure 6) at Wuhan University, China.

First, in order to obtain the shear mechanical parameters of the intact rock and the bedding plane, the intact specimens (without bedding planes) and the specimens with a single planar bedding plane were prepared for direct shear tests under various normal stress conditions. The test results are shown in Figure 7, which were used to calibrate the microscale parameters of the numerical model. Then, the direct shear tests of the specimens with a single bedding plane with various inclination angles were carried out, and the results are discussed in the next section.

4. Establishment and Verification of the Numerical Models

4.1. Numerical Simulation of the Direct Shear Tests

In this paper, the numerical models were established in . The specimen size was 100 m00 mm. The radius of the circular particles was given by a Gaussian distribution: the particle radius ratio was 1.66 and the minimum radius was 0.35 mm. A single specimen contained 13206 total particles. Normal stress was applied vertically to the top of the specimen, and the normal stress was kept constant during the shearing process by using a servo-control mechanism [39]. The loading rate of shear box was set as 0.1 m/s. The calculation is based on a complex time-stepping algorithm; during the numerical shearing process, the time-step was approximately 5.8 × 10−8 s/step. Thus, the shear velocity in the numerical model was equal to 5.8 × 10−9 m/step, which was low enough to ensure that the specimen remained in a quasistatic equilibrium. In addition, the processing time for a shear displacement of 3 mm was approximately around 6 h using a 64-bit PC, with a speed of 3.6 GHz and 8.0 GB RAM.

The contacts in the intact rock were determined with the parallel bond model, while the contacts within the bedding planes were determined with the smooth joint model. The numerical model and shear box are shown in Figure 8.

If the upper and lower shear boxes were set too close, the particles at the presupposed shear plane would be pushed during the direct shear process when the inclined angle θ = 0°. Then the “climbing” effect might occur, resulting in the irregular shear stress jump. Thus, comparing to the maximum particle = 0.581 mm, the 2 mm interval between the shear boxes was set to avoid the above phenomenon during the shear process as shown in Figure 8.

4.2. Calibration of the Microscale Parameters

As one of the most rapidly developing areas of computational mechanics, the discrete element method (DEM) for particulate and blocky systems has undergone a significant development in recent decades [41]. As a discrete element approach, the major advantage of the bonded-particle model (BPM) is the replacement of the complex empirical constitutive behavior with simple particle contact logic [42]. Since the microscale parameters in the DEM do not directly correspond to the macroscale parameters in real materials, choosing reasonable microscale parameters is a challenge [43].

4.2.1. Calibration of the Bonded-Particle Model

The strength microscale parameters were calibrated by using direct shear tests with intact specimens identifying the shear strength characteristics (cohesive strength and friction angle ). In addition, the microscale deformation parameters were calibrated under uniaxial compression tests to match the elastic modulus and Poisson’s ratio . The uniaxial compression tests were carried out in BPM specimens with a width of 50 mm and a height of 100 mm (similar to the physical experiments). The calibrated microscale parameters are presented in Table 1.

A comparison between the results of the numerical and physical tests is presented in Table 2 and Figure 9, which show that the results of the numerical model were in good agreement with the physical tests.

4.2.2. Calibration of the Smooth Joint Model

The sensitivity analysis of each of the microscale parameters of the SJM has been conducted in order to calibrate the SJM, which can obtain the quantitative relationships between the macro- and microscale parameters. Combined with the previous studies [28, 29], the calibration of the SJM was conducted in three steps. Firstly, a normal uniaxial compression test was carried out to calibrate the microscale normal stiffness . Then, the microscale shear stiffness was calibrated with the macroscale shear stiffness by a direct shear test. Finally, the rest of the microscale parameters were calibrated with the shear peak strength .

The final selected SJM parameters and the comparison between the numerical and experimental results are presented in Tables 3 and 4, respectively.

4.3. Verification of the Numerical Model

The direct shear tests of the specimens containing a single planar bedding plane at different inclination angles (0°, 30°, 60°, and 90°) under 1.0 MPa of normal stress were conducted to verify the improved numerical model and the selection of microscale parameters. The numerical models are shown in Figure 10.

The shear stress-shear displacement curves are presented in Figure 11, in which the solid lines and points represent the numerical and physical results, respectively. The results of the comparison between the numerical and physical tests were in good agreement. The effectiveness of the numerical modeling method to simulate the direct shear tests of jointed rock mass was proven.

5. Numerical Simulation of the Shear Behavior of the Transversely Isotropic Rocks

Making full use of the advantages and convenience of the numerical methods, the direct shear tests of the specimens with various inclination angles and various bedding planes spacings were carried out by using the numerical models established above. The influences of the geometrical characteristics of the bedding planes on the direct shear mechanical properties and fracture modes of the transversely isotropic rock masses were then studied.

To simulate the low shear strength of the gypsum specimens the strength of the microscale parameters in the parallel bonds was relatively small. In this section, which emphasizes the difference between the mechanical characteristics of the bedding planes and the intact rock, the strength of the microscale parameters of the parallel bonds was increased ( = = 27.0 MPa). Additionally, the microscale stiffness of the SJM was relatively small in this new numerical model, leading to a large deformation and a longer computational runtime during shear failure. Thus, the microscale stiffness of the SJM was increased accordingly in the following numerical experiments ( = = 35 GPa/m).

With the above-mentioned modification, the direct shear test of the specimens with three different bedding planes spacings (spacing was 100 mm, 20 mm, and 10 mm) and seven different inclination angles (θ = 0°, 15°, 30°, 45°, 60°, 75°, and 90°) under five different normal pressures (1.0 MPa, 1.5 MPa, 2.0 MPa, 3.0 MPa, and 5.0 MPa) was carried out.

5.1. Influence on the Equivalent Macroscale Shear Mechanical Parameters

The equivalent macroscale shear mechanical parameters (equivalent internal friction angle and equivalent cohesion ) can be obtained by the linear fit with the Mohr-Coulomb criterion, presented in Figures 12 and 13.

5.1.1. Influence on Equivalent Internal Friction Angle

As seen in Figure 12, the internal friction angles were equivalent when the inclination angle is at θ = 0° and 90° under different bedding plane spacings . This is because the fitting results only reflect the shear mechanical parameters of the bedding planes when θ = 0°. In addition, when θ = 90°, the shear failure of the specimen was almost entirely in the intact rock, and the frictional effect of the specimen during the shear process was mostly borne by the intact rock. Thus, the bedding plane spacing had little effect on .

At the same time, the maximum occurred at about θ = 30° for all three different spacings . In addition, when θ = 30° increased to 90°, as the spacing decreased, the trends of and inclination angle gradually change from linearly decreasing to a parabolic curve that first decreased and then increased. This is because when the number of bedding planes was increased, the specimens at θ = 45° to 60° became more prone to fracture by slipping during the shearing process along the bedding planes. At this time, the shear strength of a specimen was primarily dominated by the bedding planes, and the equivalent internal friction angle was approximately the internal friction angle of the bedding planes, which was at θ = 0°.

Therefore, the influence of the bedding plane spacing on of under various bedding plane inclination angles can be divided into two stages, as described below. One stage was when inclination angle was θ = 0° to 30°. The influence of the spacing on was very small, and increased with the increasing inclination angles θ. The other stage was when θ = 30° to 90°. It can be seen that when bedding plane spacing is smaller, the weakening of can be observed more clearly.

5.1.2. Influence on the Equivalent Cohesion

Figure 13 shows that, under the same bedding plane spacing, the maximum equivalent cohesion was reached at approximately θ = 15°, while the minimum was reached at approximately θ = 45°. This is because when the inclination angle θ = 15°, the number of the bedding planes that intersected the preexisting shear plane in the direct shear tests (the AB plane as shown in Figure 8) was the least for this angle, of the seven different conditions of bedding plane inclination angles. Therefore, the shear resistance was primarily controlled by the intact rock in the direct shear test, which was clearly proved by the fracture pattern observed after the direct shear test, as shown in Figure 15.

With an increasing bedding plane inclination angle, the shear resistance was primarily controlled by the bedding planes during the shear process and slip along the bedding planes created the main fracture pattern, resulting in the decrease in the equivalent cohesion .

After the inclination angle θ became greater than 60°, the cracking or slipping of the bedding planes would not cause the overall failure of the specimens but made the specimens behave as a number of sheets parallel to each other. This resulted in an increase in the equivalent cohesion .

Furthermore, the trend of and inclination angle curve was similar under various bedding plane spacings. In addition, the decrease in the bedding plane spacing caused the equivalent cohesion to decrease. At different inclination angles the degree of cohesion weakening was approximately same.

5.2. Influence on the Shear Stress-Shear Displacement Curves

Due to space limitations, only the shear stress-shear displacement curves under a normal stress of 2.0 MPa were discussed here, as shown in Figure 14.

The direct shear tests of the specimens with inclination angle θ = 0° were completely dominated by the bedding planes. Thus, the trends of the curves were similar under different bedding plane spacings. However, the direct shear peak strength increased slightly with the decrease in the bedding plane spacing. This may be because the decrease in the bedding plane spacing caused the increase in the number of bedding planes. Thus, the total normal deformation increased under the same normal pressure, causing a slight reduction of the height of the shearing box. In the numerical calculation, the shear stress was calculated by dividing the force on the loading wall (wall 4 in Figure 8) by half of the height of the shearing box, which was slightly greater than the shear peak strength.

In general, with the decrease in the bedding plane spacing, the number of bedding planes increased, and the direct shear peak strength decreased. Therefore, the existence of the bedding planes had a weakening effect on the strength characteristics of the specimens.

In addition, the specimens showed different degrees of brittleness under different bedding planes geometries. When θ = 15°, the curves of all three bedding plane spacings were similar, with a similar shear peak strength and trend. All the shear stress-shear displacement curves had a sharp decline after the shear peak strength, which suggested strong brittle characteristics. This is also because the shear strength was mainly dominated by the intact rock, and a thoroughgoing fracture plane was easily formed during the shear process, resulting in a sharp decline in the shear strength.

With intermediate inclination angles (θ = 30° to 60°), the shear peak strength for the three bedding plane spacings varied greatly, indicating that the weakening effect of the bedding planes on the shear strength was becoming significant. When the bedding plane spacing was 10 mm or 20 mm, the shear stress-shear displacement curves remained stable after the peak strength and indicated fairly ductile behavior. This behavior was most clearly defined at θ = 45°. Even for the specimens with a single bedding plane, the curves had approximately no stress drop.

At θ = 75° or 90°, the curves of the three different bedding plane spacing tests had a small stress drop, and a large residual strength remained with the bedding plane spacing of 10 mm or 20 mm. In addition, the lower bedding plane spacing (10 mm) had a higher residual strength, which was contrary to the results of the intermediate inclination angles (θ = 30° to 60°).

5.3. Influence on the Fracture Pattern

The fracture patterns of the specimens under a normal stress of 2.0 MPa are showed in Figure 15, in which the red dots indicate the fractures from the SJM and the black dots indicate the fractures from the parallel bonds, corresponding to the fracture in the bedding planes and intact rock, respectively. Figure 16 shows the fracture pattern of the physical direct shear tests with the bedding plane spacing of 20 mm which was roughly the same as the results of the numerical models and verified the numerical method.

As shown in Figure 15(a), for the specimens with a single bedding plane, no thoroughgoing fracture plane formed when the inclination angle was θ = 45° or 60°. When θ = 45°, the specimens only broke in the lower left and upper right sections. In addition, the middle of the specimen was intact, resulting in a stable shear stress state without a stress dropping. When θ = 60°, the specimens still did not form a thoroughgoing fracture plane, but fracturing along the preexisting shear plane was observed. Therefore, the corresponding shear stress-shear displacement curves showed a stress drop, followed by a large postpeak residual strength. For the other inclination angles (θ = 15°, 30°, 75°, and 90°), all the specimens formed a thoroughgoing fracture plane in the shear process. Thus, the stress drop was nearly complete, and the residual strength was close to zero, as echoed in the shear stress-shear displacement curves presented in Figure 14.

The fracture patterns were similar for the bedding plane spacings of 10 mm and 20 mm. When the inclination angle was θ = 15°, a smooth connected fracture plane formed, making the postpeak residual strength small and supporting the analysis described above.

With intermediate inclination angles (θ = 30° to 60°), the specimens first broken along the lower left and upper right corners perpendicular to the bedding planes, starting from the A and B points of the loading walls (as seen in Figure 8). Next, the bedding planes failed, turning the specimen into a number of parallel sheets. In addition, the shearing force of the shear box on the sheets could be decomposed into a squeezing effect and slipping effect. The specimens could still withstand a certain shearing load, which was reflected in ductility in the shear stress-shear displacement curves.

With increasing inclination angles, the fractures in the lower left and upper right sections were propagating closer to one another, and a thoroughgoing fracture plane eventually formed, as shown in the fracture patterns at θ = 75° and 90°. This was reflected in the shear stress-shear displacement curves by a large stress drop after the peak strength. However, unlike the fracture pattern when the inclination angles were θ = 0° and 15°, the fractures in the specimens created a thoroughgoing fracture plane that was uneven and saw-toothed. Thus, the specimens maintained a certain residual strength, as seen in the shear stress-shear displacement curves. Furthermore, with the decrease in the bedding plane spacing, the saw-toothed nature of the thoroughgoing fracture plane was enhanced, resulting in a greater postpeak residual strength. This explained the analysis of the shear stress-shear displacement curves when θ = 75° and 90°.

6. Discussion

The actual rock masses rarely included the regular crystal or grains. And the irregular grain size has a certain impact on the strength of the numerical sample. The irregular grain will form a certain degree of interlocking effect. The more irregular the grain is, the greater the interlocking effect will be, which will result in the increase of the shear strength of the numerical samples. The “clump” model was used to replace the regular circle balls to assemble the numerical models which could create any different irregular elements we want to carry out simulation. With the same microscopic parameters, the numerical sample assembled with simple clump model was conducted by direct shear tests under different normal stress. The clump model and clump template were seen in Figure 17.

The results of the direct shear tests with irregular clump model were compared with the circle ball model. The strength curves of direct shear tests and the results of the linear fitting were shown in Figure 18 and Table 5, separately.

It can be seen that, under the same normal stress, the direct shear strength of the clump model was obviously increased, compared with the ball model. And, from the result of linear fitting, we can know that the irregular grain has small effect on the cohesion of the numerical model, with a slight increase. And it has a great influence on the internal friction angle φ.

In this paper, the macroscopic shear mechanics properties of the actual rock-like materials can be well fitted by the circle ball assembly. The rationality of the numerical models was verified by the comparison with the direct shear test results of the rock-like materials in the laboratory. However, how to show the superiority of the irregular grains relative to regular particles and which kind of irregular grains can better simulate the strength of the actual rock mass need to be further studied.

7. Conclusions

In this study, the bonded-particle element method with an improved embedded smooth joint model was applied to simulate the direct shear tests of transversely isotropic rock mass.

The sensitivity analysis of the microscale parameters of the SJM was carried out to determine the qualitative or quantitative relationships between the microscale parameters of the SJM and the macroscale mechanical characteristics of the bedding planes. In addition, the calibration sequence and method of the SJM parameters were designed. To better simulate the jointed rock mass, the SJM modeling method was improved. In addition, the physical experiments of the artificial transversely isotropic rocks were compared to validate the reason of the calibration method and modeling method.

Numerous numerical model tests were carried out to study the influence of the geometric characteristic of bedding plane geometry on the direct shear mechanical strength parameters, shear stress-shear displacement curves, and fracture patterns in the transversely isotropic rock masses.

The influence of the geometric properties of the bedding planes on the equivalent internal friction angle of the transversely isotropic rocks can be divided into two stages. When the inclination of bedding planes is less than approximately 30°, the spacing of the bedding planes has a small effect on . Otherwise, decreases with the decrease in the bedding planes spacing. With the decrease in the bedding plane spacing, the equivalent cohesion also decreased for various inclination angles. Furthermore, under same bedding plane spacing, as the inclination angle increased, first increased and then decreased and then increased; the maximum value of was reached at approximately θ = 15° and the minimum value of was reached at approximately θ = 45°.

The shear failure of transversely isotropic rocks with various geometrical characteristics can be divided into two patterns. The first pattern is a thoroughgoing fracture plane near the preexisting shear fracture plane, which creates a stress drop after the stress peak in the shear stress-shear displacement curve, when the inclination angle is near either 0° or 90°. The second pattern is damage perpendicular to the direction of the bedding plane near the both sides of the loading surface. In addition, when the inclination angle has intermediate values (θ = 30° to 60°), the intact rock is broken into many parallel sheets, which were squeezed and slipped past each other during the direct shearing process. The shear stress can maintain a stable residual strength after the peak strength.

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 they have no conflicts of interest.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (41272342, 51409013, and 41772308).