Modified Bishop method for stability analysis of weakly sloped subgrade under centrifuge model test

Abstract

The sliding forms of weak sloped and horizontal subgrades during the sliding process differ. In addition, the sliding form of weakly sloped subgrades exhibits considerable slippage and asymmetry. The accuracy of traditional slice methods for computing the stability safety factor of weakly sloped subgrades is insufficient for a subgrade design. In this study, a novel modified Bishop method was developed to improve the accuracy of the stability safety factor for different inclination angles. The instability mechanism of the weakly sloped subgrade was considered in the proposed method using the “influential force” and “additional force” concepts. The “additional force” reflected the weight effect of the embankment fill, whereas the “influential force” reflected the effect of the potential energy difference. Numerical simulations and experimental tests were conducted to evaluate the advantages of the proposed modified Bishop method. Compared with the traditional slice method, the error between the proposed method and the exact value is less than 32.3% in calculating the safety factor.

Graphical abstract

Keywords

weakly sloped subgrade / stability analysis / additional force / influential force / modified Bishop method

1 Introduction

Slope stability is a crucial topic of highway and railway construction, particularly in south-west China, where weakly sloped subgrades widely exist. The weakly sloped subgrade has two characteristics: (1) the bottom of the embankment or the soft soil layer is inclined, and (2) the subgrade soil is loose and exhibits a relatively low strength and high compressibility.

The slice method based on the limit equilibrium is widely used to analyze the safety factor of the slope stability. Bishop [1] and other researchers [2–5] have made notable contributions to this field and have developed various analytical methods.

A weakly sloped subgrade is prone to stability problems, such as slope collapse and an excessive lateral deformation, owing to the inclination angle [6–11]. According to the literature, research on the subgrade stability can be divided into theoretical and experimental studies.

Recent developments based on a theoretical analysis can be summarized as follows. With a constant modification and improvement of the limit equilibrium method, the main research directions focus on the sliding surface shape, safety factor, and three-dimensional limit equilibrium method. The method of analyzing the sliding surface is first investigated [12,13]. Next, the constraint of the sliding surface problem is transformed to determine the dynamic upper and lower bounds of the control variables [14]. To improve the accuracy of safety factor calculation, the probability and sensitivity of the stability of weakly sloped composite subgrades have been analyzed [15–17]. Meanwhile, some scholars have extended two-dimensional solutions into three dimensions [18]. In addition, a stability model of the asymmetric slope was established using different traditional slice methods, but only a single-sliding direction was specified for the entire failure mass [19]. Based on Cheng’s stability model, the applicability of a dynamic programming method of slope stability analysis was evaluated [20]. However, none of the methods mentioned above can be directly applied to the stability analysis of weakly sloped subgrades.

A different group of researchers focused on experimental studies. An integrated analysis of the deformation and failure processes was applied to investigate the behavior and mechanism of a slope failure [21–24]. The effects of the slope angle and dynamic force on the slope deformation characteristics have been determined through laboratory model tests [25]. The effects on the gravitational deformation of the resulting slopes were investigated based on centrifuge model tests [26]. The influence of the slope height on the role of vegetation for improving the seismic slope stability was investigated through a centrifuge test [27]. The variation trends of the horizontal displacement, vertical settlement, and stability of weakly sloped subgrades during the filling process have also been compared and analyzed [28]. Numerical simulations are also important research methods in subgrade stability research. Meshfree methods [29] and DEM [30–33] are popular approaches to solving large deformation problems. Peridynamics [34,35] combines the advantages of both meshfree methods and molecular dynamics method. However, its dynamic equation parameters are difficult to obtain. These numerical methods are more suitable for studying the process of subgrade crack development. To obtain higher calculation accuracy, a smooth extended finite element method [36] and a phantom node method [37] can be adopted. Some useful conclusions have been drawn through numerical simulation studies, and the stability of the weakly sloped subgrade was evaluated based on the experimental or numerical simulation results. However, the effects of the potential difference and the lateral embankment of the sliding body were ignored. In this study, two novel concepts, i.e., an “influential force” and an “additional force,” were established through experimental tests and numerical simulations to reflect such effects. The influential force and additional force reflect the effects of the potential difference and the lateral embankment of the sliding body, respectively. The modified Bishop method was established by introducing two forces on the basis of the Bishop method. Furthermore, the modified Bishop method was validated in the real engineering example.

The contributions of this study can be summarized as follows: (1) two novel concepts, i.e., the influential force and additional force, are introduced, (2) the acting position and direction of the two forces are determined, (3) the development of a modified Bishop method is described, and (4) the safety factor calculated using the modified Bishop method is shown to be more accurate.

2 Classification of weakly sloped subgrade

A weakly sloped subgrade is a common unfavorable geological condition that exists in mountainous and hilly areas. As shown in Fig. 1, {\theta }_{\mathrm{1}} is the inclination angle of the embankment bottom, {\theta }_{\mathrm{2}} is the inclination angle of the bottom of the soft soil layer, and h is the thickness of the soft soil layer at the center of the subgrade. Based on different inclinations {\theta }_{\mathrm{1}} and {\theta }_{\mathrm{2}}, weakly sloped subgrades can be divided into three forms.

(1) The subface of the embankment is horizontal, and the subface of the soft soil layer is inclined, i.e., {\theta }_{\mathrm{1}}=\mathrm{0}, and {\theta }_{\mathrm{2}}\ne \mathrm{0}.

(2) The subface of the embankment is inclined, and the subface of the soft soil layer is horizontal, i.e., {\theta }_{\mathrm{1}}\ne \mathrm{0}, and {\theta }_{\mathrm{2}}=\mathrm{0}.

(3) Both the subfaces of the embankment and soft soil layer are inclined, i.e., {\theta }_{\mathrm{1}}\ne \mathrm{0}, and {\theta }_{\mathrm{2}}\ne \mathrm{0}.

To capture the instability process of a weakly sloped subgrade, the centrifuge model test and numerical simulations were conducted to demonstrate the rationality of both the influential force and additional force.

3 Centrifuge model tests

Centrifuge model tests were conducted in this study to investigate the instability of the weakly sloped subgrade prototype. The original size of the under-investigated prototype was 73.3 m × 40 m, and the prototype was scaled down by 67-fold for applying the centrifuge model tests.

A centrifuge model test is an effective approach for analyzing the failure mechanisms of weakly sloped subgrades, although only somewhat limited results have been reported. A centrifuge model test can compensate for the loss of soil weight during the model rescaling process.

3.1 Test device

A series of centrifuge model tests were conducted at Nanjing Hydraulic Research Institute, and the final effective size of the model container (Fig. 2) was 1100 mm × 400 mm × 600 mm (length × width × height) based on the test results. The subgrade model was placed in the model container. To observe the deformation of the subgrade directly, a 100-mm thick transparent organic glass was installed in front of the model container. The model side was pasted with a 40 mm × 40 mm deformed mesh strip (Fig. 3). Through the monitoring system (organic glass), the deformation of the model during the tests could be recorded and the soil displacement vector diagram could be plotted.

3.2 Rationality of centrifuge model test

A dimensional analysis was conducted to verify the rationality of the model test. In the model and the prototype of weakly sloped subgrade, the relationships of dimensionless size are as follows:

(1){\theta }_{\mathrm{2}}{}_{\mathrm{m}}={\theta }_{\mathrm{2n}},{\mu }_{\mathrm{m}}={\mu }_{\mathrm{n}},

where “m” represents the model and “n” represents the prototype.

The relationships of the dimensional size are defined as follows:

(2)h_{\mathrm{m}}=N{h}_{\mathrm{n}},a_{\mathrm{m}}=N{a}_{\mathrm{n}}=Ng,

where h_{\mathrm{m}} is the thickness of the soft soil layer of the model, h_{\mathrm{n}} is the thickness of the soft soil layer of the prototype, a_{\mathrm{m}} is the acceleration of the model, a_{\mathrm{n}} is the acceleration of the prototype, g is the gravitational acceleration, and N is the model ratio.

Thus, the centrifuge model test can better reflect the instability mechanism of the weakly sloped subgrade.

3.3 Test method

For the centrifuge model test, the first form ({\theta }_{\mathrm{1}}=\mathrm{0}, {\theta }_{\mathrm{2}}\ne \mathrm{0}) was selected. The model container size and the inclination effect were considered, and the inclination angle of the subface of the soft soil layer was set as 15°. The filling height of the prototype embankment was 8 m, and that of the model embankment was 120 mm, considering a model ratio N of 1:67. According to previous studies [16], the thickness of the soft soil layer h_m is the key parameter determining the type of slip surface. Thus, h was no less than 125 mm. The lateral profile of the model is shown in Fig. 2.

The acceleration and deceleration of the centrifuge during the start and braking times caused a change in the tangential acceleration. The start and braking times were controlled using \tau /Ng\le \mathrm{0.01} [38], where \tau is the horizontal tangential force, and the calculated start and braking times are 45 and 60 s, respectively.

3.4 Material parameters

For the centrifuge model tests, the main parameters of the materials in each soil layer were as follows.

(1) The underlying hard soil layer consisted of silty clay. Its average density was 2.6 g/cm³, the average moisture content was 24.1%, the cohesion was 41 kPa, and the angle of internal friction was 46°.

(2) Clay was used as the material of the soft soil layer. Its average density was 1.75 g/cm³, the average moisture content was 47.04%, the cohesion was 14 kPa, and the angle of internal friction was 10°.

(3) The material of the embankment fill was silty clay. The average density was 1.90 g/cm³, the average moisture content was 21.03%, the cohesion was 25 kPa, and the angle of internal friction was 25°.

3.5 Procedure

Figure 3(a) shows that cracks appeared on the embankment top at the 20 g level. With an increase in the centrifugal acceleration, the number of cracks increased gradually. Although the cracks widened and extended downward, the distribution range was always within the upper half of the embankment. In comparison, no cracks appeared on the top of the embankment of the weakly horizontal subgrade under the same conditions.

When the centrifugal acceleration increased to 40 g, a through crack appeared at the upper edge of the embankment. Meanwhile, the left embankment settlement and the uplift at the slope foot could be observed during the tests. The weakly sloped subgrade model formed a continuous slip surface. In the weakly horizontal subgrade model, the crack was formed near the center of the top surface of the embankment.

The weakly sloped subgrade became unstable and failed at the 60 g level (Fig. 3(b)). By comparison, the weakly horizontal subgrade lost its slope stability at the 70 g level.

3.6 Centrifuge model test results

Based on the plane strain theory [39], three dimensional experiments can be converted into two-dimensional cases. The total displacement and deformation of the weakly sloped subgrade after the tests are shown in Fig. 4. The slip surface is approximately circular in Fig. 4. The direction of the arrow represents the displacement direction, and the arrow length represents the displacement magnitude. Point A represents the maximum horizontal displacement, point B represents the maximum settlement, and point C represents the maximum uplift. In addition, point D represents the slope foot, point E represents the top of the lateral slope, point F represents the starting point of the slip surface, and point G represents the ending point of the slip surface.

Compared to the weakly horizontal subgrade, the deformation and instability characteristics of the weakly sloped subgrade were asymmetrical. The maximum horizontal displacement point existed at the slope foot, and the depth from point A was approximately 180 mm. The maximum vertical displacement occurred at point C. The horizontal distance between point C and the foot of the slope was approximately 160 mm. The vertical displacement occurred as a settlement within the width range of the embankment, and point B was close to the slope foot. The sliding point F crossed the centerline of the embankment and was located on the top surface of the upper half of the embankment. However, the sliding point G was further away from the slope foot compared to the weakly horizontal subgrade. Cracks initially appeared on the embankment top and gradually widened to the slope shoulder. The failure patterns of the weakly sloped subgrade and the weakly horizontal subgrade were different in four aspects: the slip surface moved up, the superficial area of the sliding body increased, the width of the cracks on the subgrade clearly increased, and the uplift height of the slope decreased. In summary, as the angle of inclination of the soft soil layer increased, the entire embankment load significantly affected the sliding instability of the weakly sloped subgrade. The failure patterns of the weakly sloped subgrade and the weakly horizontal subgrade were clearly different.

4 Numerical simulation of stability

Although centrifuge tests have certain advantages, they cannot be widely used in engineering practice because of their high cost. A numerical simulation was conducted to quantitatively compare the stability between the weakly horizontal subgrade and the weakly sloped subgrade, as described in this section.

4.1 Calculation methods and parameters

The strength-reduction finite element method is more accurate than the traditional slice method in calculating the slope stability [40]. Therefore, finite element software (PLAXIS) was used in this study to analyze the slope stability of the subgrades. The calculation models of the horizontal and weakly sloped subgrades are shown in Figs. 5 and 6, respectively. The horizontal and vertical displacements of the bottoms of the models were fixed, as were the horizontal and vertical displacements of the bottoms of the models. The horizontal displacement of both sides of the models were limited, whereas the vertical displacement of both sides of the models were not.

To ensure the calculation accuracy, the Mohr–Coulomb model was used for the soil material. The specific calculation parameters are listed in Table 1.

4.2 Analysis of calculation results

The total displacement contour maps of the weakly horizontal and weakly sloped subgrades are shown in Figs. 7 and 8, respectively. Based on the displacement contour maps, the deformation of the weakly sloped subgrade was asymmetrical, whereas that of the weakly horizontal subgrade was symmetric. The simulation results were consistent with the experimental results presented in Section 2.

The conditions of weakly sloped subgrade were divided into two types to investigate the safety factor of the weakly sloped subgrade. Under the first condition, the embankment height was constant, and the angle of inclination of the soft soil layer was changed. Under the second condition, the angle of inclination of the soft soil layer was constant, and the embankment height of different subgrades was varied.

Condition 1: The embankment height was 4 m, and {\theta }_{\mathrm{1}} and {\theta }_{\mathrm{2}} ranged from 0° to 15°. The calculated results are listed in Table 2.

The calculation results show that the safety factor has a negative correlation with the inclination of the soft soil layer when the embankment height remains unchanged.

Condition 2: The inclination angle of the soft soil layer remained unchanged, while the embankment height was varied from 2 to 8 m. A total of 96 cases were calculated, and the calculated results for {\theta }_{\mathrm{1}}=\mathrm{0} and {\theta }_{\mathrm{2}}\ne \mathrm{0} are listed in Table 3.

The results showed that the safety factor negatively correlated with the embankment height when the inclination angle of the soft soil remained unchanged.

5 Determination of influential and additional forces

Based on the centrifuge model test and numerical analysis results, the weakly sloped subgrade was more prone to instability than the weakly horizontal subgrade under the same loading conditions. In addition, the weakly sloped subgrade had an asymmetrical structure. Hence, the two novel concepts, i.e., influential force and additional force, were introduced in this study to accurately model the stability behavior of weakly sloped subgrades.

5.1 Influential force

Potential energy is the energy of an interacting object owing to its relative position. Because of the inclined plane, the weakly sloped subgrade loses its symmetry, and a potential energy difference exists in the soil above the inclined plane. The slip surface of the weakly sloped subgrade is assumed to be circular, and the sliding body satisfies the assumptions for the slice method. Therefore, it is necessary to introduce the influential force concept to evaluate the effect of the potential energy difference.

The lowest point of the circular slip surface is assumed to have zero potential energy. This zero potential energy point is denoted as “O” , and the barycenter point of the soil slice is M (Fig. 9). The potential energy E_P of the soil slice is expressed as follows:

(3)E_{\mathrm{p}}=-{\displaystyle {\int }_{\stackrel{\cap }{OM}}F\mathrm{d}r},

where F is the “influential force” acting on the soil slice, and Fdr is the elementary work applied along the direction of the micro-virtual displacement dr. Owing to the effects of gravity, the potential energy of the soil slice can be expressed as follows:

(4)E_{\mathrm{p}}=mg{h}_{\mathrm{E}},

where m is the mass of the soil slice, g is the gravitational acceleration (9.8 m/s²), and h_E is the horizontal height difference between M and O.

Through a deviation of Eqs. (3) and (4) and combining them into force F, we have the following:

(5)-{\displaystyle \frac{\mathrm{d}{E}_{\mathrm{p}}}{\mathrm{d}r}}=-F\Rightarrow mg{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{E}}}{\mathrm{d}r}}=-F.

In this paper, the “influential force” is denoted as F_E. The sliding body is assumed to be a rigid body, and the gravitational potential energy is a scalar quantity. Therefore, F_E can be considered to act on the slip surface, and its direction is consistent with that of the micro-virtual displacement dr. The direction of action is downward along the inclined plane {\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}.

5.2 Additional force

Owing to the inclined plane, the lateral embankment of the sliding body has an influence on the sliding body. The additional force is based on a comparison of the weakly horizontal subgrade and weakly sloped subgrade. The distribution form, acting position, and acting direction of the additional force F_{\mathrm{\theta }} are unknown. To identify the most reasonable additional force, the corresponding calculation program is developed based on the Bishop method in calculating the safety factor.

5.2.1 Introduction of additional force

For the weakly horizontal subgrade (Fig. 10(a)), the potential energy difference did not affect the entire subgrade. Moreover, the gravity acting on the soil, which was in the shadow section, did not affect the sliding body. Only the gravity load influenced the sliding body.

For the weakly sloped subgrade (Fig. 10(b)), as the subgrade lost its symmetry, the potential energy difference caused the soil along the inclined plane to slide downward.

The component force of the gravity of the lateral embankment (i.e., the soil in the shadow part) acted on the sliding body. The component force also facilitated a sliding, except for the influence of the gravity load of the sliding body and the potential energy difference. Therefore, the introduction of an additional force was necessary.

5.2.2 Distribution of additional force

The weight of the lateral embankment of the sliding body could be obtained based on the embankment size, inclination angle {\theta }_{\mathrm{1}}, center position, radius of the circular slip surface, and subgrade unit weight. Because the component force of lateral embankment gravity facilitated the sliding of the sliding body, the action range of the component force was supposed to be at the right side of the sliding body (bounded by the lowest point of the circular slip surface).

The acting position (1/3 above the bottom of each soil slice) and direction (along the direction \theta \mathrm{=}{\theta }_{\mathrm{1}}\mathrm{+}{\theta }_{\mathrm{2}}) of the additional force remained unchanged during the determination of the reasonable distribution form. The filling height and center coordinate of the circular slip surface significantly influences the safety factor, and the embankment filling height was set to 5 m. The safety factors of the four distribution forms are listed in Table 4.

It should be noted that the center coordinates of the slip surface for the calculation were (0, 3.5 m). In comparison, the embankment filling height was 5 m.

The expression of error is expressed as follows:

(6)e={\displaystyle \frac{\left|F_{\mathrm{s}}-F_{\mathrm{s}\mathrm{1}}\right|}{F_{\mathrm{s}\mathrm{1}}}},

Where e is the error, F_{\mathrm{S}}is the safety factor, and F_{\mathrm{S1}} is the average safety factor.

The calculated results showed that different distribution forms had a minimal influence on the safety factor, and the error between the calculation results and the average value was lower than 0.5%. The additional force was assumed to be uniformly distributed to simplify the calculation procedure.

5.2.3 Acting position of additional force

The acting position \lambda is expressed as the ratio of the vertical distance from the action point to the bottom soil slice surface and the length of the long vertical side of the soil slice. To determine the acting position of the additional force, the distribution form and acting direction of the additional force were kept constant. The safety factors and errors between the calculation results and finite element results are listed in Table 5.

Based on the calculation results, the different acting positions of the additional force significantly influenced the safety factor. The relationship between the acting positions of the additional force and the safety factor was linear.

When the acting position \lambda was zero, the error between the calculation result and the finite element result (safety factor of 1.250) was 4.6%. When the acting position \lambda was 1/2, the error was 2.6%. The error was lower than 1.0% when the acting position \lambda was 1/3. The acting position \lambda was adopted as 1/3 to simplify the calculation procedure.

Remarkably, the additional force was still assumed to act on the circular slip surface in deriving the modified Bishop method. The additional force was multiplied by a correction coefficient \eta and is expressed as follows:

(7)\eta \mathrm{=}{\displaystyle \frac{\sqrt{{(\left|R-Y_i\right|-\lambda (h-Y_i))}^{{}_{\mathrm{2}}}+X_{{}_i}^{\mathrm{2}}}}{R}},

where (X_i,Y_i) is the coordinate of the intersection point between the subface of the soil slice 'i' and the long vertical side, R is the radius of the circular slip surface, \lambda is the acting position, \eta is the correction coefficient.

5.2.4 Acting direction of additional force

To determine the acting direction of the additional force, the distribution form (uniform distribution) and the acting position (\lambda = 1/3) of the additional force were fixed. The safety factor corresponding to different acting directions and the errors between the calculation and finite element results are shown in Table 6.

Based on the calculated results, the direction of the additional force significantly influenced the safety factor, and the variation showed a linear trend. When the acting direction angle was smaller than {\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}, the safety factor was lower than that of the finite element result (safety factor= 1.250). By contrast, when the acting direction angle was larger than {\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}, the safety factor was higher than that of the finite element result. The safety factor was closest to that of the finite element when the acting direction angle was {\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}, and the error was less than 1.0%.

From these results, the distribution form of the additional force was uniform, the acting position was 1/3, and the acting direction angle was {\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}. The total error was lower than 2.5%.

6 Modified Bishop method

Based on the results presented in Section 4, the force analysis of a weakly sloped subgrade was conducted. A modified Bishop method reflecting different forms of the weakly sloped subgrade could be established through a force analysis and formula derivation.

6.1 Force analysis

According to the Bishop method, the sliding body in the weakly sloped subgrade is divided into vertical soil slices (Fig. 11). Three assumptions were made, i.e., the slip surface is circular, the safety factors of the different soil slices are equal, and the lateral force on the soil slice is considered. However, the modified Bishop method also considered both the additional force and the influential force.

In Fig. 11, the lowest point is point “O.” The soil slices of the circular slip surface near the slope toe are numbered from 1 to n_m, and soil slices near the upper slope are numbered from n_{m+ 1} to n. Soil slice i represents any of the soil slices numbered from n_{m+ 1} to n, and was selected as the research object. In Fig. 12, {\alpha }_i is the normal angle at the midpoint of the bottom of soil slice i, and b_i is the secant length of the arc at the bottom of i. Points 1 and 2 denote the positions of the influential force and additional force, respectively.

For any soil slice i ranging from n_{m+ 1} to n, five forces acted on i, i.e., W_i, F_{Ei}, F_{\theta i}, T_i, and N_i. The gravity of soil slice i, W_i, included W_{\mathrm{1}i}(the embankment part) and W_{\mathrm{2}i}(the part below the original ground surface), and W_i is directed downward vertically. The influential force F_{Ei} was derived from the potential energy difference, and its direction was {\theta }_{\mathrm{1}} + {\theta }_{\mathrm{2}} downward along the inclined plane. The additional force F_{\theta i} was derived from the weight of the lateral embankment of the sliding body, and its direction was {\theta }_{\mathrm{1}} + {\theta }_{\mathrm{2}} downward along the inclined plane. The reaction force acting on the slip surface could be split into the tangential reaction force (anti-slide force) T_i and the normal reaction force N_i. The direction of T_i was opposite to the sliding direction, and N_i was directed toward the center of the circle.

For any soil slice j ranging from 1 to n_m, an upward trend was shown throughout the entire sliding process; hence, the effects of the additional force and influential force were neglected. The forces acting on the soil slice were the self-weight of the soil slice and the reaction force acting on the circular slip surface.

In calculating the safety factor, the gravity of the soil slice W_i could be determined from its geometric size and bulk density. The anti-slide force T_i could be obtained by applying the Coulomb strength theory, and the normal reaction force N_i could be obtained from the normal pressure. The methods of computing the influential force F_{Ei} and additional force F_{\theta i} are described in Sections 6.2 and 6.3.

6.2 Calculation of additional force

Based on the acting direction and position of the additional force, the additional force equation was derived as follows:

(8)F_{\theta }=W_{3k}\sin \theta ,

where W₃ is the gravity of the lateral embankment of the sliding body, and k is the form of the weakly sloped subgrade. The value of W₃ can be obtained by the embankment area S, surface angle {\theta }_{\mathrm{1}}, the center coordinates of the circular slip surface, the radius, and the bulk density of the embankment filling.

The additional force acting on soil slice i is determined using the following equation:

(9)F_{\theta i}={\displaystyle \frac{W_{\mathrm{3}k}}{n-n_{m\mathrm{+1}}}}\sin \theta ,

The tangential and normal components of the additional force along the circular slip surface were obtained using the following expressions:

(10)T_{\mathrm{\theta }i}=F_{\mathrm{\theta }i}\cos ({\alpha }_i\mathrm{-}\theta )={\displaystyle \frac{W_{\mathrm{3k}}}{n-n_{\mathrm{m+1}}}}\sin \theta \cos ({\alpha }_i\mathrm{-}\theta ),

(11)N_{\mathrm{\theta }i}=F_{\mathrm{\theta }i}\sin ({\alpha }_i\mathrm{-}\theta )={\displaystyle \frac{W_{\mathrm{3k}}}{n-n_{\mathrm{m+1}}}}\sin \theta \sin ({\alpha }_i\mathrm{-}\theta ).

6.3 Calculation of influential force

Based on the method for determining the influential force, the influential force acting on each soil slice in the sliding body is expressed as follows:

(12)F_{\mathrm{E}i}=W_{\mathrm{1}\mathrm{k}i}\sin \theta ,

where W_1ki is the gravity of soil slice i. W_1ki can be determined from the size of i, the surface inclination angle {\theta }_{\mathrm{1}}, and the bulk density of the embankment filling.

The tangential and normal components of the influential force along the circular slip surface were obtained using the following expressions:

(13)T_{\mathrm{E}i}=F_{\mathrm{E}i}\cos ({\alpha }_i\mathrm{-}\theta )=W_{\mathrm{1k}i}\sin \theta \cos ({\alpha }_i\mathrm{-}\theta ),

(14)N_{\mathrm{E}i}=F_{\mathrm{E}i}\sin ({\alpha }_i\mathrm{-}\theta )=W_{\mathrm{1k}i}\sin \theta \sin ({\alpha }_i\mathrm{-}\theta ).

6.4 Derivation of safety factor equation

The cohesive force and internal friction angle of the soft soil at the bottom of soil slice i are denoted as c_i and {\phi }_i, and r₁ and r₂ are the bulk densities of the soft soil and embankment filling, respectively.

(1) For soil slices numbered n_{m+ 1} to n, according to the vertical and horizontal static equilibrium conditions of soil slice i, the following equation can be obtained:

(15)N_i\cos {\alpha }_i=W_i+\Delta X_i-T_i\sin {\alpha }_i+(F_{\mathrm{\theta }i}\mathrm{+}{F}_{\mathrm{E}}{}_i)\sin \theta ,

(16)\Delta E_i=N_i\sin {\alpha }_i-T_i\cos {\alpha }_i+(F_{\mathrm{\theta }i}+F_{\mathrm{E}i})\cos \theta .

WhereN_iis the normal reaction force acting on the slip surface, {\alpha }_iis the normal angle at the midpoint of the bottom of soil slice i, W_iis the gravity of the soil slice, T_i is the tangential reaction force acting on the slip surface, \Delta X_iis the difference of tangential forces on both sides of the soil slice.\Delta E_iis the difference of normal forces on both sides of the soil slice.

Equation (16) is substituted into Eq. (15), as follows:

(17)\Delta E_i=(W_i+\Delta X_i)\tan {\alpha }_i-T_i\sec {\alpha }_i\mathrm{+(}{F}_{\mathrm{\theta }i}+F_{\mathrm{E}i})(\sin \theta \tan {\alpha }_i+\cos \theta ).

(2) For soil slices numbered 1 to n_m, according to the vertical and horizontal static equilibrium conditions of the soil slice j, the following equation could be obtained:

(18)N_i\cos {\alpha }_i=W_i+\Delta X_i-T_i\sin {\alpha }_i,

(19)\Delta E_i=N_i\sin {\alpha }_i-T_i\cos {\alpha }_i\mathrm{.}

Equation (19) is substituted into Eq. (18) as follows:

(20)\Delta E_i=(W_i+\Delta X_i)\tan {\alpha }_i-T_i\sec {\alpha }_i.

The safety factor of the entire slope is denoted as F_{\mathrm{s}}. For soil slices numbered n_{m+ 1} to n (the soil slices numbered 1 to n_m have the same value), the following equation could be obtained:

(21)T_i={\displaystyle \frac{T_{\mathrm{f}i}}{F_{\mathrm{s}}}}={\displaystyle \frac{b_i{\tau }_{\mathrm{f}i}}{F_{\mathrm{s}}}},

and the expression for {\tau }_{\mathrm{f}i} could be calculated using the Mohr–Coulomb strength theory:

(22){\tau }_{\mathrm{f}i}=c_i+{\displaystyle \frac{N_i}{b_i}}\tan {\phi }_i.

Substituting Eq. (22) into Eq. (21) yields the following expression:

(23)T_i={\displaystyle \frac{c_i{b}_i+N_i\tan {\phi }_i}{F_{\mathrm{s}}}}.

Equation (23) is substituted into Eq. (15) to obtain

(24)T_i={\displaystyle \frac{c_i{b}_i+(W_i+\Delta X_i+(F_{\mathrm{\theta }i}+F_{\mathrm{E}i})\sin \theta )\tan {\phi }_i}{F_{\mathrm{s}}{m}_i}}.

m_i=\cos {\alpha }_i+{\displaystyle \frac{\sin {\alpha }_i\tan {\phi }_i}{F_{\mathrm{s}}}}.

Assuming that \Delta X_i\mathrm{=}\mathrm{0}, the moment of forces acting on each soil slice in the entire sliding body about the center coordinates of the slip surface was determined as follows:

(25){\displaystyle \sum _{i=\mathrm{1}}^n{W}_{i}R\sin {\alpha }_i}+{\displaystyle \sum _{i=m+\mathrm{1}}^n（T_{\mathrm{\theta }i}\eta R+T_{\mathrm{E}i}}R)-{\displaystyle \sum _{i=\mathrm{1}}^n{T}_{i}R}=\mathrm{0.}

Equations (24) and (25) were combined to derive the safety factor equation of the modified Bishop method as follows:

(26)F_{\mathrm{S}}={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{m_i}}(c_i{b}_i+W_i\tan {\phi }_i)+{\displaystyle \sum _{i=m+\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{m_i}}\left({\displaystyle \frac{W_{\mathrm{3k}}}{(n-n_{\mathrm{m+1}})}}+W_{\mathrm{1}\mathrm{k}i}\right){\sin }^{\mathrm{2}}\theta \tan {\phi }_i}}}{{\displaystyle \sum _{i=\mathrm{1}}^n{W}_i}\cos {\alpha }_i+{\displaystyle \sum _{i=m+\mathrm{1}}^n\left({\displaystyle \frac{\eta W_{\mathrm{3k}}}{(n-n_{\mathrm{m+1}})}}+W_{\mathrm{1}ki}\right)\sin \theta \cos ({\alpha }_i\mathrm{-}\theta )}}},

Similarly, the safety factor expressed by the effective stress can be expressed by

(27)F_{\mathrm{S}}={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{m_i}}(c_i^{\text{'}}{b}_i+(W_i-u_i{b}_i)\tan {\phi }_i^{\text{'}})+{\displaystyle \sum _{i=m+\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{m_i}}\left({\displaystyle \frac{W_{\mathrm{3k}}}{(n-n_{\mathrm{m+1}})}}+W_{\mathrm{1}\mathrm{k}i}\right){\sin }^{\mathrm{2}}\theta \tan {\phi }_i^{\text{'}}}}}{{\displaystyle \sum _{i=\mathrm{1}}^n{W}_i}\cos {\alpha }_i+{\displaystyle \sum _{i=m+\mathrm{1}}^n\left({\displaystyle \frac{\eta W_{\mathrm{3k}}}{(n-n_{\mathrm{m+1}})}}+W_{\mathrm{1}ki}\right)\sin \theta \cos ({\alpha }_i\mathrm{-}\theta )}}},

6.5 Calculation steps

The process for calculating the modified Bishop is as follows.

(1) Determination of the soil parameters: The soil parameters include the weights of the embankment filling and the soft soil layer \gamma, cohesion c, and friction angle \phi.

(2) Determination of geometric dimensions: The geometric dimensions include the center coordinate of the circular slip surface, radius of the circular slip surface R, inclination angle of the embankment bottom {\theta }_{\mathrm{1}}, inclination of the bottom of the soft soil layer {\theta }_{\mathrm{2}}, thickness of the soft soil layer at the center of the embankment, and lateral embankment area of the sliding body S.

(3) Division of soil slice: This step includes dividing the sliding body into n soil slices. Here, {\alpha }_i is the angle normal to the midpoint of the soil slice bottom, and l_i is the secant length of the arc at the soil slice bottom. The lowest point of the circular slip surface was selected as the control point of the sliding body, the additional force acted on the right side of each soil slice (1/3 above the bottom surface), and the acting direction was inclined at {\theta }_{\mathrm{1}} + {\theta }_{\mathrm{2}}.

(4) Calculation processing: A force balance was established for all soil slices and the entire sliding body. According to Eq. (24), F_{\mathrm{s}}=\mathrm{1} was first assumed for the calculation, and a new value F_{\mathrm{s}} was obtained. The iteration process was continued until two contiguous F_{\mathrm{s}} values were close.

(5) Result processing: When \left|F_{\mathrm{s}i}-F_{\mathrm{s}i-\mathrm{1}}\right|\le \mathrm{0.01}, F_{\mathrm{s}} is convergent, and the final safety factor can be determined. When m_i\le \mathrm{0.2}, the obtained safety factor will have a certain error, and other calculation methods should be selected.

To avoid heavy applications of manual modeling, a computer program was written using the above steps. The computational program is expected to be convenient for creating a subgrade design.

7 Engineering example

In this section, an example is presented to illustrate the applicability of the modified Bishop method for assessing the stability of weakly sloped subgrades.

7.1 Case study

The Xingning–Shanwei Expressway (A1 section) is located in the mountainous and hilly area of Guang Dong Province. Owing to the special geological landform, most of the embankments along the Xingning–Shanwei expressway have a weakly sloped subgrade. The cross-section of a simplified calculation model for typical weakly sloped subgrades is depicted in Fig. 13.

To reduce the boundary effect on subgrade instability, the calculation width was 120 m and the calculation height was 35 m. Based on the cross-section of the simplified calculation model, the top width of the subgrade was 26 m, and the filling height at the center of the subgrade was 7.5 m. The first slope ratio of the subgrade was 1:1.75, and the second slope ratio was 1:1.5. The bottom-left corner of the calculation model was adopted as the origin of coordinates, and the horizontal coordinate X of the subgrade centerline was 44.5 m. The slope ratio of the embankment bottom was 1:8, and the slope ratio of the bottom of the soft soil layer was 1:4. Therefore, {\theta }_{\mathrm{1}} = 7.125°, and {\theta }_{\mathrm{2}} = 14.036°. The soft soil layer thicknesses varied from 3 to 18 m within the model width range and 6 to 12.7 m within the width range of the embankment bottom surface. The thickness of the soft soil layer at the central line of the subgrade was approximately 8.5 m.

According to the engineering geological survey report, the recommended values of design parameters for each soil layer and the required parameters for the finite-element stability analysis are listed in Table 7.

7.2 Results and analysis

Different slice analysis methods, the modified Bishop method, and the finite element method were used to analyze the stability of a typical weakly sloped subgrade.

(1) Slice methods and modified Bishop method

The sliding body was divided into 20 soil slices, and the lowest point of the circular slip surface was located at the bottom-right endpoint of the eighth soil slice. Calculation results obtained using different slice methods and the modified Bishop method are listed in Table 8.

(2) Finite element method

The global sparse density of the grid was set to a medium level. The horizontal and vertical displacements at the bottoms of the models were fixed, and the horizontal displacement of both sides of the models was not limited. The grids were generated automatically using PLAXIS (Fig. 14).

The soil grid unit was a 15-node triangular unit, which consisted of 856 units, 7049 nodes, and 10272 stress points. The Mohr–Coulomb model was established for the soil material.

Based on the strength reduction method, the Phi-c reduction calculation type was adopted to calculate the safety factor. The calculation step length was set to 0.1, and the number of steps was controlled based on the stability of the safety factor. In addition, the displacement generated by the embankment load was ignored. After the embankment load was calculated, the safety factor of the entire subgrade could be solved. The horizontal contour map is shown in Fig. 15, and the results obtained using the finite element method are listed in Table 8.

Based on the calculation results, the safety factor calculated using different slice methods was 32.2% and 36.5% higher than those of the modified Bishop and finite element strength-reduction approaches, respectively. The error between the modified Bishop and finite element strength-reduction methods was lower than 3.5%. The safety factors calculated using the Swedish and Bishop methods were close to 1, and the results were conservative. Moreover, the safety factors calculated using the modified Bishop and finite element strength-reduction methods were close to 0.7, and the results indicate that the stability of the weakly sloped subgrade was low. Therefore, reinforced concrete piles should be used in a subgrade design to strengthen such subgrades.

Overall, the results showed that the modified Bishop method precisely reflects the characteristics of weakly sloped subgrades, and the obtained safety factor is close to the actual value.

8 Conclusion

Based on numerical calculations, the modified Bishop method developed in this study accurately reflects the stability characteristics of weakly sloped subgrades. The conclusions of this study are summarized as follows.

(1) Centrifuge model test and numerical simulation results showed that the sliding form of the weakly sloped subgrade was different from that of the weakly horizontal subgrade. When the slice method is used to analyze the stability of weakly sloped subgrades, the introduction of an influential force and an additional force is necessary. The influential force reflects the effect of the potential energy difference, whereas the additional force reflects the asymmetry effect.

(2) The sliding body is assumed to be a rigid body, and the gravitational potential energy is a scalar quantity. Therefore, the influential force F_{\mathrm{E}} was considered to act on the slip surface in the direction of ({\theta }_{\mathrm{1}} + {\theta }_{\mathrm{2}}) along the inclined surface. The range of action was from the upper subgrade side bounded by the lowest point of the circular slip surface.

(3) From the calculation results of the self-made program, the distribution form was uniform, the acting position of additional force F_{\mathrm{\theta }} was 1/3, and the acting direction was {\theta }_{\mathrm{1}} + {\theta }_{\mathrm{2}}. The additional force acting on the upper subgrade side was bounded by the lowest point of the slip surface. Based on these conditions, the calculated safety factor was the closest to that of the finite element result, and the comprehensive error was lower than 2.5%.

(4) An engineering example of the stability analysis demonstrated that the proposed modified Bishop method is suitable for different inclined forms. Compared to the Bishop and Swedish methods, the results obtained using the modified Bishop method were closer to the finite element results. The modified Bishop method reflected the easy sliding characteristic of the weakly sloped subgrade, and the results approached the actual value.

In a future study, more sliding surfaces should be considered. Furthermore, the crack patterns formed during centrifuge model tests should be investigated and compared with those predicted using phase-field models [41–43].

References

[1] Bishop A W. The use of slip circle in the stability analysis of slopes. Geotechnique, 1955, 5(1): 7–17

[2] Morgenstern N R, Price V E. The analysis of the stability of general slip surface. Geotechnique, 1965, 15(1): 79–93

[3] Spencer E. A method of analysis of the stability of embankments assuming parallel inter-slice forces. Geotechnique, 1967, 17(1): 11–26

[4] Sarma S K. Stability analysis of embankments and slopes. Journal of Geotechnical Engineering, 1979, 105(12): 1511–1524

[5] Chen Z Y, Morgenstern N R. Extensions to the generalized method of slices for stability analysis. Canadian Geotechnical Journal, 1983, 20(1): 104–119

[6] You C L, Zhao C G, Zhang H C, Liu H W. Study on construction test of embankment on soft clay of plateau slope. Chinese Journal of Geotechnical Engineering, 2002, 24(4): 503–508 (in Chinese)

[7] Zhao S, Zheng Y, Shi W, Wang J. Analysis on safety factor of slope by strength reduction FEM. Chinese Journal of Geotechnical Engineering, 2002, 24(3): 343–347 (in Chinese)

[8] Hu Y G, Luo Q, Zhang L, Cheng Y M. Deformation characteristics analysis of slope soft soil foundation treatment with mixed-in-place pile by centrifugal model tests. Rock and Soil Mechanics, 2010, 31(7): 2207–2213 (in Chinese)

[9] Lian Z Y, Han G C, Kong X J. Stability analysis of excavation by strength reduction FEM. Chinese Journal of Geotechnical Engineering, 2001, 23(4): 407–411 (in Chinese)

[10] Liu X, Sheng K, Li Z, Gan L, Shan H, Hong B. Experimental research on foamed mixture lightweight soil mixed with fly-ash and quicklime as backfill material behind abutments of expressway bridge. Advances in Materials Science and Engineering, 2017: 1–11

[11] Liu X, Sheng K, Shan H. Discussion of “Revised soil classification system for coarse-fine mixtures” by Junghee Park and J. Carlos Santamarina. Journal of Geotechnical and Geoenvironmental Engineering, 2018, 144(8): 07018017

[12] Baker R, Garber M. Theoretical analysis of the stability of slopes. Geotechnique, 1978, 28(4): 395–411

[13] Wu J, Cheng Q, Liang X, Cao J. Stability analysis of a high loess slope reinforced by the combination system of soil nails and stabilization piles. Frontiers of Structural and Civil Engineering, 2014, 8(3): 252–259

[14] Cheng Y M. Location of critical failure surface and some further studies on slope stability analysis. Computers and Geotechnics, 2003, 30(3): 255–267

[15] Jiang X, Qiu Y J, Wei Y X. Engineering behavior of subgrade embankments on sloped weak ground based on strength reduction FEM. Chinese Journal of Geotechnical Engineering, 2007, 29(4): 622–627 (in Chinese)

[16] Jiang X, Zhu Q J, Jiang Y, Qiu Y. Stability analysis of embankment over sloped weak ground with limit equilibrium method based on reliability. Journal of Railway Science and Engineering, 2013, 10(2): 47–55 (in Chinese)

[17] Liu X, Sheng K, Hua J, Hong B, Zhu J. Utilization of high liquid limit soil as subgrade materials with pack-and- cover method in road embankment construction. International Journal of Civil Engineering, 2015, 13(3): 167–174

[18] Liu Y, He Z, Li B, Yang Q. Slope stability analysis based on a multigrid method using a nonlinear 3D finite element model. Frontiers of Structural and Civil Engineering, 2013, 7(1): 24–31

[19] Cheng Y M, Yip C J. Three-dimensional asymmetrical slope stability analysis extension of Bishop’s, Janbu’s, and Morgenstern-Price’s techniques. Journal of Geotechnical and Geoenvironmental Engineering, 2007, 133(12): 1544–1555

[20] Pham H T, Fredlund D G. The application of dynamic programming to slope stability analysis. Canadian Geotechnical Journal, 2003, 40(4): 830–847

[21] Zhang G, Hu Y, Wang L. Behaviour and mechanism of failure process of soil slopes. Environmental Earth Sciences, 2015, 73(4): 1701–1713

[22] Lü X L, Xue D W, Chen Q S, Zhai X L, Huang M S. Centrifuge model test and limit equilibrium analysis of the stability of municipal solid waste slopes. Bulletin of Engineering Geology and the Environment, 2019, 78(4): 3011–3021

[23] Fang H W, Chen Y F, Xu Y X. New instability criterion for stability analysis of homogeneous slopes. International Journal of Geomechanics, 2020, 20(5): 04020034

[24] Zhao Y, Zhang G, Hu D, Han Y. Centrifuge model test study on failure behavior of soil slopes overlying the bedrock. International Journal of Geomechanics, 2018, 18(11): 04018144

[25] Lin C, Li H, Weng M. Discrete element simulation of the dynamic response of a dip slope under shaking table tests. Engineering Geology, 2018, 243: 168–180

[26] Weng M, Chen T, Tsai S. Modeling scale effects on consequent slope deformation by centrifuge model tests and the discrete element method. Landslides, 2017, 14(3): 981–993

[27] Liang T, Knappett J. Centrifuge modelling of the influence of slope height on the seismic performance of rooted slopes. Geotechnique, 2017, 67(10): 855–869

[28] Liu J L, Chen L W, Wang D L. Characters of embankment on inclined weak foundation. Rock and Soil Mechanics, 2010, 31(6): 2006–2010 (in Chinese)

[29] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69(2): 102–109

[30] Zhou S, Zhu H, Yan Z, Ju J W, Zhang L. A micromechanical study of the breakage mechanism of microcapsules in concrete using PFC2D. Construction & Building Materials, 2016, 115: 452–463

[31] Zhou S, Zhu H, Ju J W, Yan Z, Chen Q. Modeling microcapsule-enabled self-healing cementitious composite materials using discrete element method. International Journal of Damage Mechanics, 2017, 26(2): 340–357

[32] Zhou S, Zhuang X Y. Characterization of loading rate effects on the interactions between crack growth and inclusions in cementitious material. Computers,Materials and Continua, 2018, 57(3): 417–446

[33] Zhou S, Zhuang X. Micromechanical study of loading rate effects between a hole and a crack. Underground Space, 2019, 4(1): 22–30

[34] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[35] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[36] Chen L, Rabczuk T, Bordas S P A, Liu G R, Zeng K Y, Kerfriden P. Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209–212(1): 250–265

[37] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92– 93(1): 242–256

[38] Wu B Y. Soft Soil Foundation Treatment. Beijing: China Railway Publishing House, 1995, 247–248 (in Chinese)

[39] Xu Z L. A Concise Course in Elasticity. Beijing: Higher Education Press, 2002, 11–12 (in Chinese)

[40] Ugai K. A method of calculation of total safety factor of slope by elasto-plastic FEM. Soil and Foundation, 1989, 29(2): 190–195

[41] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[42] Zhou S, Rabczuk T, Zhuang X. Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies. Advances in Engineering Software, 2018, 122: 31–49

[43] Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Engineering Geology, 2018, 240: 189–203

RIGHTS & PERMISSIONS

Higher Education Press