Abstract

The stability of rocky slope is determined mainly by discontinuities. The discrete element calculation method can be used to analyze the geometric features of rock structures and to deal with the nonlinear deformation and destruction of the rock mass that may be affected by the discontinuities. In this paper, the Barton–Bandis (B-B) nonlinear strength criterion was introduced to improve the Mohr–Coulomb (M-C) slip model of joint slope. The modified model could reflect the real-time shear strength in a changing stress state. Using the numerical calculation of the shear test of the structural plane, we found that the corrected slip curve fits well with the process before the failure of the shear test. Furthermore, the modified model can track the disintegration between blocks while sliding failure of joint slope. With an increase in the number of structural planes and the complexity of the relative location of blocks, the interaction force between blocks and the sliding failure mode of the joint slope would be more complex. Changing the nonlinear parameters with the stress state of the structure plane could effectively solve this problem.

1. Introduction

With the improvement and renewal of the application of computer technology [13], numerical calculation has become an important method to effectively solve geotechnical engineering problems [4, 5]. Specifically, the test data obtained from the experiments [69] and the related constitutive models [10, 11] are important factors to determine whether the numerical simulation is effective.

As for the rocky slope, its stability is determined mainly by the discontinuities [1214]. The discrete element calculation method can be used to analyze the geometric features of rock structures and to deal with the nonlinear deformation and destruction of the rock mass that may be affected by the discontinuities [15, 16]. For example, Wu and Hsieh [17] adopted the discrete element calculation method to analyze the debris movement and deposition of the Chiu-fen-erh-shan landslide under the impact of the earthquake. Liu et al. [18] investigated the initiation and kinetic characteristics of the Xiaogangjian rockslide failure via the three-dimensional discrete element method. Fan et al. [19] studied the jointed rock mass stability based on the tunnel seismic prediction and discrete element method.

It is essential to discuss the structural plane shear strength for the joint rocky slope. Scholars tend to use the Mohr–Coulomb (M-C) slip model when conducting discrete meta-analysis and to treat the parameters of the slip model as constants (e.g., [2022]). However, the description of slope slip failure is a significant limitation in this setting method.

This limitation of the M-C slip model could be improved effectively by introducing the Barton–Bandis (B-B) nonlinear criterion. The B-B model empirical formula was proposed by Barton and Bandis based on a large number of shear tests of natural structures. For rock engineering, the B-B model was widely used to analyze and infer the shear strength of rock joints. Recently, many experts and scholars had studied and discussed the application of the B-B model in the discrete element. By combining the B-B empirical formula with discrete elements, shear slip numerical tests [23] and double-axis compression tests [24] have been used to analyze discontinuities. Existing research studies combined M-C parameters with the empirical formula of Bandis [25] and discussed the contact problems between noncontinuous deformation analysis and the Bandis empirical formula [26]. The connection between these two models was achieved effectively using the discrete element method.

In the process of introducing and revising the B-B shear strength criterion, however, most scholars have not set the stress state of continuous change in the actual slip damage as the key point to the analysis of the shear strength of the rock mass. As a result, it is difficult to introduce nonlinear principles into actual working conditions. To apply the modified slip model to practical engineering, this paper conducts a new exploration and provides an innovation.

In this paper, by combining a discrete element numerical calculation using the B-B empirical formula with the structural plane shear test, we applied the corrected B-B parameters to the M-C parameters in the slide model. The shear strength parameters of the structural plane changed with the stress state, which we realized by FISH language programming. We then analyzed the slope slip failure before and after modification of the M-C slip model.

2. Shear Strength Transformation of B-B Formula

To obtain the connection between shear strength parameters and the stress state, we transformed and analyzed the B-B empirical formula. For the shear deformation of the structural plane, the curve before the peak intensity was satisfied as follows [27]:where is shear stress and is shear displacement. In addition, m =  and is initial shear rigidity, n = , and is the asymptotic value of the hyperbolic curve of .

The shear stiffness is the stress gradient in the unit deformation, which could be obtained as follows:where is the shear stiffness, is the normal stress, is the shear stiffness index, is the failure coefficient, and .

The modified function could be adapted to the empirical formula as follows [28]:where is the peak shear strength. This formula could be used to characterize variations in the shear stiffness of joints.

We derived formulas (1)–(3) through the shear test data. The shear strength curves of the structural plane can be obtained by the shear test of the structural plane (Figure 1). In this shear test, the normal stress was set to 50 kPa, 100 kPa, 150 kPa, and 200 kPa. Because of limited space, the experimental operation of the shear test of the structural plane is not included in this paper.

According to formulas (1)–(3), the related parameters could be extracted from the shear test data. By processing the data, we obtained a correlation formula and its linear fitting curve, as shown in Table 1 and Figure 2.

When we substituted the formulas in Table 1 into formula (3), the shear stiffness characterization formula based on the Bandis empirical formula could be obtained as follows:

Then, the normal stress and the normal deformation satisfy the following formula [28]:where p and q are constants. The normal stiffness could be obtained as follows:

By fitting results from the direct shear test of the structural plane , the formula could be obtained as follows:

By substituting formula (7) into formula (6), under a given normal stress , the structural plane stiffness could be obtained as follows:

3. The Parameter Transformation of the B-B Empirical Formula

To introduce the nonlinear parameters of the B-B model into the M-C slip model, we analyzed the nonlinear empirical formula. This formula could be used to estimate the peak shear strength of irregular and unfilled structure [29, 30]:where is the effective normal stress (MPa) and is the structural plane roughness coefficient that can be obtained by comparing it with the typical curve given by Barton. is effective compressive strength for the structural plane (MPa); if the structural plane does not have any adjacent rock weathering, is equal to the rock uniaxial compressive strength; otherwise, can be equal to a quarter of the rock’s uniaxial compressive strength; is the basic friction angle of the rock’s mass structural plane (°), which can be obtained by laboratory test. In the empirical formula, the combination could be defined as the roughness angle, which was related to the structural plane rough rolling form (°). In this paper, the values of the parameters , , and were obtained according to the literature [31].

The linear M-C formula could be written as follows:

Also, taking the derivative of on both sides, the formula would be changed as follows:

Then, taking the derivative of on both sides, formula (9) could be changed as follows:

Substituting formula (13) into formula (11), the equivalent friction angle and the equivalent cohesion value could be obtained as follows:

By referring to the relevant literature and combining this information from the literature with the actual situation of the structural plane, we took the following parameters: . We defined as . After applying these data, the formula could be obtained as follows:

By combining the B-B nonlinear relationship with the M-C surface slip model, we could obtain the normal stiffness, tangential stiffness, cohesion, and friction angle of the M-C sliding model through the normal stress of the structural plane.

4. Modification of Empirical Formula

To verify the rationality of formulas (4), (8), and (15), the shear test was numerically simulated by 3D Distinct Element Code (3DEC) software. As for the 3DEC model, we adapted the M-C constitutive model as the basis, adopted the M-C slip model as the structural plane, and then modified and replaced the strength parameters with the B-B formula. The test specimen model was simplified to a 60 mm 60 mm 120 mm rectangle, and the joint plane was the half-height of the specimen. The model and boundary conditions are shown in Figure 3, and the physical and mechanical parameters of materials are given in Table 2.

In the process of numerical tests, while maintaining the normal stress as 0.05 MPa, 0.1 MPa, 0.15 MPa, and 0.2 MPa, the shear stress uses 0.5 and under different normal stresses for the numerical calculation. In addition, when , then 0.25 and 0.75 were taken to provide basic data for the fitting of the modified parameters. We monitored the shear displacement of the feature point (0, 0, 0) on the joint surface in the numerical experiments. By comparing the numerical simulation data with the experimental data, the results of the numerical test tended to be close to the indoor test results only when . In other settings, especially when approaches , the simulation results would produce a significant difference. Therefore, we should revise the B-B empirical formula.

By changing the shear stiffness value in the numerical calculation, we obtained the shear deformation consistent with the indoor test, and then we also could obtain the corresponding value of different shear stress under  = 0.2 MPa. On the basis of these shear stiffness values, we introduced the logarithmic function to modify the B-B empirical formula [30]. The corrected function and the modified joint shear stiffness could be obtained as follows:

The data referred to in this section are given in Table 3.

As shown in Table 3 and Figure 4, the calculated displacement with the modified Bandis formula was basically consistent with the measured shear displacement. Therefore, the revised formula could be applied to the calculation of shear strength. For the analysis of structural surface deformation based on 3DEC, it would be difficult to realize the change of shear stiffness of the structural plane with the change of tangential stress . So, we calculated obtained in the condition of as the reference value of shear stiffness. Then, the shear stiffness value could be approximately simplified as a function of normal stress, as follows:where D is constant. Formula (17) could simulate the shear deformation characteristics of the structural plane before shear failure. Using this method, we realized the deformation characteristics of the semiquantitative analysis with few errors, and it was easy to implement in the numerical simulation.

By reintroducing these formulas to the numerical calculation, we changed the load method to the shear test (Figure 5) [32]. Specifically, compared with the shear test method shown in Figure 3, a fixed tangential loading speed (0.005 m/s) was used to complete the shear test shown in Figure 5. Also, the related shear numerical test curves can be obtained (Figure 6).

By comparing Figure 1 with Figure 6, we found that the shear strength curve obtained by introducing the B-B empirical formula was basically same as the measured curve. Additionally, the numerical test curve could effectively simulate the shear test curve in the process of changing from the ascent phase to the peak phase. Therefore, these formulas could be used to analyze the joint slope model for sliding failure by 3DEC.

5. Comparison and Analysis of Slope Failure

5.1. Joint Slope Model

To verify the applicability of calculating formulas, we used two sets of controlled structural planes to carry out the numerical calculation of the sliding failure. Among these planes, the joint’s strike was 0° and dip was 80°, and the rock layer’s strike was 0° and dip was 20°. We numbered and monitored the sliding block and joints in the calculation (Figure 7). The right and bottom boundaries of the joint slope model are constrained by the three-direction displacement.

According to previous discussion, by combining formulas (15) and (16) and making the appropriate reduction, we could obtain the empirical formula of the shear stiffness of the limestone joints in the condition , as follows:

Using 3DEC built-in FISH language programming, we automatically obtained the normal stress of the structural plane, and then we could calculate the cohesive force, friction angle, shear stiffness, and normal stiffness by formulas (8), (14), and (17). We automatically reset the shear strength parameters of the structural plane every 50 steps, so that the strength parameters of the structure could be constantly adjusted in the process of slippage. The related parameters in numerical calculation are listed in Table 4.

5.2. Comparison of Slip Modes

We carried out a rigid sliding failure analysis on the slopes and extracted the displacement vector diagram of the sliding zone every 200 steps until the slip block had detached from the slope (Figure 8).

According to a comparison and analysis of the slide vector graphs, when the shear strength was constant, the slip velocity of the block was consistent, and no separation occurred between blocks in the process of slippage. Relative to the corrected slip model, the start of the slide process was slower and the displacement volume was relatively smaller. The slippage disruption model was too ideal to describe the sliding failure, and its predictions of disruption were more conservative than in reality.

After introducing the B-B formula to modify the M-C slide model, the structural strength of each block was different because of their stress state. In the process of slippage, the slip rate of the blocks was different, and the tendency for disintegration and destruction between blocks was evident. Because the sliding failure mode was similar to the actual situation, it also proved the rationality of introducing the B-B model.

The monitoring of shear strength of the structural plane is shown in Figure 9. By means of real-time monitoring of the displacement of four blocks in the slip area, we obtained the calculated step-slip displacement curve, as shown in Figure 10.

Figure 9 shows that the normal stiffness and shear stiffness of the structure plane constantly changed during the sliding process. When the stress state changed, the normal stiffness and shear stiffness of each structure would change synchronously. At the steps from 0 to 300, there was an adjustment period in the stress, which corresponded to the glide phase of the displacement monitoring (Figure 10(c)). Because the initial control condition of the slip model was rock mass gravity, the shear strength for different spatial location was significantly different. Among these differences, the peak positive stress of the horizontal joint was the maximum value and its changing amplitude was the largest, whereas the value of vertical joint was relatively small.

The monitoring curves (Figure 10) show that, under the unmodified M-C slip model, the slippage rate of each block remained basically consistent, and with an increase in the calculation step (up to 800 steps), the slip rate of the block had only a small difference. In the process of modified slip, for No. 2 block, there was free surface along the slide direction, and the block was pushed by the trailing edge block; thus, its slip rate was the fastest. No. 3 block was in the trailing edge of the sliding zone, and its sliding direction was restricted by the front block; thus, its slip rate was the slowest. As for No. 1 block and No. 4 block, because of the influence of the surface and the velocity of the bottom block, the sliding speed of these blocks remained basically the same.

Taking No. 2 block as an example, we compared the slip rate under two slide models as shown in Figure 10(c). After we introduced the B-B formula, the curve pattern of the slip rate did not change, and both curves could be fitted to quadratic polynomials. In terms of slip rate, however, the joint shear strength with the uncorrected M-C slip model did not change in the process of sliding. Therefore, its slide process was ideal and its calculations were more conservative than in reality. This trend was consistent with the block slip phenomenon observed from the displacement vector diagram of slip damage.

6. Discussion and Conclusion

By introducing the B-B empirical formula to the M-C slip model, we corrected the M-C parameters. The modified model could reflect the real-time shear strength in a changing stress state. It should be noted that this method focuses on the peak strength of shear failure and does not consider the residual strength. Using the numerical calculation of the shear test of the structural plane, we found that the corrected slip curve fits well with the process before the failure of the shear test and that it could be applied to the test situation.

According to the rigid sliding failure analysis of the joint slope, the slope had been in an overall slip state with the unmodified M-C model, and the internal block displacement and the rate of the sliding area remained basically consistent. After we modified this model by introducing the B-B empirical formula, we could transform the stress state difference caused by the structural plane and the thickness of the overburdened rock layer into the nonlinear change of shear strength for each structural plane. In the sliding process using the modified M-C model, we noted the tendency for disintegration to occur between blocks, which was close to the reality.

Whether or not we had modified the M-C model was not essential to the curve form of the individual block slip rate, but the slip process was ideal, and the calculation results were overly conservative because the shear strength of the unmodified M-C slip model was constant during the process of slippage.

In practice, with an increase in the number of structural plane and the complexity of the relative location of blocks, the interaction force between blocks and the sliding failure mode of the joint slope would be more complex. Therefore, the M-C slip model would not be able to accurately reflect the sliding failure mechanism of the joint slope. Changing the nonlinear parameters with the stress state of the structure plane could effectively solve this problem.

Data Availability

The test data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the China Postdoctoral Science Foundation funded project (2021M700608), the Hebei Key Laboratory of Earthquake Disaster Prevention and Risk Assessment (FZ213203), the Postdoctoral Innovative Talents Support Program, Chongqing (CQBX2021020), the Natural Science Foundation of Chongqing, China (cstc2021jcyj-bsh0047), the Hubei Provincial Natural Science Foundation of China (2020CFB352), and the follow-up of the Geological Disaster Prevention and Control Project in the Three Gorges Area (0001212019CC60001 and 0001212021CC60001).