Abstract
In order to study the failure mechanism and assess the stability of the inlet slope of the outlet structure of Lianghekou Hydropower station, the strength reduction method considering the ubiquitous joint model is proposed. Firstly, two-dimension numerical models are built to investigate the influence of the dilation angle of ubiquitous joints, mesh discretization, and solution domain size on the slope stability. It is found that the factor of safety is insensitive to the dilation angle of ubiquitous joints and the solution domain size but sensitive to the mesh discretization when the number of elements less than a certain threshold. Then, a complex three-dimension numerical model is built to assess the stability of the inlet slope of the outlet structure of Lianghekou Hydropower station. During the strength reduction procedure, the progressive failure process and the final failure surface of the slope are obtained. Furthermore, the comparison of factors of safety obtained from strength reduction method and analytical solutions indicates that the effect of vertical side boundaries plays an important role in the stability of jointed rock slope, and the cohesive force is the main contribution to the resistant force of vertical side boundaries.
1. Introduction
The stability of jointed rock slopes for mining, hydraulic, and hydroelectric engineering is crucial for the safety of engineering work [1]. To prevent and reduce the occurrence of instability of rock slope, the stability and deformation of rock slopes during and after excavation at different scenarios are analyzed [2–12]. For slope stability analysis, the limit equilibrium method (LEM) is widely used by engineers and researchers. Usually, slope stability was evaluated by LEM in two dimensions [13]. However, the 2D limit equilibrium analyses simplify the problem into plane strain conditions without considering the true three-dimensional characteristics of slope [14]. Therefore, many researchers improved a number of 3D limit equilibrium methods for slope stability analysis based on extensions of corresponding methods [3, 15–21], which are beneficial for stability assessment of slope with complex failure surfaces.
Although researchers made many improvements to the 3D limit equilibrium method, there are still some limitations when analyzing the stability of rock slopes. Firstly, the limit equilibrium method did not take into account the stress and deformation of the rock slope. Secondly, the failure surface of the slope was predefined by engineers. Thirdly, many assumptions on the internal force distributions were used to simplify the governing equations to solve the factor of safety. In addition, the evolution process of the failure surface could not be simulated using LEM.
The strength reduction method (SRM) based on the finite element method or finite difference method can overcome those limitations of LEM in the stability analysis of slope. This method was firstly proposed by Zienkiewicz et al. [22] to analyze the stability of the slope. Recently, many researchers [23–27] improved this method. Cheng et al. [28] compared SRM with LEM and found that generally, the two methods will give a similar factor of safety (FOS). Tschuchnigg et al. [29] further discussed the influence of plasticity flow rule and mesh discretization on the FOS. The SRM method was also applied to assess the stability of practical rock slope engineering. Using the strength reduction FEM, Jiang et al. [30] analyzed the stability of the left abutment slope of Jinping-I hydropower station during the period of construction and assessed the effect of shear-resistance tunnels. Zhang et al. [31] proposed a method combining fracture mechanics and the SRM to simulate the fracture behavior in a jointed rock slope and assessed the slope stability. Schneider-Muntau et al. [32] introduced a strength reduction method for the barodetic material model and presented the comparison of the results of slope stability calculated using an elastoplastic material model (Mohr-Coulomb) with the results of an analytical analysis according to Bishop are presented. Nevertheless, most of these studies focused on either homogeneous soil or rock slopes with predefined failure surface. As a natural geological material, in actuality, rock mass contains various randomly or nonrandomly geological discontinuities of different sizes, persistence, and toughness. Moreover, different from the circular failure surface of homogeneous soil slope, the rock slope is prone to failure along the discontinuities. Generally, the failure surface of the rock slope is constituted with located discontinuities, such as faults, and random discontinuities, such as joints. In this case, the failure surface of rock slope cannot be predefined and part of it should be searched during the process of strength reduction. In this study, the strength reduction finite element analysis considering ubiquitous joint model is adopted to assess the stability of jointed rock slope. Firstly, the strength reduction method considering the ubiquitous joint model are presented. Then, one example of 2D slope containing a set of ubiquitous joints is adopted to investigate the effect of mesh discretization, dilation angle, and domain size on the failure mechanism and the stability of slope. Finally, the proposed strength reduction method is adopted to analyze the failure mechanism and the stability of the inlet slope of outlet structure of Lianghekou Hydropower station.
2. Methodology
2.1. The Ubiquitous Joint Model
The ubiquitous joint model [33] is a modified Mohr-Coulomb model considering the presence of the orientation of weak planes. The criterion for failure on the weak planes, whose orientation is given, consists of a composite Mohr-Coulomb envelope with tension cutoff. The position of a stress point on the latter envelope is controlled again by a nonassociated flow rule for shear failure and an associated rule for tension failure. The procedure of implementing the ubiquitous joint model is as follows: Step 1: Calculate the stresses using the Mohr-Coulomb model. The stresses corresponding to the elastic guess for the current step is analyzed firstly for general failure, and relevant plastic corrections are made using the Mohr-Coulomb model. The resulting stress component, labeled as is then examined for failure on the weak plane. Step 2: The resulting stress is then examined for failure on the weak plane. If the stress violates the composite yield criterion along the weak plane, then the tensor of stress corrections is calculated and added to to yield new stress values, denoted as . Step 3: If the stress is in the envelope of composite yield criterion along the weak plane, no plastic flow along the weak plane takes place, and the new stress is equal to .
The detailed procedure for implementing the ubiquitous joint model refers to the section of constitutive models in the manual of FLAC software.
2.2. Strength Reduction Principle
In the strength reduction method, the factor of safety (FOS) is defined as the ratio between the actual shear strength and the reduced shear strength for the fault, joints, and intact rock when the slope arrives at a critical state. When implementing the strength reduction procedure, the reduced shear strength parameters, cohesive force , and the fiction angle are obtained by the following:where SRF is the strength reduction factor (SRF), c and φ are the actual shear strength parameters. In the strength reduction procedure, the reduced shear strength decreases gradually until the slope arrives at the critical state, and the corresponding SRF is equal to the FOS.
The solution procedure of the SRM considering the ubiquitous joint model is as follows: Step 1: An initial value of the SRF is selected and the reduced shear strength parameters of faults, joints, and intact rock are calculated by equation (1). Step2: Implement numerical calculation for the slope using reduction shear strength parameters. If the slope fails, the SRF is regarded as the FOS. Step3: The SRF increases with an increment gradually. Go back to step 2 until the slope fails.
2.3. Failure Criterion of the Slope
There are three widely used failure criteria for assessing the stability of the slope, namely, numerical nonconvergence criterion [24, 34], the criterion of plastic yield zone connection [35–37], and displacement mutation criterion [38, 39], respectively. Regarding the first criterion, the failure state of the slope is determined by the convergence condition. However, for the stability analysis of a complex 3D slope model, the convergence condition is influenced by many factors, such as element type, yield condition, constitution model, and convergence precision, which may terminate the computation early [40]. The second criterion states that when the plastic yield zones coalesce from the foot to the top of the slope during strength reduction, the slope will be in a failure state. For the third criterion, the displacement of some monitoring points on the slope surface is firstly recorded during the strength reduction procedure. The inflection point of the curves of displacement versus the SRF is considered as a failure state of slope and the corresponding SRF is equal to the FOS. Compared with the abovementioned criteria, the displacement mutation criterion is selected in the following analyses because of its clear physical concept and numerical tractability.
3. Verification and Analysis of Influencing Factor
3.1. Problem Description and Mesh
A rock slope with a slope angle of 45° is selected a case study, as shown in Figure 1. In the slope, ubiquitous joints with a dip angle β = 30° are taken into account. The slope dimensions are shown in Figure 1. The mechanical parameters for rock and the joint materials are shown in Table 1. It is noted that the ubiquitous joints drew in Figure 1 are schematic. Hence, the spacing between lines does not represent the actual spacing of joints. The mesh discretization for this slope is shown in Figure 2, which contain 1875 zones and 1976 grid points.


3.2. Influence of Dilation Angle of Joints
The definition of the dilation angle of joints is an important issue when using the strength reduction technique considering ubiquitous joints. To study this issue, five different dilation angles are adopted to investigate the effect of flow rule, namely, 0°, 5°, 10°, 15°, and 20°. The horizontal displacement of a monitoring point at the crest of the slope, shown in Figure 2, is recorded during the strength reduction procedure. Figure 3 shows the curves of horizontal displacement of the monitor point versus the strength reduction factor (SRF) for the cases with different nonassociated flow rules. One can observe that the horizontal displacement increases with the increase of the SRF. Moreover, the curves obtained by five dilation angles are approximately the same when the SRF is no more than 1.26 but separated from each other when the SRF is larger than 1.26. According to the displacement mutation criterion, the FOS of those five dilation angles is the same and is equal to 1.26. The results indicate that the FOS is insensitive to the dilation angle of the joints.

Figure 4 illustrates the evolution of incremental shear strains during the SRM when the dilation is equal to 20°. Since the shear strain rate around the failure surface is much larger than other regions of the slope, in the following analysis, the center of the zone of incremental shear strain can be defined as the failure surface of the slope. It can be seen from Figure 4 that the failure surface firstly occurs at the toe of the slope. With the increase of the SRF, the failure surface propagated toward the crest of the slope gradually and passed through the whole slope when the SRF arrived at 1.26, as shown in Figure 4(c).

(a)

(b)

(c)
3.3. Influence of Mesh Discretization
It is well known that the element type, the mesh discretization, and the convergence tolerations have a significant impact on the FOS. In order to study the mesh dependency, the same model as shown in Figure 1 is used, but nine types of mesh discretization from coarse to very fine are defined. Their number of elements are 169, 300, 469, 675, 919, 1200, 1519, 1875, 2700, and 4219, respectively. From a practical point of view, the coarse and very fine mesh are not realistic, but these two extreme cases are included to show the effect of the mesh influence.
Figure 5 compares the curves of the FOS versus the horizontal displacement of the monitoring point for the cases with different mesh discretization. It is observed that the horizontal displacement increases with the increase of the SRF. Moreover, the curves obtained by different mesh discretization are approximately the same when SRF is no more than 1.24 but separated from each other when SRF is larger than 1.24. In order to investigate the effect of mesh discretization on the factors of safety, the curve of the FOS versus the number of elements is illustrated in Figure 6. It is seen from Figure 6 that the mesh discretization has a significant impact on the FOS. With the increase of the number of elements, the FOS decreases until the number of elements is equal to 2800, which indicates that the FOS is insensitive to the mesh discretization when the number of elements is over a certain threshold.


Figure 7 shows the contour of incremental shear strains obtained from three types of mesh discretization. It is found that the mesh discretization has a significant influence on the distribution of incremental shear strains. For a very coarse mesh, shown in Figure 7(a), the zone of incremental shear strains is wider and the failure surface does not pass through the toe of the slope. However, for a very fine mesh, shown in Figure 7(c), the zone of increment shear strain is narrow and the failure surface is clearly defined.

(a)

(b)

(c)
It also can be observed from Figure 7(c) that the central plane of the zone of incremental shear strains, namely, the failure surface, is a plane. It means that the failure mode of the slope is the plane failure mode that is somehow different from that of a homogeneous slope in which the failure surface is a circle failure mode. Therefore, there is an analytical solution to this 2D ubiquitous joints rock slope and the FOS is 1.20. As discussed above, when the mesh discretization is fine enough, the FOS obtain from SRM is 1.24, which is close to the analytical solution. Therefore, it is reasonable to assess the ubiquitous jointed rock slope using the proposed method.
3.4. Influence of Domain Size
In order to assess the size effect of the calculation domain in the strength reduction method, a length factor d, which quantifies the distance from initial boundaries to extended boundaries of the slope, is introduced and shown in Figure 8. In this study, six different cases of d are considered, namely 5 m, 15 m, 25 m, 35 m, 45 m, and 55 m. In the strength reduction procedure, the mesh discretization of the initial region of the slope is the same as that in Figure 2. Figure 9 illustrated the displacement of the monitoring point evolving with the SRF. It is clear that the curves obtained by six cases are nearly the same when the SRF is no more than 1.26 but separated from each other when SRF is larger than 1.26. According to the displacement mutation criterion, the FOS of those six cases is the same and the value is 1.26. It also indicates that the FOS is insensitive to the calculation domain size.


4. Engineering Application
4.1. Project Description
The Lianghekou Hydropower station is located on the middle reach of Yalongjiang River, 10 kilometers northern of Yajiang, Sichuan province, China. In this project, a 295 m height gravelly clay-cored rockfill dam is constructed to form a reservoir. In order to fulfill the flood discharging capability of the reservoir, four outlet structures are built, including a tunnel-spillway, deep hole flood discharging tunnel, emptying tunnel, and rotation-flow shaft flood discharging tunnel. The inlet slopes of these four outlet structures are excavated as a whole, which results in a massive engineering slope with 606 m height and 1200 m length, as shown in Figure 10. A typical geological profile 8–8 of the slope is shown in Figure 11. The rock of the slope is composed of deep grey and medium-thick layer metamorphism siltstone, silty slate, and sericite slate. This slope is a typical antidip layered slope and the orientation of the bedding plane is 200–~210°∠65~–75° (dip direction/dip angle). Besides the bedding plane, the other three sets of joints are developed in the slope and their orientation is 130~150°∠10~45°(J3), 90~120°∠70~90° (J4), 90~120°∠10~30° (J6), respectively. These joints are rigid structural plane without filling. In addition, as shown in the geological profile 8–8, a large number of faults are developed. Most of these faults are compressive faults which developed along the bedding plane, including f1, f2, f3, f4. Below the 3050 m platform, two main faults are detected, namely, Fault f1 and Fault f13-04. As shown in Figures 10 and 11, the F1, with the occurrence 180°∠40°~60°, runs through the entire engineering slope and outcrops above the inlet of the shaft flood discharging tunnel. The f13-04, with the orientation 10°∠30°, outcrops at the step between elevation 2790 m and elevation 2815 m. According to the spatial distribution of faults, a huge potential failure region sliding along the fault f13−04 is observed, where the fault f1 form the back boundary and the fourth set of joints (J4) developed ubiquitously form the two vertical side boundaries of the block. After the completion of the excavation, fault fb13−04 will be exposed on the excavated surface, which will deteriorate the stability of the slope drastically.


4.2. Description of Numerical Model and Material Property
Different from the displacement finite element analysis, the strength reduction method focuses on assessing the failure mechanism and the stability of slope rather than obtaining the distribution of displacement and stresses in the slope. Therefore, the potential failure region is the region of great concern during the strength reduction procedure. In order to simplify the numerical model, the terrain in the potential failure region shown in Figure 10 is consistent with the real terrain in the project site, while the terrain around the potential failure region is simplified. The three-dimension numerical model for strength reduction procedure is shown in Figure 12. In this model, the lowest left corner is used as the origin of the x-y plane, and the z-axis begins at an elevation of 2540 m. The size of the slope model is 840 m length in the x direction, 1000 m width in the y-direction, and 700 m height in the z direction. The numerical model is made up of 174944 gridpoints and 214823 zones. The two main faults (f1 and f13−04) and the fourth set of joints (J4) are considered in this model. The grid of potential failure region is shown in Figure 12. The mechanical parameters from the recommended values of the designer [41] are listed in Table 2.

4.3. Results
For the potential failure domain, the key block is composed of fault f1, fault f13−04, and the fourth set of joints (J4). It is noted that the fourth set of joints is ubiquitous in the slope and their locations are undetermined. The proposed SRM method is used to search for the locations of J4 forming the key block and determine the factor of safety of the key block. Two monitoring points M1 and M2, shown in Figure 12, are selected to record the displacement for different strength reduction factors. The curves of y-direction displacement of monitor points versus the strength reduction factor are illustrated in Figure 13. It can be seen that the y-direction displacement of three monitoring points is roughly proportional to the SRF. However, when the SRF increases to 1.36, the displacement curves of M1 and M2 start to deviate from the linear relationship between displacements and SRF. Therefore, SRF = 1.36 can be considered as a catastrophic failure point and the FOS of the slope.

The progressive failure process is investigated during the procedures of reducing the shear strength of the fourth set joints, Fault f13−04, and Fault f1. Figure 14 shows the evolution of the zone of incremental shear strain in the block at different SRF stages. Figure 15 illustrates the evolution of the plastic yield zone in the block during the process of strength reduction. In Figure 15, the yield state “None” refers to the elastic state of the zone and the “u: shear-n” indicates that the zone is in shear failure along the joints currently. As we can see from Figures 14 and 15, for the case that SRF is equal to 1.1, the zone of incremental shear strain and the plastic zone initiate at the left side of the outcrop of Fault f13−04 and the middle of the inlet of rotation-flow shaft flood discharging tunnel. As the SFR increases, the zone of incremental shear strain propagates along the strike of the fourth set of joints. When the SRF increases to 1.26, the plastic zone and the zone of incremental shear strain propagate through the inlet of the rotation-flow shaft flood discharging tunnel. When the SRF increases to 1.36, the plastic zone propagates though the whole potential failure region and becomes the vertical side boundaries of a new failure block.

(a)

(b)

(c)

(a)

(b)

(c)
As a comparison, the solution obtained by the LEM is also proposed. According to the failure surface obtained by the SRM, the vertical side boundaries and sliding surface of the failure block is illustrated in Figure 16. The areas of sliding surface, left side boundary, and right side boundary are = 10,389 m2, = 4909 m2, and = 605 m2, respectively. The volume of the failure block is V = 32,500 m³. If the resistant force of two vertical side boundaries is ignored, the FOS can be determined as follows: where θ is the dip angle of bottom sliding surface, , is cohesion and friction angle of bottom sliding force, is the unit weight of the potential failure block. The FOS obtained by equation (2) is 1.10, which is far less than the FOS obtained SRM. Obviously, the difference between the two safety factors is attributed to the effect of the two vertical side boundaries. This indicates that the two vertical side boundaries have a significant impact on the stability of the jointed rock slope, and ignoring the effect of the vertical side boundaries will lead to underestimation of the stability of the slope.

The resistant force of the vertical side boundaries consists of two parts, namely, the cohesive force and the frictional force acting on the vertical side boundaries. For the limit equilibrium method, it is impossible to obtain the magnitude and distribution of the normal force of the vertical side boundaries. Therefore, the FOS is calculated only considering the effect of the cohesive force of vertical side boundaries. Assuming that the direction of the cohesive force of the vertical side boundaries is along the sliding direction of the potential failure block, the resistant force can be modified as follows:where is the cohesion of the fourth set of joints. The FOS considering the cohesion force of the vertical boundaries is 1.27, which is 0.17 more than the FOS obtained in the case that the effect of vertical side boundaries is ignored and 0.09 less than the FOS obtained from the SRM. This result indicates that the cohesive force is the main contribution to the resistant force of the vertical side boundaries, while the frictional force has a slight influence on the stability of the slope.
5. Conclusions
The purpose of this paper is to provide a methodology to analyze the stability of the jointed rock slope by strength reduction method considering the ubiquitous joint model. In order to verify the proposed method, several two-dimension models are built to investigate the effect of the dilation angle of joints, mesh discretization, and boundary size of the numerical model on the stability of the jointed rock slope. The numerical results show that the factor of safety is insensitive to the dilation angle of the joints and the boundary size of the numerical model. For the influence of mesh discretization size, the width of the zone of incremental shear strain and the factor of safety are sensitive to the size of the element. When the number of elements increases to a certain threshold, the factor of safety will not be affected by the mesh discretization.
The stability of the inlet slope of the outlet structure of Lianghekou Hydropower station has been investigated using the proposed method. The stability of this slope is controlled by the huge potential failure region formed by fault f13−04, fault f1, and the fourth set of joints (J4). The final failure block is searched during the strength reduction procedure and the factor of safety is 1.36. In order to investigate the effect of vertical side boundaries, the factor of safety of the final failure block is calculated by the LEM in two different cases. The results show that the effect of the vertical side boundaries plays an important role in the stability of the jointed rock slope and the cohesive force is the dominant contribution to the resistant force of vertical side boundaries.
Data Availability
The 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
The work reported in this paper has received financial support from National Key R&D Program of China (nos. 2019YFC0605001 and 2018YFC1505005), the National Natural Science Foundation of China (no. 42077253) and the Guiding Project of Scientific Research Plan of Hubei Science and Technology Department (No. B2020408). These supports are gratefully acknowledged.