An Analysis of Slope Stability Based on Finite Element Method and Distinct Element Method

The analysis of slope stability involves complex geological and topographical boundary conditions, nonlinear behavior of material stress-strain, coupling analysis of initial in-situ stress, water pressure and seismic load, etc., and in most cases, analytical solutions cannot be obtained. Under the background of the continuous development of computer and calculation method, the numerical analysis method represented by finite element has been gradually popularized and applied in geotechnical engineering in 1970s, and has developed into a powerful calculation and analysis tool. Among them, the finite element strength reduction method and the discrete element method are the two most widely used slope numerical analysis methods. In this paper, two typical cases, Ankang reservoir landslide and Wenma Highway slope, are simulated by the two methods. Taking Ankang reservoir landslide as the research object, this paper would use MIDAS / GTS finite element analysis software, and two-dimensional finite element numerical simulation would be carried out to study the influence of reservoir water level periodic fluctuation on the reinforcement effect of anti-slide pile. Under the condition of water saturation and water loss cycle, main material of landslide body and landslide belt, namely the strong weathered phyllite, displays obvious deterioration phenomenon, showing the trend of rapid decline first and then slow decline; after the anti-slide pile is set in the middle and front of the slope, the stability of it has been greatly improved, but with the increasement of the number of water level changes, the reinforcement effect of the anti-slide pile continues to weaken, and the weakening speed is fast at first, and then slows down. Taking the bedding slope of Wenma Highway as the research object, this pater would adopt UDEC discrete element software to simulate the deformation and failure process of the slope after excavated, and analysize the failure mechanism at the same time. The failure process of bedding slope can be divided into four stages: the formation of tension cracks caused by excavation, the expansion of cracks and the formation of deformation body, the sliding of deformation body and the accumulation of damaged rock mass at the foot of slope. Tensile failure is the main failure mode, and shear failure occurs locally. The failure of bedding slope starts from the foot of slope, which is traction sliding.


Introduction
The analysis of slope stability is the prerequisite of slope design, which determines whether the slope is unstable and how much thrust exists when the slope is unstable, so as to provide scientific basis for support structure design. However, this problem has not been properly solved, because in order to solve this problem, we must first find out the geological conditions and strength parameters of the  2 slope, and at the same time, have a scientific and reasonable analysis method [1]. For homogeneous soil slope, the traditional methods mainly include: limit equilibrium method, limit analysis method, slip line field method, etc.; for the current engineering application, the main method is limit equilibrium method, but the position and shape of sliding surface need to be known in advance. For homogeneous soil slopes, various optimization methods can be used to search for dangerous sliding surfaces. However, for rock slopes, due to the fact that there are a large number of discontinuities with different structures, occurrences and characteristics, the traditional limit equilibrium method is not qualified enough to figure out the dangerous sliding surfaces and the corresponding stability safety factors.
In 1975, the British scientist Zienkiewicz proposed to use the method of increasing load or reducing geotechnical strength to calculate the ultimate load and safety factor of geotechnical engineering [1]. In the 1980s and 1990s, it was used for slope and foundation stability analysis [2]. However, due to the lack of strict, reliable and powerful large-scale finite element program at that time, as well as the selection of strength criteria and the lack of specific operation technology, the calculation accuracy was insufficient, which was not widely adopted by geotechnical engineering.
At the end of the 20th century, many papers have been published in the world to study the solution of the stability safety factor F of homogeneous soil slope by the strength reduction method of finite element method [3,4]. Because the results of some calculation examples are close to the results of traditional method, they are gradually recognized by the academic community. Some foreign scholars think that the finite element strength reduction method has brought the analysis of slope stability into a new era. Especially in 1999, D.V. Griffith and others of Colorado Institute of mining and technology used their own finite element program to analyze the stability of homogeneous soil slope [5,6], which is different from other programs in that it can simulate the influence of water level and pore water pressure, and can also analyze the slope stability under the condition of reservoir water falling.
Distinct Element Method (DEM) was proposed by Cundall in the early 1970s [7]. At first, it mainly studied the mechanical behavior of rock and other discontinuous media. Cundall proposed the first practical DEM and used it to simulate the progressive movement of rock blocks [8]. Later, Cundall and Strack put forward a two-dimensional program to simulate particles, and the results were in accordance with the dynamic phofootlastic experiments [9]. So far, through the joint efforts of many scholars, the DEM has been greatly developed. The early DEM can only deal with the distinct stiffness block system, but it has been extended to simulate the deformed block. In 1980, Itasca Consulting Group developed UDEC and put it into the market [10]. Lorig and Brady have developed a distinct element boundary element coupled calculation program [11]. Cundall and other scientists developed a 3D distinct element program (3DEC) to simulate jointed rock mass [12,13].
Finite element strength reduction method (FEM) and distinct element method (DEM) are widely used in slope numerical analysis. Therefore, in this paper, two typical cases, Ankang reservoir landslide and Wenma Highway slope, are used for numerical simulation.

Finite Element Reduction Method
With the development of computer technology, the finite element strength reduction method has been concerned at home and abroad recently. When the slope reaches the failure state through the finite element strength reduction, the displacement on the sliding surface will change suddenly and produce large and unlimited plastic flow. The finite element program cannot find one from the finite element equations which can satisfy the static balance and the stress-strain relationship. Therefore, it is reasonable to use the convergence criterion of force and displacement as the criterion of slope failure. The calculation accuracy and influential factors of the finite element strength reduction method are analyzed in details, including yield criterion, flow rule, finite element model itself and calculation parameters. The specific measures to improve the calculation accuracy are given: the error between the stability safety factor obtained by Mohr Coulomb equal area circle yield criterion and the traditional Spencer method is about 5 %, which proves the feasibility of this method in engineering; in the case of plane strain, the mole matching DP criterion can be used; and the finite element method can also be applied to the stability analysis of rock slope. Therefore, in this paper, the strength reduction program of MIDAS/GTS finite element software is used to simulate the Ankang reservoir landslide, and the sliding surface and safety factor of the landslide are obtained.

Distinct Element Method
Distinct element belongs to the category of discontinuous mechanics, which describes the continuous elements and discontinuous elements in the research medium respectively. For example, the two basic components of rock mass, rock block and structural surface are described by the law of continuous mechanics and contact law respectively. The contact (structural plane) is the boundary of continuum (rock block), and a single continuum can be treated in the process of mechanical solution. It forms an independent object and interacts with other continuum through contact. For rock slope, the block often deforms and destroys along the structural plane, while UDEC (Universal Distinct Element Code) distinct element software (a kind of analysis program based on distinct element theory) has great advantages in dealing with discontinuous medium, which is suitable for dealing with the dynamic response of solid medium under external action, and can better simulate the failure process of discrete medium along the structural plane, such as sliding or tensile crack. Therefore, this paper uses UDEC distinct element software to simulate the bedding slope of Wenma Highway to reflect the specific process of slope deformation and failure and explore the mechanism of slope failure.

Ankang Reservoir Landslide
Ankang reservoir landslide is located in the slope section of "V" valley slope on the East Bank of Ankang reservoir, which is a typical bedrock landslide. The rear edge of the landslide is clear, which is a 2 m high landslide scarp. The shear outlet is located at the front edge of the river scarp, and the main sliding direction is about 276°. The plane shape of the landslide (figure 1) is approximately long fanshaped, with a length of 50~80 m, a width of 300~320 m, and a slope of about 25° to 35°. The landslide mass is mainly strongly weathered phyllite (Od), overlying silty clay with gravel. The average thickness is about 13 m, and the total volume is about 27.94×10 4 m 3 . According to the landslide volume classification, it belongs to medium-sized landslide; the sliding zone is located at the contact surface of strongly weathered phyllite and moderately weathered phyllite, which is nearly arcshaped; the sliding bed is moderately weathered Ordovician phyllite (Od). It is easy to form a sliding surface between the weathering boundary of strong and medium phyllite. The treatment measures are to set up a row of anti-slide piles in front of the landslide mass, and arrange drainage works on the slope.

Model Building.
Madis numerical simulation is carried out for the typical section of Ankang reservoir landslide. The negative direction of X axis is the main sliding direction (horizontal), and the negative direction of Y axis is the gravity direction (vertical). The origin is the left corner of the model. The length of the model in X direction is 148.40 m, and the height in Y direction is 81.40 m. The geological conditions in the model are mainly based on field survey data, and are mainly divided into sliding mass, sliding zone, strongly weathered phyllite and moderately weathered phyllite. The antislide pile is set in the front of the slope, which is also the area with large displacement. The post grid model has a total of 1325 elements and 1478 nodes.

Boundary Condition.
Displacement boundary: the slope is a free surface, and fixed displacement constraints are set on the right side and bottom surface. Water head boundary: referring to the reasonable water level of Ankang reservoir in each month, it is proposed that the water level changes at a constant speed between 300 m and 330 m. Therefore, taking the lower left corner of the model as

Wenma Highway Slope
The slope of K28+466~K28+758 section is about 232 m high and 324 m wide, which is located at the foot of left bank of Zagunao River. The slope is a bedding steep slope with a slope of about 44° and the lithology is phyllite with thin layer. The schistosity plane is extremely developed and the occurrence is 23°∠48°. There are two groups of joints in the slope, the occurrence of which are 237°∠ 72° and 274°∠50° respectively. At the same time, there are many weathering fractures and unloading fractures on the slope. 5 interface in the model includes original slope line, schistosity surface and joint contact surface. The influence of joint cutting and slope excavation on the overall stability is mainly reflected in the sliding of the surface rock layer of the slope. In order to facilitate the calculation during the modeling, the rock strata and joints on the surface of the slope are densified, and the interior of the slope is simplified.

Boundary Condition.
In order to balance the system of the model and meet the final requirements of the simulation, the boundary conditions of the slope model need to be limited: (1) There is no external load on the slope surface, and it is free surface; (2) Only the influence of gravity field is considered in the slope, and the influence of temperature field, tectonic stress field and seepage field is ignored; (3) The bottom surface and both sides of the slope are set as fixed, and the bottom surface (X direction) is set to make the velocity in Y direction is 0; the constraint on both sides (Y direction) is set to make the velocity in X direction is 0.

Parameter Selection.
Before the calculation and simulation analysis, it is necessary to define the contact surface. In order to make the established slope calculation model consistent with the actual slope layer structure, the equivalent treatment is made when the physical and mechanical parameters of rock block and structural plane are given, and the final parameters are shown in tables 3 and 4.

Ankang Reservoir Landslide
The calculation conditions are divided into 5 cycles: 1, 5, 10, 20 and 40. The stress-strain characteristics of slope stability and stress-strain relationship of anti-slide pile under different water saturation loss cycles are analyzed.

Stability Coefficient.
The stability coefficient of the slope with anti-slide pile begins to decrease with the increase of the number of cycles (figure 1). Between 1 and 20 times, the stability coefficient decreases approximately linearly; when 20-40 times, the downward trend slows down, and the stability coefficient is 1.06, and the slope is still in a basically stable state, but the reinforcement effect of the anti-slide pile cannot meet the design requirements.  figure 2 (b)), the maximum displacement range of the slope is larger than that of the first cycle, and it is concentrated in front of the pile. After 10 cycles (figure 2 (c)), the displacement in front of the pile extends to the back of the pile. After 20 cycles (Figure 2 (d)), the maximum displacement range in front of the pile does not further expand. After 40 cycles (figure 2 (e)), the displacement of the slope changes obviously. The displacement range of the pile before and after the pile intersects, and the range of the displacement behind the pile increases, but the maximum displacement increases slowly, which is only 0.004 m higher than the maximum displacement of 20 cycles. The maximum displacement at the back of pile increases gradually at 1 to 20 cycles, and tends to be horizontal after 20 to 40 cycles (figure 3).   4 (a)), the maximum shear strain of slope body is mainly distributed at the boundary between strongly weathered and moderately weathered phyllite at the foot of slope. At 5 cycles ( figure 4 (b)), the maximum shear strain extends upward. After 10 cycles (figure 4 (c)), the shear strain is distributed at the weathering boundary of phyllite in the middle of the slope, which is not connected. After 20 cycles (figure 4 (d)), the shear strain band behind the pile is further extended, but it is still not through. After 40 cycles (figure 4 (e)), a through shear strain band is formed behind the pile, that is, the soil behind the pile is destroyed and the slope body slides along the position. The value of the post pile shear strain band is analyzed and counted. With the increase of the number of cycles (figure 5), the post pile shear strain increases continuously. When the number of cycles is 1~20, the growth rate is fast, but when the cycle number is 20~40, the growth rate is slow.

Stress Analysis.
The maximum principal stress nephogram formed after slope excavation is shown in Figure 6, which is similar to the maximum principal stress nephogram under the initial state. It is distributed nearly parallel to the slope surface, and gradually increases from the slope surface to the interior of the slope body. However, the maximum principal stress at the foot of the slope increases from -1 MPa to -2 MPa in this state, and the range is also expanded under the influence of excavation.
It is easy to increase the sliding force of rock mass and reduce the stability of slope. The tensile stress value in the minimum principal stress nephogram did not increase, but the tensile stress range at the slope began to expand to the interior of the slope (figure 7), indicating that the scope of rock mass with tensile failure trend increased, while the minimum principal stress inside the slope did not change significantly. However, the interaction between the compressive stress in the maximum principal stress and the tensile stress in the minimum principal stress on the slope surface easily leads to new tension cracks and sliding failure. Figure 8 shows the shear stress nephogram formed after excavation. Because of the influence of excavation, the upper rock mass loses support, the shear stress at the foot of slope increases to -1 MPa, and the shear stress concentration zone extends along the schistosity towards the upper part, connecting with the joint surface and tensile fracture. Due to the free surface formed at the foot of the excavation slope, the locking capacity of the rock mass is lost, and the sliding failure is easy to occur.

Displacement Analysis.
For the bedding slope, to excavate the foot of the slope is easy to cause the rock mass to slide along the layer, and the slope is also destroyed rapidly after excavation. It can be seen from the displacement vector diagram (Figure 9) that the surface and near surface rock layers after excavation tend to slide along the structural plane. Among them, the displacement vector at the foot of the slope is the largest, and the initial displacement is 0.052 m. As time goes on, the displacement change will further increase, indicating that the sliding effect brought by excavation of slope foot is very obvious.
The horizontal displacement nephogram of the slope after excavation ( Figure 10) reveals the process from the beginning of slope failure to the end of it, and shows the actual sliding mode, sliding range and sliding distance. When step=5e 4 , the rock strata near the foot of the slope slide, and the horizontal displacement of each part of the sliding is different. The top broken rock block at the foot of the slope first displaces and the displacement is the largest, which indicates that the failure occurs at the foot of the slope, and then the upper rock layer begins to slide. Therefore, the failure mode belongs to traction sliding failure; In the stage when step=1e 6~5 .5e 6 , the rock layer continues to slide, and the scope of the sliding rock layer does not expand. The rock layer sliding occurs only in the empty part of the slope foot caused by excavation, and the whole lower slope body does not have deformation phenomenon; when step=7e 6 , the maximum horizontal displacement is 48 m and the minimum horizontal displacement is 8 m. Finally, part of the sliding mass is accumulated at the foot of the slope, and part of it remains in the upper part of the slope. (c) step=2.5e 6 (d) step=4e 6 (e) step=5.5e 6 (f) step=7e 6 Figure 10. Contour map of horizontal displacement after slope excavation.
The effect of excavation at foot on the bedding slope is very obvious. The sliding process can be divided into four stages: the formation of tension cracks caused by excavation, the expansion of tension cracks and the penetration of sliding surface, the sliding of rock mass cut by structural plane, and the accumulation of materials.

Conclusion
In this paper, the finite element strength reduction method and distinct element method are employed to simulate the Ankang reservoir landslide and Wenma Highway slope respectively, and the conclusions are as follow: (1) Using the finite element method to simulate the reservoir water level landslide, the following conclusions are obtained: after the reinforcement of anti-slide pile, the slope is relatively stable in the initial stage, the displacement and shear strain are mainly concentrated in the front of the slope, and the reinforcement effect is obvious. However, with the increase of the number of water level cycles, the displacement and plastic strain band gradually move back to the back edge of the landslide, and the maximum values are concentrated in front of the pile, and the force on the pile body at the sliding zone is the largest and its value increases slowly. The reinforcement effect of anti-slide pile has experienced two stages, that is, before the number of cycles is 20, the reinforcement effect is weakened quickly; when the number of cycles is 20 to 40, the weakening effect is slow. The reinforcement effect of antislide pile is weakening. When the water level of the reservoir changes for many years, the slope still has the possibility of damage, so there should be enough safety reserve when designing it.
(2) The distinct element method is used to simulate the bedding rock slope. According to the above analysis, the following conclusions can be obtained: the failure process of bedding slope is mainly caused by tension, and shear failure occurs at the foot of slope, resulting in the loss of locking capacity and sliding of rock mass; the interaction between the maximum principal stress and the minimum principal stress leads to the formation of many tension cracks at the slope surface, which enhances the possibility of rock mass sliding; the rock mass has the characteristics of high strength, high stability and high stability The sliding range is heavily controlled by the excavation, only the free part of the slope foot slides and does not extend into the slope.