Abstract

The rock slopes with tunnels appear widely in the actual project, but there is no executable basis for the seismic stability calculation of the rock slope with tunnel. According to the upper bound theorem of plastic limit analysis and pseudostatic method, the upper bound solution of the safety factor of the rock slope with tunnel was rigorously derived under earthquake loading. This upper solution takes into account the design parameters of the slope and the tunnel, the horizontal and vertical seismic loads, and the physical and mechanical parameters of the rock mass. Comparing the calculated results with the existing results, the validity of the proposed method was verified. The sensitivity and influence of different parameters on the seismic stability of the slope were analyzed. The results show that the three factors such as the horizontal seismic force coefficient, the slope height, and the internal friction angle are the three key factors that influence the sensitivity of the safety factor and have a great effect on it.

1. Introduction

In 2017, the Jiuzhaigou earthquake with a moment magnitude of 7.0 triggered a large number of slope failures in Sichuan Province, China. The earthquake-induced slope failures led to serious damage to the highway and tunnel system. The stability of rock and soil slope [1, 2] under earthquake has become a hot research topic at present [35]. Based on the quasi-static method, the limit equilibrium analysis and upper bound approach of limit analysis have been widely utilized to evaluate the seismic stability of the slope. Deng et al. [6] applied the limit equilibrium method to analyze the stability of the slope under earthquake considering three kinds of the slip surface (line, circular, and arbitrary curves). Lu et al. [7] extended Newmark’s method to three dimensions and proposed a new method for evaluating the seismic permanent displacement of 3D slope. Liang and Knappett [8] presented an improved Newmark sliding block procedure for predicting the seismic slip of a vegetated slope. Although the limit equilibrium analysis is simple and has certain accuracy, the solution obtained by the limit equilibrium method is not the upper or lower limit of the real solution, which has some limitations in theory [9]. The upper bound method of limit analysis has made great progress in slope engineering [10] since Chen [11] introduced the plastic limit analysis into the stability analysis of the soil slope. Zhao et al. [12] applied the upper bound limit analysis theorem and the shear strength reduction technique to define the safety factor and corresponding critical failure mechanism of a layered soft-rock slope. Liu et al. [13] presented a new approach for determining the factor of safety and the corresponding critical slip surface of a layered rock slope under seismic excitations, which were obtained by the limit equilibrium method and pseudostatic approach. A 3D rotational mechanism was adopted by He et al. [14] to analyze seismic displacement of slopes reinforced with piles using limit analysis theory and Newmark’s analytical procedure. Ausilio et al. [15, 16] used the kinematic approach of limit analysis to analyze the seismic stability of slopes reinforced with geosynthetics and slopes reinforced with piles, respectively.

However, most of the published works of slope limit analysis were mainly concentrated at the research of a slope without tunnel. Few scholars have carried out theoretical research on the slope with tunnel [17]. The rock slope with tunnel, as a composed slope structure, may have a complex interaction between the tunnel and slope. The stability calculation theory under natural or earthquake condition is far from keeping pace with the engineering construction of the slope with tunnel. In the present study, the failure mechanism of the rock slope with tunnel is constructed, and the upper bound solution of the safety factor of the rock slope is deduced based on the upper bound limit analysis theorem and the pseudostatic approach. In addition, the solutions presented in this study are compared with those obtained by the pseudostatic approach. The effects of physical and mechanical parameters [18, 19] of the rock slope and design parameters of the tunnel slope with tunnel are investigated.

2. Upper Bound Theorem of Plastic Limit Analysis

The upper bound theorem of plastic limit analysis is a useful method in the stability limit analysis, which has attracted extensive attention to solving geotechnical problems [2022]. In the stability analysis of the rock and soil slope, it is necessary to understand the damage load of rock and soil when it begins to produce unrestricted plastic flow; it is not necessary to know the change process of stress and strain with the external load. The upper limit theorem of the plastic limit analysis assumes that the rock mass moves in the form of a rigid plastic body during failure. It is required that the rate of work of force on the slope surface and the physical force is not greater than the energy dissipation in the permissible velocity field for arbitrary maneuver admissible damage mechanism, which means that the external rate of work is not greater than the internal energy dissipation rate. The external rate of work generally includes the rate of work induced by the rock mass, the pore water pressure, the slope overload, and the seismic load. The internal energy dissipation rate includes the internal dissipated rate of the failure surface and the rate of work induced by the resistance of the supporting structure.

2.1. Basic Assumptions

In the stability analysis of the system of the tunnel and slope, in order to simplify the research object, the rock and soil mass and the tunnel structure are usually treated separately as two separate structures, which is called a structured method [23]. The pseudostatic method [24] is often used to analyze the dynamic effect of earthquake in engineering practice.

Based on the relevant research results [25, 26], the following basic assumptions are made for the convenience of analysis. (1) The rock mass of the slope is considered to be soft rock which is homogeneous and isotropic. (2) The tunnel surrounding rock is regarded as an ideal elastomeric body following the Mohr–Coulomb yield criterion. (3) The failure surface is a log-spiral failure surface passing through the slope toe and tunnel. (4) The stability analysis of the slope with tunnel is simplified as a plane strain problem. (5) The tunnel is simplified as a circular cross section for calculation. The vertical surrounding rock pressure acting on the tunnel vault is simplified as a linear uniform load q, and the horizontal surrounding rock pressure acting on the side wall is simplified as a linear uniform load e. (6) The shear strength parameters of the rock material do not change with the action of the earthquake, and the quasi-static method is used to analyze the seismic effect.

2.2. Failure Mechanism

Constructing reasonable failure mode is the premise and key of upper bound limit analysis. The strength reduction analysis module of MIDAS GTS/NX was used to simulate the failure mode of the slope with tunnel. The elasticity modulus E of rock mass was 0.6 GPa, the cohesion c was 0.15 MPa, the internal frictional angle φ was 20°, and the rock unit weight γ was 19 kN·m−3. The equivalent plastic strain is a measure to analyze the comprehensive deformation state of the material under the condition of complex stress compared with uniaxial tension. Figure 1 shows the equivalent strain distribution of the slope with tunnel. The sliding surface of the slope passes through the left sidewall and the right spandrel of the tunnel from the slope toe and extends upward to the top of the slope. It is similar to the most unfavorable position of the goaf in the slope, which is proposed by the literature [27]. When the geometric center of the goaf lies on the sliding belt, the stability of the slope is much lower than that of the geometric center of the goaf which lies in the interior of the sliding belt or lies in the outside of the sliding belt. Referring to related literatures [28], the log-spiral failure surface is consistent with the actual failure surface; the failure mode of kinematically permissible velocity field is close to the actual failure mode, so the log-spiral failure mechanism is used in this paper.

In order to make the log spiral as far as possible passing through the geometric center of the tunnel, the log-spiral failure mechanism used in this paper is illustrated in Figure 2, where the failure surface passes through the left sidewall point C and the right spandrel point B of the tunnel. R is the tunnel radius. The horizontal and vertical distances from point O′ to point D are d1 and d2, respectively. The area ABCDAA is regarded as a rigid plastic body, which rotates about the rotation center O with angular velocity Ω relative to the static material below the log-spiral surface AD. Therefore, the surface AD is the velocity discontinuity surface of the failure mechanism studied in this paper. As shown in Figure 2, r0 and θ0 are the length and angle of the baseline OA, respectively; rh and θh are the length and angle of the baseline OD, respectively; rB and θB are the length and angle of the beeline OB, respectively; rC and θC are the length and angle of the beeline OC, respectively. The height of the overall failure mechanism is H. The vertical distance from the top of the slope to the right spandrel point B of the tunnel is α3H. The vertical distance from point B to point C is α2H. The vertical distance from the left sidewall point C to the slope toe point D is α1H, where α1, α2, and α3 are the height coefficients which can be determined by (3). The length of the beeline AA is L. The failure mechanism is determined by θ0, θh, and β. BCT is an intersection area of the tunnel and the sliding surface.

The equation of the log-spiral curve can be expressed aswhere φ is the included angle between the strain velocity vector V at any point on the sliding surface and the slip line at that point, that is, the internal friction angle of the surrounding rock. r and θ are the length and angle of the connecting line between a point on the sliding surface and point O, respectively. The length of the baseline OD can be expressed as

According to the geometric relationship in Figure 2, we obtain

3. Calculation of Energy Consumption

The external force includes the gravity of rock mass and seismic load. The internal force includes the resisting force required by the lining to balance the pressure of the surrounding rock and the adhesion force on the sliding surface. According to the upper bound limit analysis, the rate of the external force in the strain velocity field can be expressed as , where is the rate of the rock mass gravity; and are the rate of the horizontal and vertical seismic inertia force, respectively; kh and kv are the horizontal and vertical seismic force coefficients, respectively. The rate of internal energy dissipation can be expressed as , where is the rate of lining resistance and is the energy dissipation rate on the velocity discontinuity surface.

3.1. Rate of External Force
3.1.1. Rate of Rock Mass Gravity

It is very difficult to calculate the rate of the rock mass gravity directly from the area ABCDAA by integration. Firstly, the rock mass gravity rate of each area is obtained by the indirect method. Then, the rock mass gravity rate of the area ABCDAA is calculated by simple algebraic summation, which can be expressed aswhere , , , and are the rock mass gravity rates done by areas OAD, OAA′, OAD, and BCT, respectively; f1, f2, f3, and f4 are the functions of (θh, θ0) and γ is the unit weight of rock mass. Figure 3 shows the detailed calculation of the function f.

Because the intersection area of the tunnel and the sliding surface is smaller relative to the overall rotation area, the intersection area is simplified as trapezoid to solve the rock mass gravity rate of the area BCCB′. The function f1f4 can be written as

3.1.2. Rate of Seismic Inertia Force

The rate done by seismic inertia force includes the horizontal seismic inertia rate and vertical seismic inertia rate . The rate done by seismic inertia force of each region is calculated first, and then, the superposition calculation is carried out. and can be expressed aswhere f1, f2, f3, and f4 are the same as above and f5, f6, f7, and f8 are the functions of (θh, θ0). The function f5f8 can be written as

The sum of the external force rate can be expressed aswhere

3.2. Rate of Internal Energy Dissipation
3.2.1. Rate of Energy Dissipation on the Velocity Discontinuity Surface

When the plastic failure occurs in the slope, it can be considered that there is no rate dissipation in the rigid plastic area of ABCDAA, and the local deformation is basically the same. The energy dissipation rate on the velocity discontinuity surface can be calculated aswhere c is the cohesion of rock mass and f9 is the function of (θh, θ0), as shown below. Other parameters are the same as above.

3.2.2. Rate of Lining Resistance

Because the log-spiral failure surface passes through the tunnel, the vertical and horizontal surrounding rock pressure above the log-spiral surface is only considered. The force on the lining structure can be seen as the surrounding rock pressure. Because the lining structure prevents the damage of surrounding rock, the rate of lining resistance is equal to the rate of surrounding rock pressure, and the symbol is opposite, as shown in Figure 4. It is assumed that the exerting point of vertical surrounding rock pressure is at the midpoint of the horizontal distance between B′ and C′ and of horizontal surrounding rock pressure is at the midpoint between BB′ and CC′. Hm is the buried depth of the tunnel, and Ht is the tunnel height.

The rates done by vertical and horizontal lining resistance are and , respectively, which can be expressed aswhere Fh is the horizontal surrounding rock pressure; Fvl and Fvr are the vertical rock pressure on the left and right sides of the tunnel, respectively; f10, f11, and f12 are the functions of (θh, θ0), which can be expressed as

The sum of the rate of lining resistance can be expressed as

Thus, the sum of the internal energy dissipation can be expressed as

When the external rate is equal to the internal energy dissipation rate, the upper bound solution of the stability coefficient of the slope with tunnel can be obtained as .

4. Calculation Results and Influencing Factor Analysis

4.1. Comparison of Calculation Results

The slope stability coefficient calculated by this paper is compared with the result of the limit equilibrium method [6, 29] to verify the validity of the proposed method. According to the upper bound theorem of limit analysis, the safety factor of slope stability can be expressed as

Substituting (19) and (30) into (31), (31) becomes

Because Fs is a function of θ0, θh, and β and the right side of the above equation implies a reduction factor Fs, the function has a minimal upper bound when θ0, θh, and β satisfy the following conditions:

In order to avoid lengthy calculation, the semigraphical method [30] is used to solve (33). Substituting obtained θ0, θh, and β into (33), the minimum safety factor of slope stability can be obtained. The sequential quadratic optimized iterative method [31] is used to verify the results.

In addition, the safety factor Fs should meet the following constraint conditions:

Table 1 shows the comparison of the seismic safety factor calculated by various methods. From the calculation results of Table 1, the results obtained in this study are close to the results obtained by the simplified Bishop’s method. The safety factors of the present study are smaller than those calculated by the simplified Bishop’s method. The good agreements show that the proposed method for calculating seismic stability of the slope is an effective method.

4.2. Sensitivity Analysis of Seismic Stability Factor

The main factors affecting the seismic stability of the slope include the design parameters of the slope and the tunnel, physical and mechanical parameters of the rock mass of the slope, and horizontal and vertical seismic loads. In order to simplify the analysis, seven influencing factors of the slope with tunnel are selected for analysis, including the slope height H, the slope angle β, the distance of the slope toe from the tunnel d1, the tunnel radius R, the horizontal seismic force coefficient kh, the cohesion of rock mass c, and the internal friction angle φ. According to the experience of engineering design, each factor takes three levels for orthogonal analysis, which represent high, medium, and low levels of each factor, as shown in Table 2. They are arranged in the L27 orthogonal table, as shown in Table 3. The interaction between factors is not considered. The vertical seismic force coefficient kv is 2/3 times kh. Other factors, such as γ = 20 kN·m−3, φ = 25°, α1 = 0.1, α2 = 0.2, α3 = 0.7, and d2 = 4 m, are used for calculations. Because each factor of each level corresponds to nine groups of tests, the average value of Fs corresponding to these nine groups of tests is taken as the critical value of the safety factor. The sensitivity of the factor to the safety factor of the slope can be obtained, as shown in Table 4.

The greater the range, the greater the sensitivity of the test results to this factor. From Table 4, the sensitivities of the abovementioned seven influencing factors are as follows: kh > H > φ > β > c > R > d1. The safety factor is more sensitive to the horizontal seismic force coefficient kh and the slope height H. The main reason is that the change of the height of the slope greatly changes the whole shape of the slope: as the height of the slope increases, the slope angle will increase and the stability of the slope will be lowered. In this study, the pseudostatic method is used to simplify the seismic action, and then, the pseudostatic force is applied to the research object, so the sensitivities of the horizontal seismic force coefficient and the slope height are larger. The empty space in the tunnel is small relative to the entire sliding body, so the sensitivities of the main design parameters of the tunnel are smaller.

4.3. The Influence of Factors on Safety Factor

The control variable method was used to study the influence of the slope height H on the safety factor. The main parameters of the slope are as follows: β  = 45°, γ  = 20 kN·m−3, c = 70 kPa, φ  = 20°, α1 = 0.1, α2 = 0.2, α3 = 0.7, kh = 0.4 g, kv = 0.267 g, d2 = 4 m, d1 = 9 m, and R = 6 m. Figure 5 shows the change trend of the safety factor with the slope height. The safety factor presents a significant nonlinear decreasing trend with the increase of the slope height.

The slope height is taken as H = 15 m when calculating, and the other parameters are the same as above. Figure 6 shows the change trend of the safety factor with the slope angle. The safety factor decreases with increasing slope angle, and it presents nonlinearity.

The slope height is taken as H = 15 m when calculating, and the other parameters are the same as above. Figure 7 shows the change trend of the safety factor with the distance of the slope toe from the tunnel. The safety factor decreases with increasing the distance of the slope toe from the tunnel, and it presents approximate linearity. The influence of the tunnel on the stability of the slope is related to the position of the tunnel relative to the position of the sliding surface. The increase of the distance of the slope toe from the tunnel may decrease the strength of the sliding belt. The antislide force of the corresponding sliding surface will decrease, and the safety factor of the slope will decrease.

The parameters are the same as above. Figure 8 shows the change trend of the safety factor with the tunnel radius. The safety factor increases with increasing tunnel radius, and it presents nonlinearity. The increase of the tunnel radius may decrease the mass of the sliding zone; the glide force of the corresponding sliding surface will decrease, and the safety factor of the slope will increase.

The parameters are the same as above. Figure 9 shows the change trend of the safety factor with the cohesion. The safety factor increases with increasing cohesion, and it presents linearity.

The parameters are the same as above. Figure 10 shows the change trend of the safety factor with the horizontal seismic force coefficient. The safety factor decreases with increasing horizontal seismic force coefficient, and it presents obvious nonlinearity.

The parameters are the same as above. Figure 11 shows the change trend of the safety factor with the internal friction angle. The safety factor increases with increasing internal friction angle, and it presents obvious nonlinearity.

Some parameters, such as the height of the slope H, the internal friction angle φ, and the horizontal seismic force coefficient kh, have a great effect on the safety factor; this may be because these parameters are related to the overall system of the slope. Others have a little effect on the safety factor, and the reason may be that these factors are only related to the tunnel and sliding surface inside the slope.

5. Conclusions

Based on the upper bound theorem of plastic limit analysis and the pseudostatic method, the seismic stability of the slope with tunnel was analyzed by defining the safety factor. The sensitivity analysis of the influence factors of the safety factor was carried out by the orthogonal analysis method, and the influence degree of these factors was analyzed. The main conclusions are as follows:(1)A logarithmic spiral failure mechanism was established for the slope with tunnel. The upper bound solution of the safety factor of the slope with tunnel was obtained by calculating the external work rate and the internal energy dissipation rate. The upper bound solution takes into account the tunnel location, buried depth of the tunnel, resistance of the lining structure, seismic inertia force coefficient, cohesion and internal friction angle of rock mass, and other factors. The consistency of the results between the present study and the existing research shows that the present solutions are feasible.(2)The sensibility sequence of influence factors of the seismic safety factor is listed as follows: horizontal seismic force coefficient kh, slope height H, internal friction angle φ, slope angle β, cohesion of rock mass c, tunnel radius R, and the distance of the slope toe from the tunnel d1. Among them, the sensitivities of the horizontal seismic force coefficient kh and the slope height H are larger, while the sensitivities of the tunnel radius R and the distance of the slope toe from the tunnel d1 are smaller.(3)The slope angle β, the tunnel radius R, the cohesion c, and the distance of the slope toe from the tunnel d1 have little effect on the safety factor, while the slope height H, the horizontal seismic force coefficient kh, and the internal friction angle φ have a great effect on the safety factor.(4)This paper assumes that the log-spiral failure surface passes through the foot of the slope and passes through the tunnel. Such assumptions are not appropriate in situations where the slope angle is very small or the log-spiral surface does not pass through the tunnel. The influence of the bias tunnel is not considered when calculating the surrounding rock pressure, which affects the calculation results to a certain extent. The seismic stability of the slope with tunnel requires a comprehensive research.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors are grateful for the financial support from the National Natural Science Foundation (NNSF) of China through Grant nos. 51204215 and 51404309.