Study on the Stability of Soil–Rock Mixture Slopes Based on the Material Point Strength Reduction Method

Abstract

In this paper, the material point strength reduction method is used to investigate the stability of soil–rock mixture (SRM) slopes and the whole process of large deformation occurring after destabilization. A comparative study with homogeneous soil slopes is conducted. First, a material point slope model with typical shapes, a homogeneous soil slope, and an SRM slope with stones of different sizes distributed inside is established. Next, gravity is linearly added to establish the initial state of the slopes. Then the material strength of the slope is discounted according to the criterion of strength discounting. The material point method (MPM) simulations of the two slopes are carried out separately until the slope’s displacement changes abruptly to determine the slope’s safety factor. The final accumulation form of the slope after the damage is studied. Finally, the deformation characteristics of the two slopes under extreme conditions are explored. The research shows that the stones are beneficial to the slope in maintaining slope stability, and due to the presence of stones, the slope presents different characteristics from the pure soil slope when damage occurs.

1. Introduction

The main methods of slope stability analysis are the limit equilibrium method, the finite element strength reduction method, as well as the reliability method, the fuzzy mathematical method, and the artificial intelligence method [1]. The limit equilibrium method takes into account the strength properties of the soil but does not take into account the soil’s constitutive model and cannot reflect the relationship between stress and strain within the slope. The finite element method (FEM) of strength reduction overcomes the difficulties encountered in the former method and has a high degree of acceptance in the field of slope stability analysis because it can take into account the non-linear constitutive relationships of geotechnical materials and deal with complex working conditions such as excavation and shoring [2,3,4,5]. However, the solution format of the FEM makes it difficult or impossible to solve the subsequent behaviors of the slope after instability [6]. Methods such as reliability methods, fuzzy mathematical methods, and artificial intelligence methods are less interpretable, the conclusions obtained by utilizing these methods are often not of sufficient physical significance, and they rely on the experience of the engineer to work well [7].

In the specific implementation of the strength reduction method, the judgment criteria are the key, but at present, there is a big controversy on the selection of the judging criteria. For example, Song et al. [8] used the displacement in the vertical direction at the top point of the slope and the lateral displacement at the foot point of the slope, whether there is a sudden variation, as the criterion for slope instability. Yingren et al. and Shangyi et al. [9,10] considered that the non-convergence of FEM should be used as the criterion for slope failure. Maotian et al. [11] and Jinlong et al. [12] considered that the plastic zone or the generalized plastic shear strain running from the foot to the top of the slope should be used as a sign of slope failure. In addition, there are significant differences in the magnitude of the stability coefficients obtained with different determinants. It is worthwhile to investigate whether a simple numerical calculation method can be used to improve the judgment criteria of the strength reduction method and if MPM is a more suitable algorithm. Shibutao et al. [13] have carried out more work in this area. They have proposed the material point strength reduction method and carried out a stability analysis of a slope and compared it with the finite element strength reduction method to illustrate the advantages of the material point strength reduction method. In this paper, the stability of SRM slopes is studied by the material point strength reduction method, using the occurrence of sudden changes in vertical displacement of the slope as a criterion for instability. The damage behaviors of the slope after instability are also analyzed, including displacement, kinetic energy changes et al.

Previous numerical simulations on SRM slopes can be divided into two main areas: one is based on stability analysis represented by the FEM. Xu et al. [14] established the SRM slope model using digital image processing technology. On this basis, the stability study was carried out using the finite element strength reduction method. They found that the distribution of shear bands presents prominent rock-surrounding characteristics, and the stability coefficient of the SRM slope is higher than that of an equivalent homogeneous soil slope. Lianheng et al. [3] proposed an SRM slope model which can consider stones of arbitrary shape and realize different block size distribution and stones content. Then, the FEM was used to simulate the stability of the SRM slope with different block sizes and stones content. The other is to use the discrete element method (DEM) to analyze the Micro-Mechanical Properties of SRM [9,15,16,17] carried out a two-dimensional particle flow simulation of soil–rock mixture near the foundation of a dam by biaxial test and direct shear test and analyzed the influence laws of rock content, shape of block stone, and confining pressure.

Non-homogeneous slopes are more widely distributed in nature than homogeneous soil slopes. Moreover, one of the most representative of these is the SRM slope [18,19], which consists of stones and soils with comprehensively varying properties (in this study, the effect of water is not considered). The traditional assumption of homogeneous slopes always leads to erroneous results in the solution of such problems. At present, research on SRM slope mainly includes on-site geological investigation [20], indoor experiments [21], numerical simulation [3,4,15,22] and so on. With the development of computer technology and the advantages of numerical simulation, such as low cost and high repeatability, computer simulation has gradually become an integral part of geotechnical engineering research. It has remained a research hotspot in the field.

Current numerical simulations of SRM are more concerned with their stability and small-scale model tests to determine the mechanical properties of SRM materials and less with the instability of SRM slopes and the large deformation phase after instability. However, it must be acknowledged that it is impossible to prevent landslides from occurring in the foreseeable future [23], so it is equally important to investigate the extent of the effects of slope failure as it is to investigate the stability of slopes. As mentioned above, the material point strength reduction method is adopted in this paper to study the stability and post-damage behaviors of SRM slopes. On the one hand, because of the significant advantages of the MPM in solving large deformation problems, e.g., simple imposition of boundary conditions, and higher solution efficiency compared to DEM, SPH, and so on [24]; on the other hand, the combination of the strength reduction method makes the MPM, a method for solving dynamic problems, equally capable of determining the safety factor of a slope. The MPM, first proposed by Sulsky et al. [25] uses a hybrid Eulerian–Lagrangian description, thus combining the advantages of both the Lagrangian and Eulerian descriptions. All material information is stored in the material point, the grid exists only as a tool for solving the equations, and deformed meshes are discarded, which ensures the reliability of the material point method for modeling large deformation problems. Over the past 20 years, the material point method has evolved considerably, extending a number of stable solution formats [26,27,28,29], and has been applied to a wide range of fields such as geotechnical engineering [30,31,32,33], biomechanics [6], aerospace [34,35,36] and so on.

This paper uses the MPM to simulate and analyze the SRM slopes. First of all, a material point model of the slope is established. One is a pure soil slope model without considering stones, and the other is an SRM slope model. Next, the mechanical parameters of the slope are reduced, the safety factor of the slope can be determined, and the stones’ contribution to the slope’s stability is explored. Finally, the whole process of the landslide is studied, and the different damage behaviors of the pure soil slope and SRM slope are compared.

2. Methodology

A brief overview of the MPM is presented before discussing the material point strength discounting method. The material point method originated from the Particle in Cell [6] method and was later extended to the field of solid mechanics by Sulsky et al., namely the MPM. The method is based on a coupled particles and grid framework to describe the simulated object. Specifically, the material points play the role of Lagrangian particles carrying all the material information, the background grid plays the role of an Eulerian grid to solve the control equations, impose boundary conditions, and so on, and the deformed grid is discarded at the end of the time step [23]. This combination of background grid and material points to jointly describe the object does not encounter problems such as grid distortion when solving large deformation problems, making the MPM an ideal numerical algorithm for large deformation problems [37]. The material point strength reduction method is similar to the finite element strength reduction method. The reduced mechanical parameters, i.e., the cohesion and the angle of internal friction angle, are brought into the calculation program, and when the displacement of the slope changes abruptly, the reduction factor at this time is taken as the stability factor of the slope.

In this paper, all formula notations satisfy the Einstein summation convention, with the subscripts $I$ and $p$ denoting the background grid nodes and material points, respectively. In the classical MPM, a continuous material is described as a series of discrete material points, as shown in Figure 1, and the density of the material can be approximated as follows:

$$\rho (x_i)={\displaystyle \sum _p{m}_p\delta (x_i-x_{ip})}$$

(1)

where $m_p$ denotes the mass of the matter point, $\delta$ denotes the Dirichlet function, and $x_{ip}$ and $x_i$ denote the coordinates of the matter points and grid nodes, respectively.

Based on the theory of continuous medium mechanics, the MPM uses a weak form of discretization of the equations of motion [25]. During the simulation, all material information is stored in the material point, which needs to be updated with the aid of the background grid, and thus the relationship between the material point and the background grid needs to be established during the Lagrangian phase [38]. In the MPM, this relationship is established by means of a shape function built on the background grid nodes. According to the choice of shape function, the MPM can be divided into the classical MMP (using piece-wise linear shape functions), the generalized interpolation material point method (GIMP) (using generalized interpolation shape functions) [28], the B-spline material point method (using B-spline shape functions) [39,40] and so on.

In this paper, the uGIMP shape function is employed to ensure that the calculation has sufficient stability and good computational efficiency. Assuming that there is a field function $f(x)$ distributed over the domain, the value of the function on the grid node $I$ can be expressed as:

$$f_I={\displaystyle \underset{\mathsf{\Omega }}{\int }f(x)S_{I}d\mathsf{\Omega }={\displaystyle \sum _p{f}_p{S}_{Ip}{V}_p}}$$

(2)

where $S_I$ denotes the generalized interpolation shape function, $f_p$, $S_{Ip}$, $V_p$ denote the values of the field function, shape function, and volume at the material point p, respectively. Due to the use of generalized interpolation shape functions, the cell-crossing error [27] of the MPM is significantly reduced, which is very beneficial for the stability of the numerical method.

At the beginning of each calculation step of the MPM, the information on the material points will be mapped to the grid nodes. The mass, momentum, and velocity of the grid nodes can be represented as:

$$m_I={\displaystyle \sum _p{m}_p{S}_{Ip}}$$

(3)

$$P_{Ii}={\displaystyle \sum _p{m}_p{v}_{pi}{S}_{Ip}}$$

(4)

$$v_{Ii}=\frac{P_{Ii}}{m_I}$$

(5)

where, $m_p$ and $v_{pi}$ denote the mass and velocity of the material points, respectively.

Use the node velocity gradient to solve the strain rate ${\dot{\epsilon }}_{ijp}^{n-1/2}$ and rotation rate ${\mathsf{\Omega }}_{ijp}^{n-1/2}$ of MPs, and then update the density and stress of each MP:

$${\dot{\epsilon }}_{ijp}^{n-1/2}={\displaystyle \sum _I\frac{1}{2}}({\mathrm{S}}_{Ip,j}^n{v}_{iI}^{n-1/2}+S_{Ip,i}^{n-1/2}{v}_{jI}^{n-1/2})$$

(6)

$${\mathsf{\Omega }}_{ijp}^{n-1/2}={\displaystyle \sum _I\frac{1}{2}}(S_{Ip,j}^n{v}_{iI}^{n-1/2}-S_{Ip,i}^{n-1/2}{v}_{jI}^{n-1/2})$$

(7)

$${\rho }_p^{n+1}={\rho }_p^n/(1+{\dot{\epsilon }}_{ijp}^{n-1/2}\mathsf{\Delta }t)$$

(8)

where $\mathsf{\Delta }t$ represents the increment of each time step, a constant value in the paper. In this paper, the elastic test stresses in the slope are calculated using a linear elastic constitutive model. Then a return mapping algorithm is used to calibrate the test stresses according to the D-P model and pull back the stresses above the yield surface to the yield surface to obtain the actual stresses.

The yielding criterion in the D-P model illustrated in Figure 2 is defined as:

$$f^s=\tau +q_{\varphi }{\sigma }_m-k_{\varphi }$$

(9)

$$f^t={\sigma }_m-{\sigma }^t$$

(10)

where $\tau =\sqrt{J_2}$ represents the equivalent shear stress, ${\sigma }_m=I_1/3$ is the spherical stress. $q_{\varphi }=\frac{3tan\varphi }{\sqrt{9+12ta{n}^2\varphi }}$, similar to the coefficient of friction, controls the degree of influence of pressure on the yield limit and $\varphi$ is the angle of internal friction of the material. $k_{\varphi }=\frac{3c}{\sqrt{9+12ta{n}^2\varphi }}$ is the yield stress in pure shear and $c$ is the cohesive force of the material.

Then the momentum equation can be solved on the grid nodes based on this information.

Nodal internal forces:

$$f_{iI}^{int}=-{\displaystyle \sum _p{S}_{Ip,j}}{\sigma }_{ijp}\frac{m_p}{{\rho }_p}$$

(11)

where ${\sigma }_{ijp}$ denotes the Cauchy stress of the material points, $S_{Ip,j}$ denotes the derivatives of shape functions, and ${\rho }_p$ denotes the density of the material points.

Nodal external forces:

$$f_{iI}^{ext}={\displaystyle \sum _p{S}_{Ip}}{m}_p{b}_{ip}+{\displaystyle \sum _p{S}_{Ip}}\frac{m_p}{{\rho }_p}{\overline{t}}_{ip}{h}^{-1}\frac{m_p}{{\rho }_p}$$

(12)

where $b_{ip}$ denotes the body force, such as gravity, ${\overline{t}}_i^s={\overline{t}}_i/\rho$ represents the boundary surface force, and $h$ represents the boundary layer thickness.

Total forces:

$$f_{iI}^{}=f_{iI}^{ext}+f_{iI}^{\mathrm{int}}$$

(13)

The continuous time is discretized by the central difference method so that the momentum equation can be expressed as:

$$p_{iI}^{n+1/2}=p_{iI}^{n-1/2}+f_{iI}^n\Delta t$$

(14)

where superscript represents the time step and $p_{iI}^{n+1/2}$ represents the resultant force of node force at the nth time step.

The updated momentum of the nodes will be mapped back to the material points to obtain the momentum of the material points at the next time step. This paper uses the FLIP format for momentum mapping, i.e., mapping the increment of momentum at this time step back to the material points [6,26,27].

$$v_{ip}^{n+1}=v_{ip}+\mathsf{\Delta }t{\displaystyle \sum _I\frac{f_{Ii}{S}_{Ip}}{m_I}}$$

(15)

The positions of the material points can be updated as follows:

$$x_{ip}^{n+1}=x_{ip}^n+\mathsf{\Delta }t{v}_{ip}^{n+1}$$

(16)

where, $x_{ip}^n$ and $x_{ip}^{n+1}$ denote the location of material points at n and (n + 1) time step, respectively.

3. Model

The material point model of the slope is shown in Figure 3, where Figure 3A represents the pure soil slope model, and Figure 3B represents the SRM slope model. The following simplifications have been made when modeling SRM slope, considering that SRM slopes are primarily found in reservoir areas [14], where most of the stones have been washed by river water and are primarily pebbles, thus simplifying the shape of the stones to round [4]. For both slopes, the number of material points is 84,052. For the SRM slope, the number of material points representing stones and soil is 19,077 and 64,975, respectively, and the stone content is 22.7%.

In the simulation, essential boundary conditions are applied on both sides, and the bottom of the model and other surfaces are free. The mechanical parameters of the material are shown in

Table 1. Mechanical parameters of the SRM slope.

Parameters	Soil	Stone
$\rho /\left(\mathrm{kg}/{\mathrm{m}}^3\right)$	1800	2410
$E/\mathrm{MPa}$	50	20,000
$\nu$	0.35	0.2
$\varphi /°$	24	42
$c/\mathrm{kPa}$	10	900
$\psi /°$	24	42

. The parameters are mainly from [9]. Young’s modulus of the rock is 400 times that of the soil, the ratio of the internal friction angle to the dilatancy angle is less than twice, and the ratio of the internal friction force is 90 times. The rock will not reach the yield state in the actual simulation process.

To make the simulation more stable, this paper sets a small-time step, 0.01 ms, in the simulation. The simulation process lasted for 22 s and was divided into two stages. 0~2 s was the first stage, and a stable PIC momentum mapping format [41] and elastic constitutive model were adopted. The gravity is linearly loaded until g = 9.8 m/s², and then the momentum is set to 0 to obtain the initial state of the slope. 2~22 s is the second stage, using a FLIP momentum mapping scheme with good momentum conservation characteristics [42].

4. Simulation Results

4.1. Stability of Slopes

In this paper, the reduced material strength can be expressed as [13]:

$$c^{\prime }=\frac{c}{K_{SRF}}$$

(17)

$${\varphi }^{\prime }=\arctan \frac{\tan \varphi }{K_{SRF}}$$

(18)

where $c$ and $c^{\prime }$ denote cohesion and reduced cohesion, respectively, $\varphi$ and ${\varphi }^{\prime }$ denote internal friction angle and reduced internal friction angle, respectively, and $K_{SRF}$ reparents the reduction factor, which varies from 1 and increases by 0.1 each time until damage occurs on the slope.

We monitored the variation in maximum displacement as the reduction factor increased, as shown in Figure 3.

From Figure 4, we can see that when the reduction factor equals 1.4, the pure soil slope undergoes significant displacement while the SRM slope is still in a stable state. At = 1.7, the displacement of the SRM slope remains small until SFR = 1.8, when the SRM slope undergoes significant displacement. In other words, the coefficients of safety for the pure soil slope and the SRM slope are 1.4 and 1.8, respectively, which proves that the stability of the SRM is better than that of the pure soil slope and that the stones are beneficial for maintaining the stability of the slope in terms of the composition of the slope.

In terms of the tendency of slope displacements, the characteristics of pure soil slopes and SRM slopes are also different. When a pure soil slope is unstable and has large displacements, the displacements always stabilize very quickly, i.e., the movement time is shorter, in most cases reaching the maximum displacement at 5 s, and the displacements show a linear pattern during the growth phase. The SRM slope is characterized by an increase in displacement as the time step progresses and a decrease in the rate of increase, indicating that landslides continue for a longer duration following the instability of the SRM slope. Another distinctive feature is that the displacement curves of pure soil slopes all have a consistent form, consisting of two phases, a linear growth phase, and a horizontal stability phase. The SRM slope displacement curves show two different forms on both sides of the reduction factor equal to 2.0.

In terms of maximum displacement, the maximum displacement of an SRM slope is less than that of a pure soil slope with the same reduction factor, or it is stable without any displacement at all. Even when the SRM slope is given a reduction factor of 2.5, the maximum displacement is less than that of a pure soil slope with a reduction factor equal to 2.0.

Figure 5 plots the plastic zone distribution and displacement distribution clouds for homogeneous soil slopes and SRM slopes at a reduction factor of 1.4 and 1.8, respectively.

According to Figure 5, we can find the following features: The plastic zone distribution of SRM slope and pure soil slope is different at instability. The shear zone of pure soil slope is circular, continuous, and smooth, extending from the foot to the top of the slope. The plastic zone of the SRM slope is interlaced, with one central shear zone running to the top of the slope and the rest of the shear zones interlaced with it, showing obvious rock-winding characteristics. From the displacement cloud map, the displacement distribution of the pure soil slope is also smooth and continuous, while the displacement distribution of the SRM slope shows a layered character.

4.2. Morphology of the Slope Pre and Post Instability

Figure 6 plots the shear zone distribution of the slopes when the reduction factor is close to the safety factor, taking two cases of less than the safety factor and more than the safety factor, respectively.

According to Figure 6, we can find that the plastic zone has penetrated from the foot of the slope to the top of the slope before the reduction factor reaches the safety factor, which is true in both pure soil slopes and SRM slopes. According to Figure 4, we can know that the slope is still in a stable state at this time. For engineering problems where large deformations occur, such as large displacement landslides, the strength reduction method based on the MPM can calculate the whole process entirely from before instability to after instability. Then it is more appropriate for us to judge the safety coefficient of the slope based on whether the displacement changes abruptly or not, compared with the strength reduction method based on the FEM. When the reduction factor is greater than the safety factor, both types of slopes are obviously deformed, and comparing with Figure 5, we can find that the plastic zone is further developed, and the displacement is further increased.

According to Figure 7, we can find that for both pure soil slopes and SRM slopes, the reduction factor is 0.1 larger than the safety factor, and both of them show some similarities and differences, as follows. The variation trend of the two slopes is the same, increase and then decrease. The difference is that the maximum value of the rate of the pure soil slope is larger than that of the SRM slope, which is about four times the maximum value of the rate of the SRM slope. The rate of the pure soil slope decreases rapidly, and at 4 s, its rate decreases by half. In contrast, the rate of the SRM slope only decreases by 20% at this time, which means that the SRM slope has a stronger ability to maintain its motion state after changing from the stationary state to the motion state, which is caused by the large inertia of the SRM.

Another point is that the rate distribution of the two is different. Due to the non-homogeneity of the SRM slope, the rate distribution shows a layered characteristic. The corresponding displacement distribution characteristics are consistent with the rate of both. The displacement growth of homogeneous breakthrough is mainly in the first 4 s, while the displacement of SRM slope grows slowly, and the displacement at the fourth second is only 0.3 m. When the slope regains stability again, its displacement reaches 1.5 m, which is larger than the displacement value of the homogeneous soil slope, about twice as much as the former.

4.3. Destruction of Slopes in Extreme Conditions

The displacement distribution and plastic zone distribution of pure soil slope and SRM slope at extreme conditions are plotted in Figure 8 and Figure 9, respectively.

From Figure 8, we can find that, overall, the response of the pure soil slope is more sensitive to the reduction factors, and its maximum displacement increases by 0.4 m each time the reduction factor increases by 0.1. In the case of the SRM slopes, when the reduction factor increases from 2.3 to 2.4, the maximum displacement increases by 0.2 m, and when the reduction factor increases from 2.4 to 2.5, the maximum displacement does not change. The areas where sliding occurs are identical for both slopes; the difference is that the displacement distribution is discontinuously layered due to rocks in the SRM slope.

From Figure 9, we can find that the distribution of the plastic zone of pure soil slope almost does not change with the increase in reduction factor, both the range of plastic zone and the thickness of sliding zone; the same is not true for SRM slope, the main plastic zone almost does not change with the increase in reduction factor, but one obvious phenomenon is that the plastic zone inside the broken body further develops after the increase in reduction factor, which is due to the fact that the shear capacity of soil is greatly reduced after the increase in reduction factor, and the plastic zone is easily generated in the soil between the stones inside the slope.

5. Conclusions

This paper investigates the stability of SRM slopes and the large deformation process after failure using the material point strength reduction method. In this study, the comparison with pure soil slopes is emphasized, and the exact geometry of isotropic soil slopes as SRM slopes is established. The same numerical simulation is carried out for them, and the following conclusions are obtained.

(1)

The safety coefficient of SRM slopes is higher than that of homogeneous soil slopes, i.e., stones are beneficial to maintain the stability of slopes. In this paper, the safety factor of pure soil slope is 1.4, and the maximum displacement of SRM slope at this discount factor has no sudden shift and is basically equal to 0. Until the reduction factor equals 1.7, the SRM slope has displacement, but the value is small, and the slope is still stable. Finally, the SRM slope is destabilized and destroyed when the reduction factor equals 1.8.

(2)

Before the reduction factor reaches the safety factor of the slope, a plastic zone has been formed in the slope from the bottom to the top of the slope. However, at this time, neither the pure soil slope nor the SRM slope has been damaged. When the reduction factor increases to the safety factor of the slopes or is even more significant than the safety factor, the slopes are immediately damaged, and the pure soil slopes can quickly restore stability after the failure. In contrast, the sliding of the SRM slopes lasts longer.

(3)

In extreme conditions, i.e., when the reduction factors are much more significant than the safety factors of slopes, the maximum displacement of pure soil slopes is more sensitive to this, while the maximum displacement of SRM slopes is less affected by the increase in reduction factors; in addition, the plastic zone of pure soil slopes basically does not change when the reduction factors are significant, either in terms of the extent of the plastic zone or the thickness of the sliding zone, while the plastic zone of SRM slopes develops towards the interior of the slope when the reduction factors are enormous.