# Study on Mechanical Response and Stability Algorithm of Soft and Hard Rock Interbedded Slope Excavation

Authors: Qingyang Ren; Xin Meng
Journal: Complexity
 (2022)

---

## Abstract

Based on the actual engineering geological survey, the mechanical response mechanism and sliding surface characteristics of the interbedded slope containing six layers of soft and hard rock during excavation under heavy rain conditions are studied by numerical simulation, and the composite sliding surface stability algorithm of the layered slope is established by using the penalty function. The research shows that, during the excavation, the distribution characteristics of stress, strain, plastic zone, and horizontal displacement along the bedding plane are obvious; the stress is concentrated in the hard rock layer; however, the strain and plastic zone are concentrated in the soft rock layer. When the ratio of the maximum vertical displacement to the maximum horizontal displacement is less than 1.85, the slope will slide. According to the relative position of the soft rock layer and the toe of the slope, the sliding surface of soft and hard rock interbedded slope can be divided into two types, namely, circular-straight sliding surface and circular-straight-circular sliding surface. The verification shows that the composite sliding surface algorithm objectively reflects the destruction of the slope by bedding sliding and can reasonably evaluate the safety of soft and hard rock interbedded slope.


---

## Body

## 1. Introduction

The geological structure and mechanical properties of soft and hard rock interbedded slope are very complex, and the stability problem during excavation is prominent. In order to prevent the instability of this kind of slope excavation, the mechanical response law and stability calculation method during the excavation process should be mastered. In recent years, scholars have paid more and more attention to layered slopes. In the aspect of bedding slope stability, Cao et al. revealed the influence of various human factors on the three-dimensional stability of soft rock slopes with weak bedding during the excavation process [1]. Xue-chi et al. studied the deformation law and stability of bedding rock slope under secondary excavation [2]. Li et al. revealed the control effect of double weak interlayers on slope stability [3]. Zhang et al. studied the deformation mechanism and stability of bedding soft rock slopes [4]. Liu et al. proposed that the stability of rock slopes can be evaluated by analyzing the plane properties of in situ bedding [5]. In terms of deformation and failure of bedding slopes, Fan et al. proposed that the failure modes of bedding slopes are mainly vertical tension cracks at the rear edge of the slope, bedding sliding of rock layers along the interlayer, and rock block caving at the top of the slope [6]. Zhang et al. analyzed the deformation and failure mechanism of bedding slopes with gentle dip angles [7]. Xi et al. studied the failure characteristics of bedding rock slopes with multiple weak interlayers [8]. Zhou and Tang studied the sudden slip failure mechanism of bedding slopes based on the energy principle [9]. Li and Cheng studied the creep characteristics of bedding landslides [10]. In terms of the sliding surface of bedding slopes, Zhou et al. found that the strongly weathered soft and hard rock interbedded slopes form a fold-line potential sliding surface from joints, fissures, and rock layers, and there is more than one potential sliding surface [11]. Ding et al. studied the deformation and failure process of typical bedding mining slopes and the evolution law of rock and soil properties on the sliding surface [12]. Liu et al. used a new method to predict the failure surface of bedding rock slopes [13]. In addition, some scholars have studied the failure mode of bedding rock landslides, such as Mao et al. found that the weak structure in bedding rock slopes plays a major role in the stability of the slope [14]. Ma and Zhang proposed a sliding-bending-pull-crack composite bedding landslide model and pointed out that the deformation of the slope is obviously controlled by the weak interlayer and rock mass structure [15]. Wang et al. proposed the sliding failure-shear layer failure-tensile failure mode of weak sandwich slopes [16]. Lai et al. proposed a sliding and traction-slip composite sliding mode [17].At present, most of the relevant research studies focus on the stability, deformation and failure modes, dynamic tests, shear resistance of structural planes, and slope reinforcement of bedding slopes. The mechanical response mechanism and the stability algorithm of bedding slope excavation need further study. In order to explore the mechanical response and the stability algorithm during the excavation process of the soft and hard rock interbedded slope, the numerical simulation of the actual project is carried out in this study to study the stress, strain, plastic zone, displacement, and other mechanics of the excavation under heavy rain conditions. According to the characteristics of the sliding surface of the layered slope, the slope stability algorithm is established, and the algorithm is verified by an example.


## 2. Excavation Mechanical Response Analysis

### 2.1. Engineering Overview and Numerical Model

The deep excavation cutting slope of a highway in Guizhou, China, is a typical soft and hard rock interbedded slope. The slope has six soft and hard rock layers with an inclination angle of approximately 15°. The excavation is divided into seven levels, and the maximum excavation depth is 53.6 m. The geological conditions of the excavation area are complex, and the slope sliding trend along the bed is obvious. This area belongs to the tectonic erosion of the landform of low mountains and hills, and the cut slope to be excavated is located on the mountain slope. The local rainfall is abundant, with a total rainfall of 600–700 mm in summer, and floods often occur. The strata from top to bottom in the area discovered by engineering geological survey and drilling are as follows:① Silty clay (Q4dl+el): it belongs to grade II ordinary soil, and the maximum thickness revealed by drilling is 5.93 m② Strongly weathered sandstone (J1-2zl): it belongs to grade IV soft stone, and the maximum thickness revealed by drilling is 13.32 m③ Moderately weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 7.98 m④ Strongly weathered mudstone (J1-2zl): it belongs to III hard soil, and the maximum thickness exposed by drilling is 8.6 m⑤ Medium-weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 6.65 m⑥ Strongly weathered mudstone (J1-2zl): it belongs to grade III hard soil, and the maximum thickness exposed by drilling is 8.45 m⑦ Medium-weathered sandstone (J1-2zl): it belongs to grade V subhard stone, and the maximum thickness is 9.07 m revealed by drilling (not exposed)The main unfavorable geology in the area is bedding slope containing strongly weathered rock stratum, which is prone to bedding sliding under excavation disturbance. This area belongs to the area below seismic intensity VI, the peak value of ground motion acceleration is 0.05 g, and the general characteristic period of ground motion response is 0.35 s. Based on the previous background analysis, it can be seen that the stability of the slope during the excavation process is easily affected by the heavy rain weather and the occurrence and lithology of the rock formation and is less affected by the earthquake factors, so the heavy rain conditions are mainly considered in the numerical simulation in this study. According to the field geological exploration data and experimental data, the selected geotechnical simulation parameters under heavy rain conditions are shown in Table1.Table 1
Geotechnical parameters under heavy rain conditions.
Materialsγ (kN/m3)E (MPa)μc (kPa)φSilty clay189.00E + 010.351513.3°Strongly weathered sandstone24.16.70E + 030.307818.4°Moderately weathered sandstone25.51.50E + 040.2214527.9°Strongly weathered mudstone232.80E + 030.353515.9°According to the actual size of the slope, a 150 m × 70 m numerical model is established, as shown in Figure1, and the CPE4 element type is used for mesh division. Abaqus is used for numerical simulation, the geotechnical parameters under heavy rain conditions in Table 1 are input into the model, the excavation data of slope under heavy rain conditions are obtained through numerical simulation, and the mechanical response analysis is carried out on this basis.Figure 1
Slope model.


### 2.2. Stress Response Analysis

It can be seen from Figure2 of the distribution of shear stress after excavation of slopes at all levels that the shear stress is always in a state of constant adjustment, and its distribution is closely related to the position of each rock layer and the strength of the rock mass and is mostly concentrated in the relatively high-strength moderately weathered sandstone layers. Therefore, the shear stress has obvious band distribution characteristics along the rock layer; in the later stage of excavation, the concentration of shear stress near the toe of the slope on both sides is gradually serious.Figure 2
Nephogram of shear stress distribution: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)In Abaqus,Sxy and Syx are the shear stresses in the XY plane pointing to the Y-axis and X-axis directions, respectively, with the Y-axis direction as a positive unit, kPa. Figure 3 shows the maximum shear stress in the slope body during excavation of slopes at all levels. Comparing the Sxy and Syx curves, it is found that the average shear stress Sxy is about 1.85 times that of Syx during excavation of slopes at all levels, and the former is 95 kPa larger than the latter on average. During the excavation process, the sum of the absolute values of Sxy and Syx was in an increasing trend before the first-level excavation, with a cumulative increase of 51.31%. After the second-level excavation, a single increase of 31.53% was achieved, and it reached the maximum value of 463 kPa for excavation at all levels. However, after the first-level excavation, the shear stress dropped sharply by 22.03%, which indicated that the shear failure of the hard rock layer in the slope body occurred at this time, and most of the shear stress was suddenly released.Figure 3
Variation curve of maximum shear stress of slope excavation at all levels.


### 2.3. Strain and Plastic Zone Analysis

It can be seen from Figure4 of the maximum strain cloud after excavation of the slopes at all levels that the strain during the excavation process is mainly distributed in the soft rock layer with low strength, forming a “strain zone” distributed along the soft rock layer, and the maximum value appears in the soft rock layer. The two adjacent “strain zones” have a tendency to connect up and down, and the two “strain zones” are completely integrated into one after the first-level excavation. The “strain zone” on the side adjacent to the slope surface gradually extended obliquely downward along the soft rock layer and extended to the slope surface after the fourth to the first stage of excavation.Figure 4
Maximum strain distribution nephogram: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation, (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)At the same time, the distribution of plastic zone obtained by numerical simulation is basically the same as the distribution of strain; after excavation of slopes at all levels, the plastic zone and the maximum strain zone are mainly distributed along the soft rock layer. As shown in Figure5, the maximum strain value shows a U-shaped change trend; after the excavation of the first three levels and the last four levels, it first decreased and then increased; especially, after the first level of excavation, the single increase rate was as high as 216.20%. The maximum strain value of the first three levels slope gradually decreased because there was a silty clay layer on the upper part of the slope, and the strain was large in the early stage of excavation, and then, the strain gradually decreased with the enhancement of unloading. The maximum strain value of the last four-level slope increases step by step because the strain gradually accumulates in the local soft rock layer under the combined action of structural change, unloading, and stress redistribution. After the first level of excavation, the toe of the slope is just located in the strongly weathered mudstone layer; under the combined action of the weight of the slope body, the free surface, and the soft rock layer, the strain rapidly accumulates and the maximum strain value increases by 216.20%, indicating that the plastic failure is developing rapidly, and the slope has serious stability problems.Figure 5
Maximum strain curve.


### 2.4. Displacement Response Analysis

It can be seen from the cloud (Figure6) of the horizontal and vertical displacement distribution that the horizontal displacement is mainly manifested as the movement of the layered slope to the outside of the slope along the free side, and the horizontal displacement of the intersection area between the slope toe and the soft rock layer is the most serious with the progress of excavation and forms a layered sliding area with the toe of the slope as the front end of the sliding.Figure 6
(a) Horizontal displacement and (b) vertical displacement cloud chart.
(a)(b)Under the strong unloading action of the slope body, the vertical displacement is divided into the downward settlement of the slope top and the upward uplift of the slope bottom. The influence area of vertical displacement is very wide, and almost, the whole slope moves up or down in different degrees.With the excavation at all levels, the subsidence area is gradually concentrated towards the top of the slope near the mountain body, and the uplift area is gradually concentrated towards the subgrade at the foot of the slope. The distribution of horizontal displacement along the rock layer is obvious and gradually decreases from the toe of the slope to the top of the slope, while the vertical displacement is mainly related to the depth of the rock mass and the progress of excavation.U1 and U2, respectively, represent the maximum values of horizontal displacement and vertical displacement in the model in mm. It can be seen from the curve of the maximum displacement value after excavation of each level of slope shown in Figure 7 that the maximum displacement values in the horizontal and vertical directions both decrease step by step with the excavation before the fifth level and increase with the excavation after the fifth level. The excavation increased step by step; especially, after the first-level excavation, the increase of U1 and U2 reached 130.34% and 45.41%, respectively; the former was about 2.87 times the latter, indicating that, after the first-level excavation, the amount of displacement in the horizontal direction increases rapidly.Figure 7
Variation curve of maximum displacement value.Comparing theU1 and U2 curves in Figure 7, it can be seen that the maximum displacement curves in the two directions are almost parallel, indicating that the two displacement changes are not isolated and have inherent coordination. The increase of the horizontal displacement is about 6.30 times that of the vertical displacement, but the vertical displacement is 5.29 mm larger than the horizontal displacement on average, indicating that the increase of the displacement in the slope is mainly the horizontal displacement, and the displacement value is mainly the vertical displacement. The ratio of U2/U1 is parabolic, and U2/U1 gradually decreases with the fifth stage of excavation, from a maximum of 5.99 to 1.85, indicating that the displacement values in the two directions are accelerating and approaching, and the combined displacement of the two must also be acceleration increases, predicting that the slope is about to slip.In summary, the displacement increase of the layered slope is mainly horizontal displacement, the horizontal displacement has obvious distribution characteristics along the rock layer, and the displacement value is mainly vertical displacement; the change curves of vertical displacement and horizontal displacement are almost parallel. When the ratio ofU2/U1 is less than 1.85, the risk of slippage of the slope is higher.


### 2.5. Sliding Surface Feature Analysis

Since the disturbance and bedding sliding trend caused by the excavation of the first level to fourth level slopes to the soft rock is obviously larger than that of the other levels, therefore, the sliding surface characteristics of the first level to fourth level slopes are mainly analyzed.When the strength of slopes at all levels is reduced to the limit state, the sliding surface has obvious arc-shaped sliding and bedding sliding. When the slope toe is above the soft rock layer, the sliding surface will form an arc-shaped sliding surface between the bottom layer of the soft rock layer and the slope toe, as shown in Figures8(a)–8(d). When the slope toe passes through or under the soft rock layer, the slope body will slide out directly along the bottom layer of the soft rock layer, as shown in Figures 8(b) and 8(c).Figure 8
Schematic diagram of sliding surface of slopes at all levels: (a) fourth stage excavation, (b) third stage excavation, (c) second stage excavation, and (d) first stage excavation.
(a)(b)(c)(d)It can be seen from this that the sliding surface of the bedding slope between soft and hard rocks is complex. According to the relative position of the soft rock layer and the toe of the slope, the sliding surface of soft and hard rock interbedded slope can be divided into two types: arc-shaped—bedding sliding surface and arc-shaped—bedding—arc-shaped sliding surface. If the sliding surface is assumed to be curved surface, it is not consistent with the actual situation.


### 2.6. Mechanical Response Analysis of Strength Reduction Process

It can be seen from Figure9 that, in the process of slope strength reduction at all levels, the responses of stress, strain, and displacement show a certain regularity.Figure 9
Mechanical response curve of slope during strength reduction: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation, (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)In terms of maximum shear stress,Sxy is generally larger than Syx, and Sxy first rises and then falls. It can be seen that the maximum shear stress Sxy has an obvious process of first accumulation and then release, and the shear stress is continuously transferred from the soft rock layer to the hard layer. When the shear stress in the hard rock layer accumulates to its shear strength, shear failure occurs, resulting in the release of shear stress.The change curves of the maximum strainE and the maximum displacements U1 and U2 are basically the same; in the early stage of reduction, they remain basically unchanged or slightly increase, and they suddenly increase rapidly when approaching the final value of reduction. The sudden accelerated development of strain and displacement indicates that a large mechanical change has occurred inside the slope body. Combined with the deformation cloud map of the numerical software, it can be seen that the slope is in a state of imminent sliding at this time. At the same time, comparing the sudden release of the maximum shear stress with the reduction time points of the maximum strain and displacement mutation, it can be seen that the former generally occurs earlier than the latter, so the sudden release of the maximum shear stress does not indicate that the slope is in the limit stable state, while the accelerated changes in maximum strain and displacement reveal that the slope is about to fail.


## 2.1. Engineering Overview and Numerical Model

The deep excavation cutting slope of a highway in Guizhou, China, is a typical soft and hard rock interbedded slope. The slope has six soft and hard rock layers with an inclination angle of approximately 15°. The excavation is divided into seven levels, and the maximum excavation depth is 53.6 m. The geological conditions of the excavation area are complex, and the slope sliding trend along the bed is obvious. This area belongs to the tectonic erosion of the landform of low mountains and hills, and the cut slope to be excavated is located on the mountain slope. The local rainfall is abundant, with a total rainfall of 600–700 mm in summer, and floods often occur. The strata from top to bottom in the area discovered by engineering geological survey and drilling are as follows:① Silty clay (Q4dl+el): it belongs to grade II ordinary soil, and the maximum thickness revealed by drilling is 5.93 m② Strongly weathered sandstone (J1-2zl): it belongs to grade IV soft stone, and the maximum thickness revealed by drilling is 13.32 m③ Moderately weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 7.98 m④ Strongly weathered mudstone (J1-2zl): it belongs to III hard soil, and the maximum thickness exposed by drilling is 8.6 m⑤ Medium-weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 6.65 m⑥ Strongly weathered mudstone (J1-2zl): it belongs to grade III hard soil, and the maximum thickness exposed by drilling is 8.45 m⑦ Medium-weathered sandstone (J1-2zl): it belongs to grade V subhard stone, and the maximum thickness is 9.07 m revealed by drilling (not exposed)The main unfavorable geology in the area is bedding slope containing strongly weathered rock stratum, which is prone to bedding sliding under excavation disturbance. This area belongs to the area below seismic intensity VI, the peak value of ground motion acceleration is 0.05 g, and the general characteristic period of ground motion response is 0.35 s. Based on the previous background analysis, it can be seen that the stability of the slope during the excavation process is easily affected by the heavy rain weather and the occurrence and lithology of the rock formation and is less affected by the earthquake factors, so the heavy rain conditions are mainly considered in the numerical simulation in this study. According to the field geological exploration data and experimental data, the selected geotechnical simulation parameters under heavy rain conditions are shown in Table1.Table 1
Geotechnical parameters under heavy rain conditions.
Materialsγ (kN/m3)E (MPa)μc (kPa)φSilty clay189.00E + 010.351513.3°Strongly weathered sandstone24.16.70E + 030.307818.4°Moderately weathered sandstone25.51.50E + 040.2214527.9°Strongly weathered mudstone232.80E + 030.353515.9°According to the actual size of the slope, a 150 m × 70 m numerical model is established, as shown in Figure1, and the CPE4 element type is used for mesh division. Abaqus is used for numerical simulation, the geotechnical parameters under heavy rain conditions in Table 1 are input into the model, the excavation data of slope under heavy rain conditions are obtained through numerical simulation, and the mechanical response analysis is carried out on this basis.Figure 1
Slope model.


## 2.2. Stress Response Analysis

It can be seen from Figure2 of the distribution of shear stress after excavation of slopes at all levels that the shear stress is always in a state of constant adjustment, and its distribution is closely related to the position of each rock layer and the strength of the rock mass and is mostly concentrated in the relatively high-strength moderately weathered sandstone layers. Therefore, the shear stress has obvious band distribution characteristics along the rock layer; in the later stage of excavation, the concentration of shear stress near the toe of the slope on both sides is gradually serious.Figure 2
Nephogram of shear stress distribution: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)In Abaqus,Sxy and Syx are the shear stresses in the XY plane pointing to the Y-axis and X-axis directions, respectively, with the Y-axis direction as a positive unit, kPa. Figure 3 shows the maximum shear stress in the slope body during excavation of slopes at all levels. Comparing the Sxy and Syx curves, it is found that the average shear stress Sxy is about 1.85 times that of Syx during excavation of slopes at all levels, and the former is 95 kPa larger than the latter on average. During the excavation process, the sum of the absolute values of Sxy and Syx was in an increasing trend before the first-level excavation, with a cumulative increase of 51.31%. After the second-level excavation, a single increase of 31.53% was achieved, and it reached the maximum value of 463 kPa for excavation at all levels. However, after the first-level excavation, the shear stress dropped sharply by 22.03%, which indicated that the shear failure of the hard rock layer in the slope body occurred at this time, and most of the shear stress was suddenly released.Figure 3
Variation curve of maximum shear stress of slope excavation at all levels.


## 2.3. Strain and Plastic Zone Analysis

It can be seen from Figure4 of the maximum strain cloud after excavation of the slopes at all levels that the strain during the excavation process is mainly distributed in the soft rock layer with low strength, forming a “strain zone” distributed along the soft rock layer, and the maximum value appears in the soft rock layer. The two adjacent “strain zones” have a tendency to connect up and down, and the two “strain zones” are completely integrated into one after the first-level excavation. The “strain zone” on the side adjacent to the slope surface gradually extended obliquely downward along the soft rock layer and extended to the slope surface after the fourth to the first stage of excavation.Figure 4
Maximum strain distribution nephogram: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation, (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)At the same time, the distribution of plastic zone obtained by numerical simulation is basically the same as the distribution of strain; after excavation of slopes at all levels, the plastic zone and the maximum strain zone are mainly distributed along the soft rock layer. As shown in Figure5, the maximum strain value shows a U-shaped change trend; after the excavation of the first three levels and the last four levels, it first decreased and then increased; especially, after the first level of excavation, the single increase rate was as high as 216.20%. The maximum strain value of the first three levels slope gradually decreased because there was a silty clay layer on the upper part of the slope, and the strain was large in the early stage of excavation, and then, the strain gradually decreased with the enhancement of unloading. The maximum strain value of the last four-level slope increases step by step because the strain gradually accumulates in the local soft rock layer under the combined action of structural change, unloading, and stress redistribution. After the first level of excavation, the toe of the slope is just located in the strongly weathered mudstone layer; under the combined action of the weight of the slope body, the free surface, and the soft rock layer, the strain rapidly accumulates and the maximum strain value increases by 216.20%, indicating that the plastic failure is developing rapidly, and the slope has serious stability problems.Figure 5
Maximum strain curve.


## 2.4. Displacement Response Analysis

It can be seen from the cloud (Figure6) of the horizontal and vertical displacement distribution that the horizontal displacement is mainly manifested as the movement of the layered slope to the outside of the slope along the free side, and the horizontal displacement of the intersection area between the slope toe and the soft rock layer is the most serious with the progress of excavation and forms a layered sliding area with the toe of the slope as the front end of the sliding.Figure 6
(a) Horizontal displacement and (b) vertical displacement cloud chart.
(a)(b)Under the strong unloading action of the slope body, the vertical displacement is divided into the downward settlement of the slope top and the upward uplift of the slope bottom. The influence area of vertical displacement is very wide, and almost, the whole slope moves up or down in different degrees.With the excavation at all levels, the subsidence area is gradually concentrated towards the top of the slope near the mountain body, and the uplift area is gradually concentrated towards the subgrade at the foot of the slope. The distribution of horizontal displacement along the rock layer is obvious and gradually decreases from the toe of the slope to the top of the slope, while the vertical displacement is mainly related to the depth of the rock mass and the progress of excavation.U1 and U2, respectively, represent the maximum values of horizontal displacement and vertical displacement in the model in mm. It can be seen from the curve of the maximum displacement value after excavation of each level of slope shown in Figure 7 that the maximum displacement values in the horizontal and vertical directions both decrease step by step with the excavation before the fifth level and increase with the excavation after the fifth level. The excavation increased step by step; especially, after the first-level excavation, the increase of U1 and U2 reached 130.34% and 45.41%, respectively; the former was about 2.87 times the latter, indicating that, after the first-level excavation, the amount of displacement in the horizontal direction increases rapidly.Figure 7
Variation curve of maximum displacement value.Comparing theU1 and U2 curves in Figure 7, it can be seen that the maximum displacement curves in the two directions are almost parallel, indicating that the two displacement changes are not isolated and have inherent coordination. The increase of the horizontal displacement is about 6.30 times that of the vertical displacement, but the vertical displacement is 5.29 mm larger than the horizontal displacement on average, indicating that the increase of the displacement in the slope is mainly the horizontal displacement, and the displacement value is mainly the vertical displacement. The ratio of U2/U1 is parabolic, and U2/U1 gradually decreases with the fifth stage of excavation, from a maximum of 5.99 to 1.85, indicating that the displacement values in the two directions are accelerating and approaching, and the combined displacement of the two must also be acceleration increases, predicting that the slope is about to slip.In summary, the displacement increase of the layered slope is mainly horizontal displacement, the horizontal displacement has obvious distribution characteristics along the rock layer, and the displacement value is mainly vertical displacement; the change curves of vertical displacement and horizontal displacement are almost parallel. When the ratio ofU2/U1 is less than 1.85, the risk of slippage of the slope is higher.


## 2.5. Sliding Surface Feature Analysis

Since the disturbance and bedding sliding trend caused by the excavation of the first level to fourth level slopes to the soft rock is obviously larger than that of the other levels, therefore, the sliding surface characteristics of the first level to fourth level slopes are mainly analyzed.When the strength of slopes at all levels is reduced to the limit state, the sliding surface has obvious arc-shaped sliding and bedding sliding. When the slope toe is above the soft rock layer, the sliding surface will form an arc-shaped sliding surface between the bottom layer of the soft rock layer and the slope toe, as shown in Figures8(a)–8(d). When the slope toe passes through or under the soft rock layer, the slope body will slide out directly along the bottom layer of the soft rock layer, as shown in Figures 8(b) and 8(c).Figure 8
Schematic diagram of sliding surface of slopes at all levels: (a) fourth stage excavation, (b) third stage excavation, (c) second stage excavation, and (d) first stage excavation.
(a)(b)(c)(d)It can be seen from this that the sliding surface of the bedding slope between soft and hard rocks is complex. According to the relative position of the soft rock layer and the toe of the slope, the sliding surface of soft and hard rock interbedded slope can be divided into two types: arc-shaped—bedding sliding surface and arc-shaped—bedding—arc-shaped sliding surface. If the sliding surface is assumed to be curved surface, it is not consistent with the actual situation.


## 2.6. Mechanical Response Analysis of Strength Reduction Process

It can be seen from Figure9 that, in the process of slope strength reduction at all levels, the responses of stress, strain, and displacement show a certain regularity.Figure 9
Mechanical response curve of slope during strength reduction: (a) before excavation, (b) seventh stage excavation, (c) sixth stage excavation, (d) fifth stage excavation, (e) fourth stage excavation, (f) third stage excavation, (g) second stage excavation, and (h) first stage excavation.
(a)(b)(c)(d)(e)(f)(g)(h)In terms of maximum shear stress,Sxy is generally larger than Syx, and Sxy first rises and then falls. It can be seen that the maximum shear stress Sxy has an obvious process of first accumulation and then release, and the shear stress is continuously transferred from the soft rock layer to the hard layer. When the shear stress in the hard rock layer accumulates to its shear strength, shear failure occurs, resulting in the release of shear stress.The change curves of the maximum strainE and the maximum displacements U1 and U2 are basically the same; in the early stage of reduction, they remain basically unchanged or slightly increase, and they suddenly increase rapidly when approaching the final value of reduction. The sudden accelerated development of strain and displacement indicates that a large mechanical change has occurred inside the slope body. Combined with the deformation cloud map of the numerical software, it can be seen that the slope is in a state of imminent sliding at this time. At the same time, comparing the sudden release of the maximum shear stress with the reduction time points of the maximum strain and displacement mutation, it can be seen that the former generally occurs earlier than the latter, so the sudden release of the maximum shear stress does not indicate that the slope is in the limit stable state, while the accelerated changes in maximum strain and displacement reveal that the slope is about to fail.


## 3. Research on Stability Algorithms

### 3.1. Slope Stability Algorithm for Conventional Arc-Shaped Sliding Surface

When the slope only bears gravity, its stress field can be obtained by using the elastic solution of elastic wedge. If there is no special geological structure in the slope, its sliding surface is generally arc-shaped.As shown in Figure10, if there is a sliding surface SAB passing through point A at the toe of the slope with a height of h, we take the coordinate system xoy, and the projection of SAB on the plane coordinate system xoy is LAB. Let the expression of LAB be y = f(x), and the rock and soil parameters of the slope are c, φ, and γ, respectively. The physical force component X = 0, Y = γ, the angle between the slope surface and the x-axis is θ, the coordinates of the toe of the slope are A(g, h), and the coordinates of the intersection point between the sliding surface and the top of the slope are B(xb, 0), where g, h, and θ are known quantities.Figure 10
Solution diagram of slope stress field.Referring to the method of Gu and Han [18], it is assumed that the slope stress function ψ is a pure cubic expression of x and y, namely,(1)ψ=ax3+bx2y+cxy2+dy3.From the elastic mechanics, it can be known that the expression of the stress component at a point in the slope is(2)σx=∂2ψ∂y2−X⋅xσy=∂2ψ∂x2−Y⋅yτxy=−∂2ψ∂x∂yX=0,Y=γ.We substitute equation (1) into equation (2) to obtain(3)σx=2cx+6dyσy=6ax+2by−γyτxy=−2bx−2cy.The stress components satisfy the following boundary conditions on the boundary(4)lσx+mτxy=X¯mσy+lτxy=Y¯.Among them, there is no stress on the horizontal boundary (y = 0), i.e., X¯=Y¯=0, and the direction cosine of its outer normal is l=cosπ/2=0m=cosπ=−1; at the same time, there is no stress on the slope boundary y=−xtanθ, i.e., X¯=Y¯=0, and the direction cosine of its outer normal is l=cosπ/2+θ=−sinθm=cosπ−θ=−cosθ, so, it can be obtained that(5)a=b=0c=−γ2tanθd=−γ3tan2θ.We substitute (5) into equation (3) to obtain(6)σx=−1tanθγx−2tan2θγyσy=−γyτxy=1tanθγy,where σx, σy, and τxy are the normal stress and shear stress at a point in the slope, respectively.Regarding the determination of the sliding surface, different methods have different assumptions, and the obtained sliding surface shapes are also different. For example, the simplified Bishop method assumes the sliding surface to be an arc surface. The mechanical analytical rule assumes that the expression of each sliding surface is a quadratic polynomial. The critical sliding surface obtained by the finite element method by connecting the failure points with local stability coefficient less than 1 is also close to the arc surface.As shown in Figure11, there is no layered structure between soft and hard rock in the slope. At this time, the curve segment LAB can be completely regarded as an arc segment, and the normal stress and tangential shear stress at any point on LAB can be expressed as(7)σn=l2σx+m2σy+2lmτxyτn=lmσy−σx+l2−m2τxy,where the expressions of l and m are(8)l=−dyds=−11+dy/dx2dydxm=dxds=11+dy/dx2.Figure 11
Sectional calculation diagram of stability coefficient of sandwich slope: (a) bedding sliding type I and (b) bedding sliding type II.
(a)(b)(14)Fsx=∑i=0n∫xixi+1ci−γitanϕi/1+f′x2f′x2x/tanθ+2fx/tan2θ+1+2f′x/tanθfx1+f′x2dx∑i=0n∫xixi+1γi/1+f′x2f′xfx−x/tanθ−2fx/tan2θ+f′x2−1fx/tanθ1+f′x2dx,On the sliding plane,x and y satisfy y=fx, substitute equation (8) and y=fx into equation (7), and obtain(9)σn=11+dy/dx2dydx2σx+11+dy/dx2σy−21+dy/dx2dydxτxy=11+f′x2f′x2σx+σy−2f′xτxyτn=11+dy/dx2dydxσy−σx+11+dy/dx2dydx2−1τxy=11+f′x2f′xσx−σy+f′x2−1τxy.The Mohr−Coulomb constitutive model is used for slope rock mass, and the stability coefficient of slope mass can be expressed in the following equation:(10)Fs=∫ABc+σntanϕds∫ABτnds,where c and ϕ are the cohesion force and internal friction angle of the rock and soil, respectively, and the numerator and denominator represent the total antisliding force and total sliding force of the slope on the sliding surface SAB under the limit equilibrium state, respectively.We substitute equation (9) into equations (10) and (11) to obtain the following:(11)Fsx=∫abc+1/1+f′x2f′x2σx+σy−2f′xτxytanϕ1+f′x2dx∫ab1/1+f′x2f′xσx−σy+f′x2−1τxy1+f′x2dx.We substitute equation (6) into (11) to obtain(12)Fsx=∫abc−γtanϕ/1+f′x2f′x2x/tanθ+2y/tan2θ+1+2f′x/tanθy1+f′x2dx∫abγ/1+f′x2f′xy−x/tanθ−2y/tan2θ+f′x2−1y/tanθ1+f′x2dx.We substitutey=fx into (12) to obtain(13)Fsx=∫abc−γtanϕ/1+f′x2f′x2x/tanθ+2fx/tan2θ+1+2f′x/tanθfx1+f′x2dx∫abγ/1+f′x2f′xfx−x/tanθ−2fx/tan2θ+f′x2−1fx/tanθ1+f′x2dx.Equation (13) is the calculation formula of stability coefficient of slope with conventional arc-shaped sliding surface, where c, ϕ, and θ are known, and y=fx is obtained. Then, the corresponding stability coefficient can be obtained.


### 3.2. Stability Algorithm of Soft and Hard Rock Interbedded Slope

Based on the analysis of sliding surface in practical engineering, it can be known that the sliding surface of soft and hard rock interbedded slope cannot simply be assumed as a curved surface. Therefore, combined with the engineering practice, the sliding surface of this kind of slope is assumed to be two combination types of arc-shaped—bedding sliding surface (I) and arc-shaped—bedding—arc-shaped sliding surface (II), as shown in Figure11. The coordinates of the intersection point between the sliding plane and the rock surface or slope surface are set as xi,yi. For the convenience of calculation, the coordinates of point A are set as x0,y0, and the coordinates of the boundary points between the bedding plane and the arc plane are Cx1,y1 and Dx2,y2.The soft and hard rock interbedded slope is heterogeneous, and the parametersc, ϕ, and γ in adjacent rock strata differ greatly. Therefore, equation (13) needs to be solved and accumulated according to lithology stratification. Therefore, after the transformation of equation (13), the sliding surfaces of the previous two combination types (I, II) can be expressed bywhere n is the number of soft and hard rock strata in the sliding zone, n = 2, 3, ... , c0=c1, ϕ0=ϕ1, and γ0=γ1, and f(x) is the expression of composite sliding surface.Equation (14) is the calculation formula of stability coefficient of soft and hard rock interbedded slope, and the key to solve the formula lies in the determination of f(x). If an expression of f(x) can be found so that Fs(x) in equation (14) achieves the minimum, then the minimum value is the slope stability coefficient, and f(x), in this case, is the composite sliding surface expression. If the piecewise function f(x) is calculated in strict accordance with the arc plane or along plane, it will be very tedious. Therefore, f(x) is approximately fitted into a multiple continuous function expression.Taking the sliding surface after excavation of the first stage slope obtained by numerical simulation as an example, the curves of the sliding surface were fitted with second-, third-, fourth-, and fifth-order polynomial, respectively. By seeing Figure12, it is found that the higher the order of fitting, the closer the curve is to the sliding surface projection obtained by numerical simulation. When the order of fitting is 4 or above, a relatively ideal effect can be achieved. However, if the order of fitting is too high, it is not conducive to calculation, so this study adopts the fourth-order polynomial for fitting. At the same time, it can be seen that when f(x) is a quadratic polynomial, the fitting curve differs greatly from the projection of the sliding plane, which also confirms that f(x) cannot be assumed as a quadratic polynomial only.Figure 12
Fitting curve of sliding surface of first grade slope.Based on the previous analysis results,y=fx=ax4+bx3+cx2+dx+e can be assumed. The key of the problem is to obtain a set of (a, b, c, d, e) in f(x) so that the minimum value of Fs(x) can be obtained. In fact, it becomes a constrained optimization problem, for which the penalty function method is introduced to solve it.Referring to the practice of Aizhong and Jia [19], the projection LAB of slope sliding surface SAB is isometrically divided into m segments along the vertical direction, as shown in Figure 13. Assuming that the expression of each segment is a quartic polynomial, the curve segment LAB can be expressed aswhere ai, bi, ci, di, ei (i = 1, 2, …, m) is the coefficient to be determined.Figure 13
Schematic diagram of sliding surface segmentation.(15)y=f1x=a1x4+b1x3+c1x2+d1x+e1,x0≤x≤x1∗,f2x=a2x4+b2x3+c2x2+d2x+e2,x1∗≤x≤x2∗,…………fix=aix4+bix3+cix2+dix+ei,xi−1∗≤x≤xi∗,Byy0=f1x=a1x04+b1x03+c1x02+d1x0+e1,(16)y0−a1x04−b1x03−c1x02−d1x0−e1=0.According to the continuity and smoothness off(x), the following can be obtained:(17)fixi=fi+1xixi=fi+1′xi.According to equations (15) and (17), it can be concluded that(18)aixi4+bixi3+cixi2+dixi+ei=ai+1xi4+bi+1xi3+ci+1xi2+di+1xi+ei+14aixi3+3bixi2+2cixi+di=4ai+1xi3+3bi+1xi2+2ci+1xi+di+1.It can be obtained from the known quantities in coordinates of points A and B that(19)x0=g,y0=hym∗=fmxm∗=0.Combined with equations (16), (18), and (19), we can obtain the following:(20)hix=y0−a1x04−b1x03−c1x02−d1x0−e1=0,ai+1−aixi4+bi+1−bixi3+ci+1−cixi2+di+1−dixi+ei+1−ei=0,4ai+1−aixi3+3bi+1−bixi2+2ci+1−cixi+di+1−di=0,x0=g,y0=h,ym∗=fmxm∗=0.In order to ensure that the sliding surface is always inside the slope body during the optimization process, constraint conditions are added to theith subsection point as follows:(21)gix=xi−−fxitanθ=xi+fxitanθ≥0,where 1<i≤m, xi, and −fxi/tanθ are the abscissa values of the sliding surface projection LAB and slope surface projection LAO on the same segment line, respectively.According to the types of the previous problems and characteristics of penalty functions, the mixed penalty function was selected for optimization, and the stability algorithm based on the composite sliding surface search model was established as follows:(22)Fsxmin,x∈Rns.t.hix=0,i=1,2,…,mgix≥0,i=1,2,…,mHx,μ=Fsx+12μ∑i=1mhi2x−μ∑i=1mlngix,where Fsxmin is the objective function, s.t. is the constraint condition, Hx,μ is the mixed augmented objective function, μ is the penalty parameter, and μ⟶0. The calculation process can be realized by MATLAB programming.


## 3.1. Slope Stability Algorithm for Conventional Arc-Shaped Sliding Surface

When the slope only bears gravity, its stress field can be obtained by using the elastic solution of elastic wedge. If there is no special geological structure in the slope, its sliding surface is generally arc-shaped.As shown in Figure10, if there is a sliding surface SAB passing through point A at the toe of the slope with a height of h, we take the coordinate system xoy, and the projection of SAB on the plane coordinate system xoy is LAB. Let the expression of LAB be y = f(x), and the rock and soil parameters of the slope are c, φ, and γ, respectively. The physical force component X = 0, Y = γ, the angle between the slope surface and the x-axis is θ, the coordinates of the toe of the slope are A(g, h), and the coordinates of the intersection point between the sliding surface and the top of the slope are B(xb, 0), where g, h, and θ are known quantities.Figure 10
Solution diagram of slope stress field.Referring to the method of Gu and Han [18], it is assumed that the slope stress function ψ is a pure cubic expression of x and y, namely,(1)ψ=ax3+bx2y+cxy2+dy3.From the elastic mechanics, it can be known that the expression of the stress component at a point in the slope is(2)σx=∂2ψ∂y2−X⋅xσy=∂2ψ∂x2−Y⋅yτxy=−∂2ψ∂x∂yX=0,Y=γ.We substitute equation (1) into equation (2) to obtain(3)σx=2cx+6dyσy=6ax+2by−γyτxy=−2bx−2cy.The stress components satisfy the following boundary conditions on the boundary(4)lσx+mτxy=X¯mσy+lτxy=Y¯.Among them, there is no stress on the horizontal boundary (y = 0), i.e., X¯=Y¯=0, and the direction cosine of its outer normal is l=cosπ/2=0m=cosπ=−1; at the same time, there is no stress on the slope boundary y=−xtanθ, i.e., X¯=Y¯=0, and the direction cosine of its outer normal is l=cosπ/2+θ=−sinθm=cosπ−θ=−cosθ, so, it can be obtained that(5)a=b=0c=−γ2tanθd=−γ3tan2θ.We substitute (5) into equation (3) to obtain(6)σx=−1tanθγx−2tan2θγyσy=−γyτxy=1tanθγy,where σx, σy, and τxy are the normal stress and shear stress at a point in the slope, respectively.Regarding the determination of the sliding surface, different methods have different assumptions, and the obtained sliding surface shapes are also different. For example, the simplified Bishop method assumes the sliding surface to be an arc surface. The mechanical analytical rule assumes that the expression of each sliding surface is a quadratic polynomial. The critical sliding surface obtained by the finite element method by connecting the failure points with local stability coefficient less than 1 is also close to the arc surface.As shown in Figure11, there is no layered structure between soft and hard rock in the slope. At this time, the curve segment LAB can be completely regarded as an arc segment, and the normal stress and tangential shear stress at any point on LAB can be expressed as(7)σn=l2σx+m2σy+2lmτxyτn=lmσy−σx+l2−m2τxy,where the expressions of l and m are(8)l=−dyds=−11+dy/dx2dydxm=dxds=11+dy/dx2.Figure 11
Sectional calculation diagram of stability coefficient of sandwich slope: (a) bedding sliding type I and (b) bedding sliding type II.
(a)(b)(14)Fsx=∑i=0n∫xixi+1ci−γitanϕi/1+f′x2f′x2x/tanθ+2fx/tan2θ+1+2f′x/tanθfx1+f′x2dx∑i=0n∫xixi+1γi/1+f′x2f′xfx−x/tanθ−2fx/tan2θ+f′x2−1fx/tanθ1+f′x2dx,On the sliding plane,x and y satisfy y=fx, substitute equation (8) and y=fx into equation (7), and obtain(9)σn=11+dy/dx2dydx2σx+11+dy/dx2σy−21+dy/dx2dydxτxy=11+f′x2f′x2σx+σy−2f′xτxyτn=11+dy/dx2dydxσy−σx+11+dy/dx2dydx2−1τxy=11+f′x2f′xσx−σy+f′x2−1τxy.The Mohr−Coulomb constitutive model is used for slope rock mass, and the stability coefficient of slope mass can be expressed in the following equation:(10)Fs=∫ABc+σntanϕds∫ABτnds,where c and ϕ are the cohesion force and internal friction angle of the rock and soil, respectively, and the numerator and denominator represent the total antisliding force and total sliding force of the slope on the sliding surface SAB under the limit equilibrium state, respectively.We substitute equation (9) into equations (10) and (11) to obtain the following:(11)Fsx=∫abc+1/1+f′x2f′x2σx+σy−2f′xτxytanϕ1+f′x2dx∫ab1/1+f′x2f′xσx−σy+f′x2−1τxy1+f′x2dx.We substitute equation (6) into (11) to obtain(12)Fsx=∫abc−γtanϕ/1+f′x2f′x2x/tanθ+2y/tan2θ+1+2f′x/tanθy1+f′x2dx∫abγ/1+f′x2f′xy−x/tanθ−2y/tan2θ+f′x2−1y/tanθ1+f′x2dx.We substitutey=fx into (12) to obtain(13)Fsx=∫abc−γtanϕ/1+f′x2f′x2x/tanθ+2fx/tan2θ+1+2f′x/tanθfx1+f′x2dx∫abγ/1+f′x2f′xfx−x/tanθ−2fx/tan2θ+f′x2−1fx/tanθ1+f′x2dx.Equation (13) is the calculation formula of stability coefficient of slope with conventional arc-shaped sliding surface, where c, ϕ, and θ are known, and y=fx is obtained. Then, the corresponding stability coefficient can be obtained.


## 3.2. Stability Algorithm of Soft and Hard Rock Interbedded Slope

Based on the analysis of sliding surface in practical engineering, it can be known that the sliding surface of soft and hard rock interbedded slope cannot simply be assumed as a curved surface. Therefore, combined with the engineering practice, the sliding surface of this kind of slope is assumed to be two combination types of arc-shaped—bedding sliding surface (I) and arc-shaped—bedding—arc-shaped sliding surface (II), as shown in Figure11. The coordinates of the intersection point between the sliding plane and the rock surface or slope surface are set as xi,yi. For the convenience of calculation, the coordinates of point A are set as x0,y0, and the coordinates of the boundary points between the bedding plane and the arc plane are Cx1,y1 and Dx2,y2.The soft and hard rock interbedded slope is heterogeneous, and the parametersc, ϕ, and γ in adjacent rock strata differ greatly. Therefore, equation (13) needs to be solved and accumulated according to lithology stratification. Therefore, after the transformation of equation (13), the sliding surfaces of the previous two combination types (I, II) can be expressed bywhere n is the number of soft and hard rock strata in the sliding zone, n = 2, 3, ... , c0=c1, ϕ0=ϕ1, and γ0=γ1, and f(x) is the expression of composite sliding surface.Equation (14) is the calculation formula of stability coefficient of soft and hard rock interbedded slope, and the key to solve the formula lies in the determination of f(x). If an expression of f(x) can be found so that Fs(x) in equation (14) achieves the minimum, then the minimum value is the slope stability coefficient, and f(x), in this case, is the composite sliding surface expression. If the piecewise function f(x) is calculated in strict accordance with the arc plane or along plane, it will be very tedious. Therefore, f(x) is approximately fitted into a multiple continuous function expression.Taking the sliding surface after excavation of the first stage slope obtained by numerical simulation as an example, the curves of the sliding surface were fitted with second-, third-, fourth-, and fifth-order polynomial, respectively. By seeing Figure12, it is found that the higher the order of fitting, the closer the curve is to the sliding surface projection obtained by numerical simulation. When the order of fitting is 4 or above, a relatively ideal effect can be achieved. However, if the order of fitting is too high, it is not conducive to calculation, so this study adopts the fourth-order polynomial for fitting. At the same time, it can be seen that when f(x) is a quadratic polynomial, the fitting curve differs greatly from the projection of the sliding plane, which also confirms that f(x) cannot be assumed as a quadratic polynomial only.Figure 12
Fitting curve of sliding surface of first grade slope.Based on the previous analysis results,y=fx=ax4+bx3+cx2+dx+e can be assumed. The key of the problem is to obtain a set of (a, b, c, d, e) in f(x) so that the minimum value of Fs(x) can be obtained. In fact, it becomes a constrained optimization problem, for which the penalty function method is introduced to solve it.Referring to the practice of Aizhong and Jia [19], the projection LAB of slope sliding surface SAB is isometrically divided into m segments along the vertical direction, as shown in Figure 13. Assuming that the expression of each segment is a quartic polynomial, the curve segment LAB can be expressed aswhere ai, bi, ci, di, ei (i = 1, 2, …, m) is the coefficient to be determined.Figure 13
Schematic diagram of sliding surface segmentation.(15)y=f1x=a1x4+b1x3+c1x2+d1x+e1,x0≤x≤x1∗,f2x=a2x4+b2x3+c2x2+d2x+e2,x1∗≤x≤x2∗,…………fix=aix4+bix3+cix2+dix+ei,xi−1∗≤x≤xi∗,Byy0=f1x=a1x04+b1x03+c1x02+d1x0+e1,(16)y0−a1x04−b1x03−c1x02−d1x0−e1=0.According to the continuity and smoothness off(x), the following can be obtained:(17)fixi=fi+1xixi=fi+1′xi.According to equations (15) and (17), it can be concluded that(18)aixi4+bixi3+cixi2+dixi+ei=ai+1xi4+bi+1xi3+ci+1xi2+di+1xi+ei+14aixi3+3bixi2+2cixi+di=4ai+1xi3+3bi+1xi2+2ci+1xi+di+1.It can be obtained from the known quantities in coordinates of points A and B that(19)x0=g,y0=hym∗=fmxm∗=0.Combined with equations (16), (18), and (19), we can obtain the following:(20)hix=y0−a1x04−b1x03−c1x02−d1x0−e1=0,ai+1−aixi4+bi+1−bixi3+ci+1−cixi2+di+1−dixi+ei+1−ei=0,4ai+1−aixi3+3bi+1−bixi2+2ci+1−cixi+di+1−di=0,x0=g,y0=h,ym∗=fmxm∗=0.In order to ensure that the sliding surface is always inside the slope body during the optimization process, constraint conditions are added to theith subsection point as follows:(21)gix=xi−−fxitanθ=xi+fxitanθ≥0,where 1<i≤m, xi, and −fxi/tanθ are the abscissa values of the sliding surface projection LAB and slope surface projection LAO on the same segment line, respectively.According to the types of the previous problems and characteristics of penalty functions, the mixed penalty function was selected for optimization, and the stability algorithm based on the composite sliding surface search model was established as follows:(22)Fsxmin,x∈Rns.t.hix=0,i=1,2,…,mgix≥0,i=1,2,…,mHx,μ=Fsx+12μ∑i=1mhi2x−μ∑i=1mlngix,where Fsxmin is the objective function, s.t. is the constraint condition, Hx,μ is the mixed augmented objective function, μ is the penalty parameter, and μ⟶0. The calculation process can be realized by MATLAB programming.


## 4. Example Verification of Stability Algorithm

### 4.1. Stability Algorithm Calculation Based on the Composite Sliding Surface Search Model

The previous composite sliding surface model algorithm is used to focus on the calculation of the sliding surface of the slope from the first stage to the fourth stage with obvious bedding sliding, and the results are shown in Figure14. The expressions y=fx and Fsxmin of potential sliding surface after slope excavation are obtained by calculation; the results are shown in Table 2.Figure 14
The slope sliding surface calculated by the numerical simulation and model. (a) The third- and fourth-grade sliding surfaces. (b) First- and second-order sliding surfaces.
(a)(b)Table 2
Calculation results of composite sliding surface search model algorithm.
Series of slopefxFsxminFourthf(x) = −4.2 ∗ 10−6x4 − 3.3 ∗ 10−4x3 − 0.02x2 − 0.9x + 201.14Thirdf(x) = −3.4 ∗ 10−6x4 − 4.9 ∗ 10−4x3 − 0.03x2 − x + 17.61.05Secondf(x) = −3.4 ∗ 10−6x4 − 5.1∗10−4x3-0.03x2 − x + 17.41.04Firstf(x) = −2.4 ∗ 10−6x4 − 2.9 ∗ 10−4x3 − 0.015x2 − 0.75x + 33.50.85


### 4.2. Algorithm Verification Based on the Composite Sliding Surface Search Model

Firstly, we calculate the numerical simulation results of slope stability coefficient, and at the same time, we calculate the slope stability coefficient that only considers the sliding surface as a quadratic arc-shaped surface under the conventional algorithm (y=fx is assumed to be a second order polynomial) to test the calculation effect of the slope stability algorithm based on the composite sliding surface search model.In the Abaqus strength reduction process, when the strain does not converge and the calculation is forced to stop, it indicates that the slope is greatly deformed and unstable; at this time, the strength reduction coefficient is the slope stability coefficient.Firstly, the calculated data in Abaqus are processed, and the relationship between strain and reduction coefficient during excavation is shown in Figure15.Figure 15
Strain curve during reduction.The ratioK of slope strain increment and reduction coefficient increment corresponding to time t in the calculation of strength reduction is shown in the following equation:(23)K=ΔEtΔFVt=Et−Et−1FVt−FVt−1,where Et and ΔEt are, respectively, strain variables and strain increment at time t and FVt and ΔFVt are, respectively, reduction coefficient and reduction coefficient increment at time t. In (23), assuming that the maximum value Kmax is obtained when t=t′, the corresponding reduction coefficient FVt′ is the slope stability coefficient, namely,(24)Fst=FVt′=Et′−Et′−1Kmax+FVt′−1,where Fst is slope stability coefficient. Equation (24) is used to calculate the stability coefficient Fst on the curve in Figure 15, as shown in Figure 16.Figure 16
Stability coefficient curve obtained by numerical simulation.In the calculation of slope stability coefficient under the arc-shaped sliding surface algorithm,Fsx is expressed in equation (13). In Figure 14, it is assumed that the expression of each segment is a quadratic polynomial; that is, ai=bi=0 and fix=cix2+dix+ei in equation (15). The slope stability coefficient under the quadratic arc-shaped sliding surface can be obtained by optimizing the mixed penalty function. The stability coefficients of the three methods are shown in Table 3.Table 3
Comparison of slope stability results of different methods.
Excavation seriesNumerical simulationArc-shaped sliding surface algorithmComposite sliding surface algorithmFourth1.18 stable1.25 stable1.14 basically stableThird1.06 basically stable1.14 basically stable1.05 basically stableSecond1.07 basically stable1.13 basically stable1.04 unstableFirst0.89 unstable1.04 unstable0.85 unstableAs can be seen from Table3 and the stability coefficient of slope at all levels, the result of the arc-shaped sliding surface algorithm is the largest, and the maximum deviation from the numerical simulation result is 14.42%. The maximum deviation between the composite sliding surface algorithm and the numerical simulation results is only 4.49%. Therefore, the composite sliding surface algorithm is closer to the numerical simulation results than the arc-shaped sliding surface algorithm.Because the arc-shaped sliding surface algorithm only regards the sliding surface as the secondary arc surface and ignores the destruction of bedding sliding on the slope, the calculation results of slope stability coefficient are on the high side, which will overestimate its stability. However, the composite sliding surface algorithm has a relatively lower grade on slope stability, which is easy to draw people’s attention to slope safety. In summary, the composite sliding surface algorithm objectively reflects the destruction of bedding sliding on the slope and is more reasonable than the arc-shaped sliding surface algorithm in the safety assessment of such slopes.


## 4.1. Stability Algorithm Calculation Based on the Composite Sliding Surface Search Model

The previous composite sliding surface model algorithm is used to focus on the calculation of the sliding surface of the slope from the first stage to the fourth stage with obvious bedding sliding, and the results are shown in Figure14. The expressions y=fx and Fsxmin of potential sliding surface after slope excavation are obtained by calculation; the results are shown in Table 2.Figure 14
The slope sliding surface calculated by the numerical simulation and model. (a) The third- and fourth-grade sliding surfaces. (b) First- and second-order sliding surfaces.
(a)(b)Table 2
Calculation results of composite sliding surface search model algorithm.
Series of slopefxFsxminFourthf(x) = −4.2 ∗ 10−6x4 − 3.3 ∗ 10−4x3 − 0.02x2 − 0.9x + 201.14Thirdf(x) = −3.4 ∗ 10−6x4 − 4.9 ∗ 10−4x3 − 0.03x2 − x + 17.61.05Secondf(x) = −3.4 ∗ 10−6x4 − 5.1∗10−4x3-0.03x2 − x + 17.41.04Firstf(x) = −2.4 ∗ 10−6x4 − 2.9 ∗ 10−4x3 − 0.015x2 − 0.75x + 33.50.85


## 4.2. Algorithm Verification Based on the Composite Sliding Surface Search Model

Firstly, we calculate the numerical simulation results of slope stability coefficient, and at the same time, we calculate the slope stability coefficient that only considers the sliding surface as a quadratic arc-shaped surface under the conventional algorithm (y=fx is assumed to be a second order polynomial) to test the calculation effect of the slope stability algorithm based on the composite sliding surface search model.In the Abaqus strength reduction process, when the strain does not converge and the calculation is forced to stop, it indicates that the slope is greatly deformed and unstable; at this time, the strength reduction coefficient is the slope stability coefficient.Firstly, the calculated data in Abaqus are processed, and the relationship between strain and reduction coefficient during excavation is shown in Figure15.Figure 15
Strain curve during reduction.The ratioK of slope strain increment and reduction coefficient increment corresponding to time t in the calculation of strength reduction is shown in the following equation:(23)K=ΔEtΔFVt=Et−Et−1FVt−FVt−1,where Et and ΔEt are, respectively, strain variables and strain increment at time t and FVt and ΔFVt are, respectively, reduction coefficient and reduction coefficient increment at time t. In (23), assuming that the maximum value Kmax is obtained when t=t′, the corresponding reduction coefficient FVt′ is the slope stability coefficient, namely,(24)Fst=FVt′=Et′−Et′−1Kmax+FVt′−1,where Fst is slope stability coefficient. Equation (24) is used to calculate the stability coefficient Fst on the curve in Figure 15, as shown in Figure 16.Figure 16
Stability coefficient curve obtained by numerical simulation.In the calculation of slope stability coefficient under the arc-shaped sliding surface algorithm,Fsx is expressed in equation (13). In Figure 14, it is assumed that the expression of each segment is a quadratic polynomial; that is, ai=bi=0 and fix=cix2+dix+ei in equation (15). The slope stability coefficient under the quadratic arc-shaped sliding surface can be obtained by optimizing the mixed penalty function. The stability coefficients of the three methods are shown in Table 3.Table 3
Comparison of slope stability results of different methods.
Excavation seriesNumerical simulationArc-shaped sliding surface algorithmComposite sliding surface algorithmFourth1.18 stable1.25 stable1.14 basically stableThird1.06 basically stable1.14 basically stable1.05 basically stableSecond1.07 basically stable1.13 basically stable1.04 unstableFirst0.89 unstable1.04 unstable0.85 unstableAs can be seen from Table3 and the stability coefficient of slope at all levels, the result of the arc-shaped sliding surface algorithm is the largest, and the maximum deviation from the numerical simulation result is 14.42%. The maximum deviation between the composite sliding surface algorithm and the numerical simulation results is only 4.49%. Therefore, the composite sliding surface algorithm is closer to the numerical simulation results than the arc-shaped sliding surface algorithm.Because the arc-shaped sliding surface algorithm only regards the sliding surface as the secondary arc surface and ignores the destruction of bedding sliding on the slope, the calculation results of slope stability coefficient are on the high side, which will overestimate its stability. However, the composite sliding surface algorithm has a relatively lower grade on slope stability, which is easy to draw people’s attention to slope safety. In summary, the composite sliding surface algorithm objectively reflects the destruction of bedding sliding on the slope and is more reasonable than the arc-shaped sliding surface algorithm in the safety assessment of such slopes.


## 5. Conclusion

(1)
The stress, strain, plastic zone, and horizontal displacement of bedding slope have obvious distribution characteristics along the rock surface. The stress, strain, and plastic zone are concentrated in hard rock and soft rock, respectively. The increase of slope displacement is mainly horizontal displacement, and the larger displacement value is mainly vertical displacement.(2)
The strain and displacement of bedding slope will accelerate when the slope is close to the ultimate stability state, and the sudden release of shear stress usually occurs earlier than the ultimate stability state.(3)
According to the relative position of soft rock layer and slope foot, the sliding surface of bedding slope can be divided into two types, namely, arc-shaped—bedding sliding surface and arc-shaped—bedding—arc-shaped sliding surface.(4)
The composite sliding surface algorithm objectively reflects the destruction of the slope by bedding sliding and can reasonably evaluate the safety of soft and hard rock interbedded slope.



---

*Source: 3331097-2022-11-17.xml*