Next Article in Journal
Multifactor Mathematical Modeling and Analysis of the Impact of Extreme Climate on Geological Disasters
Previous Article in Journal
Trimethoprim Removal from Aqueous Solutions via Volcanic Ash-Soil Adsorption: Process Modeling and Optimization
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Finite Element Method-Peridynamics Coupled Analysis of Slope Stability Affected by Rainfall Erosion

1
Department of Engineering Mechanics, Hohai University, Nanjing 211100, China
2
Hunan Institute of Water Resources and Hydropower Research, Changsha 410007, China
*
Author to whom correspondence should be addressed.
Water 2024, 16(15), 2210; https://doi.org/10.3390/w16152210
Submission received: 4 June 2024 / Revised: 20 July 2024 / Accepted: 29 July 2024 / Published: 5 August 2024
(This article belongs to the Section Soil and Water)

Abstract

:
Rainfall is a pivotal factor resulting in the cause of slope instability. The traditional finite element method often fails to converge when dealing with the strongly nonlinear fluid–solid coupling problems, making it impossible to fully analyze the sliding process under the state of slope instability. Therefore, this paper uses the coupling of peridynamics (PD) and the finite element method (FEM) to propose a data exchange mode between the seepage field and the deformation field. The influencing factors of fine particle erosion during rainfall are further considered. According to the damage mechanism of the slope sliding process to the original structure of the soil, a modified erosion constitutive relationship is proposed, which takes into account the destructive effect of plastic deformation on coarse particles. Then, the influence of rainfall duration, rainfall intensity, erosion, and initial saturated permeability coefficient on slope stability was simulated and analyzed. This paper provides a novel concept for slope stability analysis and safety evaluation under rainfall conditions.

1. Introduction

A large amount of statistical data show that rainfall, especially heavy rain, is a pivotal cause of geological disasters such as landslides [1,2]. Therefore, slope instability analysis under rainfall conditions is of great significance for construction safety and geological disaster protection. Rainwater soaking reduces soil cohesion, matrix suction, and shear strength, making it more susceptible to damage. Moreover, with rainfall infiltration, there will be erosion and migration of fine particles in granular soils, leading to an increase in porosity and corresponding permeability coefficient, which in turn leads to a further increase in pore water pressure. All of these are the causes of slope instability. Therefore, to accurately simulate the possible instability of soil slopes under rainfall infiltration conditions, it is necessary to reflect the physical coupling of the seepage, erosion, and deformation involved.
Rainfall-induced slope instability is a typical soil–pore fluid coupling problem with strong nonlinearity. The finite element model (FEM) is still the mainstream numerical method for the seepage–stress coupling problem [3,4,5,6]. For the cofferdam of tidal flats behind Changxing Island, Li established its seepage–stress coupling FEM considering the soil strength degradation with the wave cyclic loading, and explored the seepage–stress coupling properties of the cofferdam on a soft clay foundation under a storm surge attack [3]. Xu established a seepage/stress-damage coupling FEM of a deep excavated canal in the Xichuan Section, and explored its long-term performance, including settlement and damage [4]. Zhou used a seepage–stress coupling FEM model to study the rainfall-induced loess landslide including the hazard’s occurrence and evolution [5]. Liu et al. [6] established a finite element numerical model for slope under rainfall load, which realized the unidirectional coupling of the hydro-mechanical behavior of the unsaturated soil slope; then, they studied the influence of rainfall history and saturated permeability coefficient on the internal stress and pore pressure of the slope. The finite element method together with strength reductions is usually used to evaluate slope stability. The slope is judged to be unstable when the slope experiences phenomena such as sudden changes in displacement or non-convergence of the calculation [7]. However, the instability of slopes under rainfall infiltration shows a significant nonlocal landslide morphology, that is, an inconspicuous sliding band, while the results of finite element plastic or damage model simulation show a narrow sliding band morphology, which is not consistent with reality.
When faced with discontinuous media containing crack extension, holes, etc., classical continuum mechanics will fail. To solve this problem, researchers have proposed different calculation methods such as the extended finite element (XFEM) [8,9] and phase field [10,11]. The XFEM avoids the mesh meshing continuously to fit the crack extension by introducing the enriched shape functions reflecting the displacement jump at the crack and the stress singularity at the crack tip into the standard FEM displacement model. In the XFEM, the crack topology is represented implicitly by level sets, which enables the cracks to propagate completely independent of the fixed mesh. However, the large number of enriched nodal degrees of freedom, and complex crack propagation processing involving the determination of the extension direction of crack and iteratively converging the propagation length, greatly increased the computational time. For the phase field method, it is based on the variational approach of fracture mechanics and introduces phase field scalars characterizing the material damage state to simulate discontinuities in a smeared form in narrow banded regions. However, for problems involving multiple physical fields, the fracture phase field method will involve more fields, greatly increasing the difficulty of solving the problem. Different from other methods, Silling proposed peridynamics based on the concept of a nonlocal interaction. It describes material deformation by solving spatial integral equations instead of partial differential equations [12]. It is usually solved using the meshless particle method, which not only has good applicability in simulating damage evolution and crack propagation but also has the potential to calculate problems with spatial characteristic scales or significant nonlocal effects. Because of the above advantages, many scholars have applied peridynamics to the deformation and failure analysis of rock and soil with nonlocal characteristics. Jabakhanji [13] established a peridynamics model for coupling unsaturated soil deformation with transient water–gas flow, and verified the simulation accuracy by comparing it with experimental results. Gu et al. [14] proposed a fluid–structure coupling scheme based on peridynamics to simulate the interaction between soil and pore fluid in the liquefaction analysis of saturated granular soil, and used this method to analyze the liquefaction process of liquefiable soil on the horizontal and inclined ground under seismic action. Zhou et al. [15] established a one-dimensional saturated soil frost heave model based on peridynamics and accelerated the calculation of convolution in the model through a fast Fourier transform. This method greatly improved the calculation efficiency while ensuring the calculation accuracy. Menon and Song [16] established an updated Lagrangian peridynamics model to simulate large deformations of unsaturated soils under drainage conditions. The model introduced the concept of “bond-related secondary neighborhood”, which greatly enhanced the stability and accuracy of the soil skeleton under large deformations. The effectiveness and robustness of the method were verified through numerical examples. Song and Silling [17] used the energy balance of unsaturated porous media to derive the effective stress state and suction state of the porous media skeleton that is conjugated to the nonlocal deformation state. They established the multiphase constitutive correspondence principle between classical unsaturated seepage mechanics and peridynamics through energy equivalence, providing a new way to introduce various constitutive relations into peridynamics. Song and Hossein [18] combined the Cosserat continuum theory with peridynamics and introduced the length scale related to the microstructure as a nonlocal parameter in the model to simulate the shear band bifurcation and cracking phenomena in dry porous media. Menon and Song [19] proposed a nonlocal fluid model for state-based peridynamics, introduced the classical local constitutive model of the solid skeleton and the generalized unsaturated Darcy law into peridynamics, and simulated the shear band of unsaturated soil. Liu et al. [20] proposed an improved bond-based peridynamics moisture–mechanical coupling model to simulate the shrinkage cracking of soil and performed a detailed simulation of the moisture migration and deformation of the soil ring and soil strip. Ren et al. [21] proposed a peridynamics-smoothed particle hydrodynamics coupling strategy to simulate the damage of buried explosions to soil. Liu et al. [22] used fully coupled hygro-mechanical ordinary state-based peridynamics to establish a model of soil strip desiccation deformation and curling, successfully capturing the entire curling process and exploring the influencing factors of the liquid limit, additional evaporation surface, and thickness on the curling performance. Artificial viscosity and virtual particle methods were used in the simulation to improve the calculation accuracy and eliminate the numerical instability caused by shock wave propagation. Sedighi et al. [23] proposed a nonlocal clay erosion formula based on peridynamics, integrating clay swelling, particle separation, and transport into one model. The final simulation results were in good agreement with the experimental data.
In slope instability simulation, even if the slope enters a critical failure state, peridynamics can continue the calculation, which is of great significance for analyzing the failure process of the slope. Lai et al. [24] applied peridynamics to slope stability analysis and used non-ordinary state-based peridynamics combined with the Drucker–Prager constitutive model to model a failure process of a two-dimensional soil slope. The results were consistent with the finite element method and the dynamic sliding process of slope failure was given. Zhang and Zhang [25] used ordinary state-based peridynamics combined with the Drucker–Prager yield criterion to successfully simulate the localized deformation of the slope and locate the critical failure surface of the slope. This method successfully avoided the phenomenon of non-physical plastic property changes under extreme non-uniform deformation. Wang et al. [26] combined non-ordinary state-based peridynamics with random field models to study the influence of spatial variability of soil parameters on slope strength. They also compared the effects of pulse and non-pulse earthquakes on the landslide process. Zhou et al. [27] used peridynamics to research the relationship between shear bandwidth and horizon size, and combined non-ordinary state-based peridynamics with the strength reduction method to conduct safety analysis on intact slopes and slopes with defects, respectively, providing guidance for real engineering. Although peridynamics has many advantages over the traditional finite element method, since it considers nonlocal effects, it requires more computing resources than traditional methods during discretization, which greatly prolongs the calculation time. It also limits the application of peridynamics in practical problems. To solve this problem, many scholars have conducted research. Yang and Liu [28] established a multiscale method based on the combination of a boundary element method and peridynamics. By introducing boundary element nodes at the same position as the PD material points on the interface, a shared node coupling model was established to give full play to the advantages of both methods. While accurately simulating object deformation and crack propagation, higher computational efficiency was achieved. Ni et al. [29] used PD and FEM coupling to simulate the propagation of hydraulic fractures in saturated porous media. Liu et al. [30] proposed an element-based coupling method of peridynamics and the finite element method. Sun and Fish [31] used non-ordinary state-based peridynamics coupled with the finite element method combined with Boit theory to simulate crack propagation in saturated porous media. Jin et al. [32] established a coupling model of non-ordinary state-based peridynamics and the finite element method, connecting the two methods by embedding peridynamics material points in the interface elements, and proposed two coupling force distribution schemes: volume coupling and contact coupling. Although some researchers have incorporated peridynamics into the slope stability analysis, the integration of erosion effects with peridynamics remains an area underexplored. Consequently, the simulation of slope instability under the influence of rainfall could benefit from further refinement.
Fine particle erosion is one of the important factors affecting slope instability. The basic mathematical model of erosion is mainly developed based on experimental simulation and porous continuum mechanics. Dahaghi et al. [33] established an expression for the permeability during the erosion process based on the mass conservation equation for particle transport in porous media, derived an analytical model for particle movement and deposition, and considered the effects of the inflow fluid rate and the concentration of particles it carries on the erosion process. Yang et al. [34] divided the soil into four components: a stable skeleton composed of coarse particles, fine particles that can be eroded, liquefied particles that have been eroded, and seepage liquid. They then established a mass conservation equation to describe the migration of liquefied particles. Golay et al. [35] studied the erosion between fluid and soil particle surfaces at the pore scale and described the evolution of the water–soil interface using a level set function. Chang [36] conducted experiments to study the corrosion process under complex conditions, analyzed the influence of hydraulic gradient and stress state on erosion, and derived the control equation of the critical failure hydraulic gradient. Sterpi [37] experimentally studied the erosion of fine particles in soil samples by controlled seepage and established an FEM model to analyze the effect of fine particle loss on the stress–strain distribution in the soil. Zhang et al. [38] established an unsaturated seepage and erosion coupled model through a study of the interaction between erosion and infiltration and analyzed the impact of erosion on soil slope stability.
Considering the above statement, this paper attempts to propose a PD-FEM coupling model and a modified erosion constitutive model to simulate and analyze the slope instability caused by heavy rainfall, and explore the influence of rainfall duration, rainfall intensity, erosion amount, and initial saturated permeability coefficient on slope safety. Specifically, Section 2 gives the basic formula of non-ordinary state-based peridynamics and establishes a numerical model of the slope. Section 3 proposes a coupling scheme of non-ordinary state-based peridynamics and finite element method, and explores the effects of rainfall duration and rainfall intensity on slope stability. Section 4 further considers the erosion of soil, proposes a modified erosion constitutive model, and explores the effects of erosion and initial saturated permeability on slope stability. Finally, Section 5 summarizes the main conclusions of this paper.

2. Elastoplastic Analysis Based on Non-Ordinary State-Based Peridynamics

2.1. Non-Ordinary State-Based Peridynamics

Peridynamics analyzes the deformation and motion process of the object by studying the interaction of a material point with other points within its horizon. For particles beyond the horizon, peridynamics assumes that they have no interaction, as shown in Figure 1. x and x are the two material points interacting in the initial configuration within the horizon, and y and y are the material points interacting within the original horizon in the deformed configuration. From this, the initial deformation vector state and the current deformation vector state can be defined X _ ξ = x x , Y _ ξ = y y . Among them, ξ represents the relative position of the material point in the initial configuration, and   represents the vector corresponding to the state field.
Deduced from the principle of virtual work, the motion equations of peridynamics can be obtained as
ρ ( x ) u ¨ ( x , t ) = H { T _ [ x , t ] x x T _ [ x , t ] x x } d V x + b ( x , t )
where ρ and u ¨ represent the material density and acceleration, respectively; T _ is the force vector state, which represents the force interaction between the material point and points within its horizon. Let H { T _ [ x , t ] x x T _ [ x , t ] x x } d V x be denoted as L ( x , t ) , which represents the total force density at point x ; d V x is the differential volume of material point x ; b ( x , t ) represents the body force density; δ is the size of the horizon, and its value is related to the characteristic length of the material and the phenomenon to be studied. Using H = H ( x , δ ) : = { x R : { x x δ } } can represent all material points within the horizon of x .
Finally, the calculation formula for the force state can be obtained as follows:
T _ = ω _ ( | ξ | ) P K 1 ξ
At the same time, a Gaussian influence function as shown in Formula (3) is used in the force vector state to reduce the calculation error of the deformation gradient tensor and force density, effectively improve the accuracy of applying boundary conditions, and improve calculation stability,
ω _ ( | ξ | ) = e | ξ | 2 δ 2
In the theory of non-ordinary state-based peridynamics (NOSBPD), the nonlocal shape tensor can be defined using the union of the relative position vectors of the material points in the initial configuration:
K [ x , t ] = H ω _ ( | ξ | ) ( ξ ξ ) d V x
where represents the dyadic product of two vectors. ω _ ( | ξ | ) is the influence function, and their interaction size can be controlled according to the distance between material points.
Classical continuum mechanics defines the deformation gradient tensor by establishing the relationship between the initial configuration and the current configuration, i.e., F = y ( x , t ) x . Similarly, peridynamics uses the nonlocal shape tensor to obtain the nonlocal deformation gradient tensor.
F [ x , t ] = H ω _ ( Y _ ξ ξ ) d V x K 1
Then, represent the first Piola–Kirchhoff stress using the deformation gradient
P = det ( F ) σ F T
where P is the first Piola–Kirchhoff stress, det ( ) is the determinant of the matrix, and σ is the real Cauchy stress in the current configuration. It can be known from Formula (7):
σ = R τ R T
where τ is the irrotational Cauchy stress, which can be obtained through the stress–strain relationship, R is the orthogonal rotation tensor.
According to continuum mechanics, the Lagrangian-Green strain tensor can be expressed by the deformation gradient as E = 1 2 ( F T F I ) , and for small deformation problems, it can be simplified to E = ( F + F T ) 2 I .
To suppress the numerical oscillation existing in NOSBPD, a penalty force density f hg [39] is added to the force vector state, which is defined as follows:
f hg = C hg C modulus h proj | x j x i | y j y i | y j y i |
where C hg is the hourglass constant coefficient, C modulus is the micromodule of the prototype microelastic brittle model, h proj is the projection of the hourglass vector h on the vector x x . The expression of the hourglass vector is h = y + F ( x x ) y .
Adding an artificial damping term C to the original peridynamics motion equation and specifying an appropriate damping value can quickly obtain a quasi-static solution to the equation when solving steady-state problems, improving the calculation speed. By discretizing the modified motion equation, we can obtain the following:
ρ i u ¨ i + C u ˙ i = j = 1 N f hg V j + b ( x i ) + j = 1 N [ T _ x j x i T _ x i x j ] V j
where N is the total number of material points in the horizon of material point i , and V j is the volume of material points. When the displacement and velocity of the material point at the current time step are known, the acceleration can be obtained by solving the equation. Use the Velocity–Verlet integration method to perform explicit integration, calculate the velocity and position of the material point at the next time step, and complete the state update of the deformed body.
u ˙ m i d = u ˙ i n + Δ t 2 ρ ( L + b ) n u ˙ i n + 1 = u ˙ m i d + Δ t 2 ρ ( L + b ) n + 1 u i n + 1 = u i n + u ˙ i n Δ t + ( Δ t ) 2 2 ρ ( L + b ) n + 1

2.2. Elastoplastic Constitutive Model Updating

The irrotational strain rate tensor d can be decomposed into the spherical strain rate ε ˙ m , the elastic deviator strain rate e ˙ e and the plastic deviator strain rate e ˙ p :
d = ε ˙ m + e ˙ e + e ˙ p
where the spherical strain rate elasticity ε ˙ m = 1 3 tr ( d ) I , the elastic deviator strain rate can be obtained from Hooke’s law e ˙ e = S ˙ 2 G , S ˙ is the deviator stress rate, and G is the shear modulus.
The soil constitutive uses the perfect elastoplastic constitutive with the Drucker–Prager yield criterion. When using the associated flow law, the yield function can be written as
f = J 2 + A I 1 B
where I 1 is the first invariant of stress, J 2 is the second invariant of deviator stress, A and B are the parameters that control the yield criterion. For the two-dimensional plane strain problem, the calculation formulas of the two are A = tan φ 9 + 12 tan 2 φ and B = 3 c 9 + 12 tan 2 φ . φ and c are the friction angle and cohesion of the soil, respectively. According to the calculated yield function, when f < 0 , the soil is in the elastic state, and when f 0 , the soil enters the plastic state.
The direction of the plastic deviatoric strain rate should satisfy the orthogonal flow law. Due to the use of the associated flow law, the plastic potential function is equal to the yield function, i.e., g = f . The expression of the plastic strain rate is as follows:
e ˙ p = λ ˙ g τ         = λ ˙ ( S 2 | S | + A I )
The plastic multiplier λ ˙ can be calculated from the consistency condition:
λ ˙ = tr ( T ) tr ( T )
In the formula, = 2 S | S | + 2 A I , = 2 d + ( 3 K G 2 ) ε ˙ .
After obtaining the plastic multiplier, the plastic deviator strain rate can be obtained by bringing back Equation (13). Then, for the perfect elastoplastic model, the relationship between the stress rate is
τ ˙ = 3 K ε ˙ m + 2 G ( d ε ˙ m λ ˙ ( S 2 | S | + A I ) )
Update the stress at the current time step by the stress rate
τ n + 1 = τ n + Δ t τ ˙ n + 1

2.3. Case Verification and Analysis

A slope as shown in Figure 2 is established, and the instability process of the slope under the action of gravity is analyzed. The slope soil has a mass of 20 kN/m3, an elastic modulus of 100 Mpa, and a Poisson’s ratio of 0.3. The Drucker–Prager yield criterion and associated flow law are used to describe the plastic characteristics of the soil. The internal friction angle is 21 and the cohesion is 13.03 kPa. The AB and CD edges constrain displacement in the x direction, the BC edge constrains displacement in both the x and y directions, and the other boundaries are unconstrained.
As shown in Figure 3, when the slope enters a stable state under the action of gravity, the equivalent plastic strain at the slope corner is the largest and the soil damage is the most serious. The plastic zone tends to penetrate from the foot of the slope along the arc to the top of the slope, which is consistent with the shape and damage trend of the slope slip zone observed in Figure 4. It also reflects the nonlocal phenomenon in the real landslide process and verifies the reliability of the model. Although the slope is still in a safe state at this time, if there are external factors such as heavy rainfall, landslide damage is very likely to occur, and the stability of the slope under extreme conditions must be further analyzed.

3. Slope Stability Analysis under Rainfall Conditions with NOSBPD-FEM Modeling

3.1. Fluid–Structure Interaction Solution

Unsaturated soil is composed of three phases: soil particles, water, and air. Its physical properties will change significantly when affected by the complex external environment. Among them, water content has the strongest impact on unsaturated soil. During rainfall, rainwater falls on the slope surface and is absorbed by the surface soil particles under the action of matrix suction. When the absorbed water reaches an extreme value, the microcracks in the soil will become infiltration channels for water, resulting in unsaturated seepage, which will increase the volume moisture content of the lower soil layer and cause the pore water pressure u w of the soil to increase. Since the pores in the soil are assumed to be always connected to the outside atmosphere, the pore air pressure u a will not change significantly, resulting in a continuous decrease in the matrix suction ( u a u w ) , which in turn affects physical parameters such as soil permeability coefficient, cohesion, and internal friction angle [40].
According to the Terzaghi principle of soil mechanics, saturated soil can be regarded as consisting of two phases: soil particle skeleton and pore water. When the soil pores are filled with water, the pore water pressure will bear part of the stress, while the deformation of the soil is mainly controlled by the effective stress of the soil skeleton, and the pore water pressure does not contribute to the strength and deformation of the soil. The situation of unsaturated soil is more complicated. It is composed of three phases: soil particle skeleton, pore water, and gas in the pores. The pore gas pressure also plays a role in sharing some of the stress. To describe the effective stress of unsaturated soil, Zhao et al. [41] proposed a specific expression for the effective stress of unsaturated soil based on the work-based effective stress principle of unsaturated soil proposed by Dean and Houlsby [42]:
σ = σ [ S r u w + ( 1 S r ) u a ] I
where σ is the effective stress; σ is the total stress; S r is the saturation; and I is the second-order unit tensor.
To consider the effect of seepage on the original stress state of the soil, effective stress needs to be used to correct the peridynamic equation of motion. Decompose the Cauchy stress in Equation (6) into two parts: spherical stress and deviatoric stress, and obtain the Cauchy spherical stress term reflecting the pore water pressure, and obtain the corrected Cauchy stress. Use the corrected Cauchy stress to continue calculating the force density state.
Soil is a loose structure with many pores inside. When rainwater from outside flows in and displaces the original air in the pores, the density of the soil will increase, increasing the sliding force of the slope and harming an adverse effect on the stability of the slope. The updated expression of soil density is
ρ = ρ 0 + θ w ρ w
where ρ 0 is the dry soil density; θ w is the current volume moisture content of the soil; ρ w is the density of water.
During rainfall, the infiltration water flow will also drag the soil particles, generating a seepage volume force f s in the same direction as the seepage velocity, and its magnitude is proportional to the water head gradient.
[ f s ] = [ γ H x γ H y ]
where γ is the soil specific weight; H is the total water head, which can be obtained by solving the seepage control equation [43]
x ( K x ( θ ) H x ) + y ( K y ( θ ) H y ) = θ t
where K x and K y are the permeability coefficients in the direction x and the direction y respectively; θ is the volume water content.
This paper selects the Van-Genuchten water–soil characteristic curve equation and the unsaturated permeability coefficient equation to describe the seepage characteristics of the soil. The Van-Genuchten water–soil characteristic curve expression is [44]
θ w = { ( θ s θ r ) { 1 + [ α w ( u a u w ) ] n w } m w + θ r ( u w < 0 ) θ s ( u w 0 )
where θ s is the saturated volume moisture content; θ r is the residual volume moisture content; u a is the pore gas pressure; α w is a parameter related to the inverse of the air intake value; n w is a parameter related to the slope of the water–soil characteristic curve; m w is a parameter related to the residual state.
The unsaturated permeability coefficient equation is
K = K s { 1 ( α w u w ) n w 1 [ 1 + ( α w u w ) n w ] m w } 2 [ 1 + ( α w u w ) n w ] m w 2
where K s is the saturated permeability coefficient.
After substituting the seepage volume force from Equation (19), it can be added to the peridynamic equation of motion as an external body force density term.
Soil deformation can also affect seepage. After calculating the volume strain of the soil through peridynamics, the porosity after deformation can be updated according to soil mechanics as
n = n 0 ε v ( 1 + n 0 )
where n 0 is the initial porosity of the soil and ε v is the volume strain. The change in porosity leads to a decrease in the cross-sectional area available for water flow, and the saturated permeability of the soil decreases. The relationship between the saturated permeability and porosity is as follows:
K s = n 3 ( 1 n ) 2 ( 1 n 0 ) 2 n 0 3 K s 0
where K s 0 is the initial saturated permeability coefficient. The change in saturated permeability coefficient will cause the change in soil permeability function, thus affecting the seepage distribution of slope.
In summary, after establishing the fluid domain and solid domain models, the pore water pressure distribution is obtained according to the initial and boundary conditions of the seepage, and the corresponding seepage volume force and density change are calculated. Considering the effect of the pore water pressure on the solid domain that causes the effective stress to decrease, the slope displacement and stress under the influence of seepage are calculated. The volume strain generated in the solid domain will cause the porosity to change, changing the permeability coefficient of the fluid domain, thereby affecting the pore water pressure distribution of the seepage field in the next time step, thus achieving the coupling between the solid domain and the fluid domain. Compared with other researchers on fluid–structure coupling [45,46,47], although the coupling model in this paper adopts a nonlocal method, its purpose is to capture the nonlocal characteristics of the slope slip zone and does not consider the influence of crack damage. At the same time, the soil uses a perfect elastoplastic model, and the influence of parameters such as soil saturation is not reflected in the constitutive model, which greatly simplifies the soil model.
In the coupling process, since peridynamics uses uniform meshfree particle discretization, while the finite element method uses four-node quadrilateral isoparametric element discretization, the data storage points of the two do not correspond one to one. Usually, the distribution of peridynamics material points is denser than that of FEM nodes. Therefore, the element number and local coordinates of each material point within the element, as well as the number of material points near each node, should be calculated in advance. When the coupled calculation is started, the required data for the peridynamics material points can be obtained from the element node data using shape function interpolation. Since the distribution of material points is denser, each FEM node can be considered to be located in the area formed by the nearest four PD material points. If this area is regarded as an element, the data required by the node can still be obtained by interpolating the data of nearby material points and shape functions. Thus, a data exchange model between the deformation field and the seepage field is established.
In terms of program architecture, since the peridynamics part uses explicit iteration to solve the motion equations, the maximum time step that satisfies the iterative convergence is much smaller than the timescale of the rainfall process. It is unrealistic to calculate the entire rainfall process with the peridynamics time step. Therefore, when writing a fluid–solid coupling program, the time step taken by the finite element method is used as the real time. The peridynamics only imports the pore water pressure and other data and calculates them until convergence after the FEM module calculates one or several time steps. The program flowchart Figure 5 is as follows:

3.2. Transient Flow in Unsaturated Soil Pillar

To verify the accuracy of the finite element method for simulating unsaturated and unstable seepage problems, the soil column model [48] shown in Figure 6 was established. The soil column height H = 1.0 m, width W = 1.25 × 10 2   m , saturated permeability coefficient K s = 1.0 × 10 6   m / s . The parameters of the VG model are θ s = 0.363 , θ r = 0.186 , a = 1.0 m−1, and n = 1.53. The initial pressure head of the soil column is h = 8.0   m .
The soil column FEM model uses a four-node quadrilateral element, which is evenly divided into 500 elements, with a time step of Δ t = 0.1 s , and a calculation time of 468,000 time steps. The AB and CD sides of the soil column are constant head boundaries, with head values of 0 m and −8 m respectively, and the AD and BC sides are impermeable boundary conditions.
Figure 7 shows the comparison results of the FEM solution and Warrick’s analytical solution [49] at different instants. The results show that as the seepage progresses, the upper part of the soil column gradually becomes saturated, and the pressure head below the saturated zone decreases sharply, while the pressure head below the soil column still maintains the initial state and is not affected by the seepage. The simulation results of the finite element can correspond well with the analytical solution, verifying the accuracy of the model.

3.3. Numerical Validation of Fluid–Structure Interaction Scheme

To verify the feasibility and accuracy of the above fluid–solid coupling model, a soil column model with a height of 1m and a width of 1.25 × 10 2   m was established. Its mass density is 1.72 × 10 3   kg / m 3 , the permeability coefficient is 1.75 × 10 7   m / s , and the compression modulus is 3.4 Mpa. The sides and bottom of the soil column are impermeable, and displacement constraints in the x direction are applied to the sides of the column, and displacement constraints in the x and y directions are applied to the bottom of the column. Water seepage and soil compression only occur in the vertical direction. The top of the soil column is drained on one side and is subjected to a uniformly distributed downward load of 100 kPa.
Reference [50] gives an analytical solution for the pore water pressure at any depth z of a soil column at different times t.
p = 4 π P 0 m = 1 1 m ( sin m π z 2 H ) exp ( m 2 π 2 4 T ν )
T v = C v t H 2 C v = E s K γ w
where H is the depth of soil column, m is a positive odd number, C v is the consolidation coefficient of soil, γ w is the density of water, T v is the dimensionless time factor, E s is the compression modulus, P 0 is the uniformly distributed load applied on the top.
Figure 8 shows a comparison between the analytical solution and the numerical solution of the pore water pressure distribution at t = 100 s. The pore water pressure in the lower section of the soil column is linearly related to the depth, while the pore water pressure in the upper section decreases slowly. It can be seen from the images that the numerical solution is in good agreement with the analytical solution, verifying the accuracy of the coupling model.
Figure 9 is a comparison of the Abaqus finite element solution and the homemade finite element solution of the vertical displacement distribution at t = 100 s. Since the pore water pressure bears part of the external load, the vertical displacement of the soil column regarded as an elastic body is not linearly distributed along the height, and the maximum settlement value of the soil column is reduced. It can be seen from the image that the development trends of the two curves are the same, which further verifies the accuracy of the coupling model.

3.4. Effect of Rainfall Duration on Slope Stability

Select the slope model in the example in Section 2 for further simulations. The two side boundaries AB and CD of the slope constrain displacement in the x direction; the bottom BC of the slope constrains displacement in the x and y directions; AB and BC are impermeable boundaries. The constant head boundary conditions are applied on the CD and DE sides, where the pressure head on the CD side is linearly distributed h c = ( 3 Y )   m and the pressure head on the DE side is always 0. Flow boundary conditions are applied on the AF and FE edges to simulate rainfall infiltration. Heavy rainfall often forms runoff and water accumulation on the slope surface. Since the depth of water accumulation is difficult to determine, this paper assumes that there is no surface water accumulation during the rainfall process, and the rainwater that does not infiltrate will be drained away immediately. The initial saturation distribution of the slope is shown in Figure 10. The groundwater level is 3.25 m high and parallel to the bottom of the slope.
If the flow boundary condition is simply applied, pressure infiltration will occur on the slope surface, which is unreasonable if there is no surface water accumulation. Therefore, this paper reduces the flow rate by adding a proportional coefficient α to avoid pressure infiltration, and the reduced flow rate is regarded as discharge from the slope. The initial value of the proportional coefficient is 1. When the pressure head of the flow node is greater than zero, let α = α β , β be the single reduction size, and it can be taken as 0.5 for the first reduction. If the pressure head is still greater than zero after the reduction, it will continue to be reduced with the current reduction size until the pressure head is less than zero. At this time, it is determined that the true proportional coefficient exists in the interval composed of the coefficients before and after the reduction. The average value of the coefficients before and after the reduction is used as the search result for this time and substituted into the trial pressure head. If the accuracy requirement is met, it will be used as the flow boundary condition to continue the seepage calculation; otherwise, the reduction coefficient will be halved and the binary search will continue.
To observe the effect of rainfall duration on slope stability, the pressure head distribution cloud diagrams before rainfall, 6 h, 12 h, and 24 h of rainfall are given when the rainfall intensity is 5 mm/h, namely Figure 11. When there is no rainfall, the pressure head is distributed linearly with height, decreasing from 3 m to −10 m from bottom to top. As the rainfall progresses, the top and surface of the slope quickly reach saturation 6 h after rainfall, the pressure head becomes zero, a transient saturated zone is formed, and a wetting front is generated at the junction with the unsaturated zone. The maximum negative pressure head appears below the wetting front. Since the rainfall boundary at the foot of the slope is close to the groundwater level, rainwater can quickly flow into the groundwater, causing the water level below it to rise slightly; while at the top of the slope, the seepage distance is long and there is a moist front, and there is still a large unsaturated area below, so the corresponding groundwater level does not change significantly.
When the rainfall reaches 12 h, the wetting front descends, the slope matrix suction decreases, the groundwater level rises further, and the maximum pressure head appears on the right side of the slope bottom. After 24 h of rainfall, the saturated zone above the moist front further developed, wrapping the unsaturated zone in the middle of the slope. The groundwater level increased significantly and gradually developed into the interior of the slope with the continuous influence of seepage.
The equivalent plastic strain distribution before rainfall, 6 h of rainfall, 12 h of rainfall, and 24 h of rainfall is shown in Figure 12. The equivalent plastic strain is the largest at the slope toe. As rainfall continues, a transient saturated zone is formed at the top and below the slope, weakening the total stress of the soil. At the same time, the dragging effect of the volume force generated by rainwater seepage on soil particles will also harm the slope stability. Since the density of rainwater is much greater than that of air, the weight of soil in the transient saturated zone increases significantly and the sliding force increases. Under the influence of the above factors, the high plastic strain zone at the toe of the slope continues to develop forward, the plastic strain zone gradually penetrates and expands, and eventually leads to instability of the slope along the slip zone.

3.5. Effect of Rainfall Intensity on Slope Stability

To observe the effect of rainfall intensity on slope stability, the slope pressure head and equivalent plastic strain after 24 h of rainfall with intensities of 1 mm/h, 2.5 mm/h, 5 mm/h, and 7.5 mm/h were compared. The pressure head distribution is shown in Figure 13. As the rainfall intensity increases, the transient saturated zone will develop deeper into the slope, the pressure head at the wetting front will increase significantly, and the matrix suction will dissipate faster. When the rain intensity is 1 mm/h, there is no significant change in the pressure head at the foot of the slope, and when the rainfall intensity increases to 5 mm/h, the pressure head at the slope foot increases significantly, weakening the total stress of the soil, making the slope foot more susceptible to damage. There is no obvious change in the pressure head for rainfall intensities of 7.5 mm/h and 5 mm/h. This is because the rainfall surface permeability has reached its limit. Although the rainfall intensity of 7.5 mm/h is higher, the flow rate higher than the surface permeability does not participate in the seepage process.
The equivalent plastic strain change is shown in Figure 14. As rainfall intensity increases, the equivalent plastic strain zone gradually penetrates to the top of the slope. However, no matter how the rainfall intensity changes, the plastic strain at the bottom of the slope is always much larger than that at the top of the slope. On the one hand, this is because the bottom of the slope is subjected to a large gravity load from the soil above, and there is a geometric mutation at the toe of the slope that leads to stress concentration. On the other hand, since the ED side is a constant head boundary condition and the EF side connected to it is a flow boundary condition, part of the rainwater will flow out from the ED side through the toe of the slope as soon as it infiltrates into the slope, resulting in a large head gradient here. The resulting penetration force drags the soil particles, which also increases the plastic strain. At the same time, the groundwater level makes the pressure head at the bottom of the slope much greater than that at the top of the slope, further weakening the effective stress of the soil and harming the stability of the slope.
The equivalent plastic strain at the slope foot changes with time under different rainfall intensities, as shown in Figure 15. The growth rate of equivalent plastic strain was high before 5 h, and then it gradually slowed down and tended to be stable. The reason is that 5 h ago, a part of the soil near the slope toe was still unsaturated. As the rainfall progressed, the weight of this part of the soil increased, and the pore water pressure increased, which jointly affected the development of equivalent plastic strain. When all the nearby soils reached saturation, the weight tended to be stable. At this time, only the pore water pressure increased with time, resulting in a slowdown in the growth rate.

4. NOSBPD-FEM Analysis of Slope Stability under Coupling Effect of Erosion and Seepage

4.1. Seepage–Erosion Coupling Theory

The soil skeleton can be divided into coarse particles and fine particles according to the size of the particles that make up it. The coarse particles contact each other to form the basic structure of the soil skeleton, and the fine particles fill the pores between the coarse particles. Due to the low degree of consolidation of fine particles, under the action of seepage, some fine particles will migrate with the water flow to form liquefied particles. These fine particles can be transported through the seepage channel and redeposited, while coarse particles are relatively stable and will not be eroded by water flow. This physical process is called erosion. If a large amount of fine particles filling the pores are lost, the basic structure composed of coarse particles will become unstable. Therefore, erosion poses a serious threat to slope stability.
The process of fine particle erosion from initiation to development is shown in Figure 16. When erosion just starts, there are only liquefied particles that migrate from other locations to the microelement in the pores. As erosion further develops, the fine particles originally attached to the coarse particles in the microelement detach from the microelement to form liquefied particles. Therefore, the reasons for the change in the density of liquefied particles can be divided into three parts: the liquefied particles are produced by the stripping of fine particles from the soil skeleton; some liquefied particles are redeposited into the soil skeleton; and the inflow and outflow of liquefied particles. Based on the above analysis, the mass conservation equation of liquefied particles in the microelement can be written as
ρ tr t + ( ρ tr v x ) x + ( ρ tr v y ) y = q er q dp
where ρ tr is the density of liquefied particles, where the density is obtained by dividing the mass of liquefied particles by the unit soil volume; v is the flow velocity of liquefied particles; q er and q dp are the masses of liquefied particles produced by erosion in a unit volume of soil per unit time and the masses of liquefied particles deposited on the soil skeleton, respectively.
Since the liquefied particles are carried by the water flow, it is assumed that the velocity of the liquefied particles is consistent with the seepage velocity of the water. Substituting the Darcy seepage velocity q into the mass conservation equation, we can obtain
ρ tr t + x ( ρ tr q x n ) + y ( ρ tr q y n ) = q er q dp
where n is the soil porosity.
The erosion rate of fine particles in the soil skeleton is related to the shear force exerted by rainwater seepage, the content of fine particles that can be eroded, and the erosion resistance of the soil. This paper improves the soil erosion constitutive model proposed by Cividini et al. [51] to simulate the erosion process. The liquefied particle rate generated by erosion can be obtained from Formula (24):
q er ( t , v ) = β er | v | [ ρ f ( t ) ρ ( v ) ]
where β er is a non-negative material constant, reflecting the sensitivity of the soil to the infiltration flow; | v | is the real flow rate of water; ρ f is the density of fine particles in the soil skeleton that can be eroded at the current moment, and it can be expressed by the density of fine particles in the soil skeleton that can be eroded in the initial state and the density of liquefied particles in the soil: ρ f = ρ f 0 ρ tr ; ρ is the steady-state density of fine particles in the soil skeleton at the current flow rate, which can be calculated by the Formula (25):
ρ ( v ) = { ρ f 0 ( ρ f 0 ρ f * ) v v * ( 0 v ( t ) v * ) ρ f * α er lg ( v v * ) ( v * v ( t ) )
where v * is the critical velocity of pore water when fine particles are eroded; ρ f * is the steady-state density of fine particles in the soil skeleton at the critical velocity; α er is the parameter that controls the rate of change in the steady-state density of fine particles with flow velocity. Since the original erosion constitutive model does not consider soil deformation, ρ f is only affected by the liquefaction particle rate during the erosion process. When the slope becomes unstable due to external factors, a slip zone will be generated, which will change the original structure of the soil and cause some coarse particles to be destroyed into fine particles. Therefore, it is not completely correct to use the original erosion constitutive model. Therefore, this paper makes corrections to the erosion model.
The failure of coarse particles mainly comes from shearing action, so the equivalent plastic strain is used as the criterion for the degree of failure of coarse particles in this paper. Although the slope soil studied in this paper appears to be homogeneous at the macro level, the skeleton composed of coarse particles at the micro level still has a random distribution phenomenon. Under the same plastic strain, whether the coarse particles at a certain location are destroyed is related to the geometric structure and can be regarded as a probabilistic event. This paper uses the Gaussian distribution probability model to simulate the relationship between equivalent plastic strain and coarse particle failure.
ρ f = ρ f 0 + Φ ( ρ fs ρ f 0 ) ρ tr
Φ ( x ) = 1 2 ( 1 + erf ( x μ σ 2 ) )
where ρ fs is the maximum density that fine particles in the soil can reach; Φ is the cumulative distribution function of the Gaussian distribution; erf ( · ) is the error function; μ and σ are the mathematical expectation and variance, respectively. The liquefied particle deposition rate q dp is not considered in this paper.
Rainwater seepage, erosion, and slope deformation will affect the soil porosity. According to the definition of porosity, the porosity at the initial moment and the next moment are
n 0 = V b 0 V s 0 V b 0
n = ( V b 0 + Δ V b ) ( V s 0 + Δ V s ) V b 0 + Δ V b
where V b 0 is the initial volume of soil; V s 0 is the initial volume of soil skeleton; Δ V b is the change in soil volume; Δ V s is the change in soil skeleton volume. Since water displaces the air in the pores during rainfall and does not affect the total volume of the soil, the change in soil volume is mainly caused by the deformation of the soil skeleton. The volume change in the soil skeleton is composed of the volume lost due to the liquefaction of fine particles caused by erosion of fine particles. The volume strain ε v induced by the soil skeleton can be calculated by peridynamics, and the volume changes of soil and soil skeleton and can be written as
Δ V b = V b 0 ε ν
Δ V s = q er Δ t V b 0 ρ s
where ρ s is the soil particle density; Δ t is the time step. Substituting Equations (35) and (36) into Equation (34), we can obtain
n = n 0 + ε v + ς 1 + ε v
where ς = q er Δ t ρ s .

4.2. Stability Analysis of Rainfall Slope Considering Erosion

The slope model remains consistent with the example in Section 2. The initial density of the soil is 2 g/cm3, the mathematical expectation and variance of coarse particle failure are and respectively, and the maximum density that fine particles can reach is 0.8 g/cm3, the remaining parameters are given in Table 1. Taking the rainfall intensity as 10 mm/h, the saturated permeability coefficient of the example in this chapter is relatively large, and there will be no runoff or water accumulation. The rainfall intensity does not need to be reduced, and all rainfall will immediately penetrate into the slope.
The trend of liquefied particle density over time is shown in Figure 17. It can be seen that as the rainfall continues, erosion gradually develops inward along the top and surface of the slope, among which the erosion develops most rapidly at the slope foot. The reason is that the CD and DE edges are constant head boundary conditions. During heavy rainfall, the infiltrated rainwater is continuously discharged from the slope, resulting in a large head gradient, high seepage velocity, and fast liquefied particle generation rate. At the same time, rainwater also transports the liquefied particles generated elsewhere to the slope corners, resulting in a higher density of liquefied particles here.
By observing Figure 18, it can be found that the density of liquefied particles in the downslope increases with the slope height at different rainfall durations. Taking the AB side on the left side of the slope as the reference line, the distribution curve of liquefied particle density along the depth under different rainfall durations is drawn, as shown in Figure 16. The density of liquefied particles is distributed linearly along the height. As the rainfall duration increases, the curve begins to bend, the liquefied particle generation rate at the top gradually exceeds that at the bottom, and the erosion speed accelerates.
The change in equivalent plastic strain is shown in Figure 19. As the rainfall continues, the equivalent plastic strain zone gradually expands, and the plastic strain at the slope foot tends to develop downward. This is closely related to the greater erosion at the toe of the slope. Subsidence erosion increases the porosity of the area, affects the permeability coefficient, allows rainfall to fall to the slope foot faster, and weakens the total stress of the soil. Compared with the case when erosion is not considered, the loading history at the slope toe changes, resulting in a different trend in the development of plastic strain.

4.3. Effect of Initial Saturated Permeability Coefficient on Slope Stability

The liquefied particle density distribution of soil slopes with different initial saturated permeability coefficients after 24 h of rainfall is shown in Figure 20. The initial saturated permeability coefficient has a great influence on the generation rate of liquefied particles. When the initial saturated permeability coefficient increases from 5 × 10−5 m/s to 1.5 × 10−4 m/s, the peak density of liquefied particles at the slope angle increases by 2.5 times. A larger initial saturated permeability coefficient will make rainwater infiltrate faster, increase the drag force on fine particles attached to the soil skeleton, increase the rate of liquefied particle generation, and accelerate soil erosion; at the same time, the porosity of the soil increases under the action of erosion, the saturated permeability coefficient further increases, and the infiltration rate of rainwater is also accelerated. These two are in a positive correlation relationship. Therefore, the rate of liquefied particle generation is sensitive to changes in the saturated permeability coefficient.
Figure 21 shows the distribution of liquefied particle density along the height on the right side of the slope under different saturation coefficients. As the initial saturated permeability increases, the density of liquefied particles in the entire reference line increases, but the growth rate gradually slows down. When the initial saturated permeability coefficient is 5 × 10 5 m/s, the degree and velocity of erosion on the slope surface are significantly greater than those inside the slope. However, as the initial saturated permeability coefficient increases, the erosion velocities on the slope surface and inside the slope gradually become consistent.
As shown in Figure 22, there is no obvious difference in the final equivalent plastic strain distribution under different permeability coefficients. Due to the same rainfall intensity and large saturated permeability coefficients under the three conditions, the pressure head distribution in the slope is almost the same after 24 h. In the process of slope instability, the weight of the soil and the change in pore water pressure play a key role, so the final distribution of equivalent plastic strain is less different. From the above analysis, it can be seen that in the coupling model established in this paper, erosion only promotes seepage and has little effect on soil deformation.

5. Conclusions

In this paper, the NOSBPD-FEM hybrid analysis method for the seepage–stress coupling problems is established and applied to the slope stability analysis under rainfall conditions. By analyzing the evolution of the pore water pressure and equivalent plastic strain of the slope, the influence of different rainfall durations and rainfall intensities on the slope stability is systematically studied. Based on an improved erosion constitutive mode, the harm of rainfall erosion to slope safety under the effect of fluid–solid coupling is considered, and the influence of erosion and initial saturated permeability on slope stability is explored. This method provides a novel idea for slope stability analysis under rainfall conditions and is of great significance for further research on the process and mechanism of landslides. The main conclusions of this paper are as follows:
(1)
Under heavy rainfall conditions, transient saturated zones are easily formed on the slope surface. The resulting pore water pressure and the increased soil weight due to rainwater accumulation will harm slope stability. As the rainfall continues, the moist front moves downward, the transient saturated zone expands, and the slope stability will further decline. The foot of the slope is close to the rainfall infiltration boundary, and the pressure head rises faster than inside the slope, so damage will occur faster.
(2)
Rainfall intensity is also an important factor affecting slope stability. The higher the rainfall intensity, the faster the moist front moves downward. The seepage rainwater will be transferred to the deep part of the slope faster. The infiltration volume force, pore water pressure, and soil weight increase rate of the soil particles will increase, causing the slope to become unstable faster. Although the change in saturated permeability coefficient caused by soil deformation has an inhibitory effect on seepage, it has little effect on slope stability.
(3)
During rainfall, erosion gradually develops inward along the top and surface of the slope, and the erosion develops most rapidly at the slope toe. The equivalent plastic strain zone gradually expands. Due to the influence of erosion, the plastic strain at the slope toe tends to develop downward compared with the case where erosion is not considered. The initial saturated permeability coefficient has a great influence on the liquefied particle generation rate, and the two are in a positive correlation relationship.
It is worth noting that although the Gaussian probability model is used to modify the erosion model, it is still difficult to reflect the real random erosion process. In subsequent studies, a more accurate probability model can be considered. At the same time, the idealized slope drainage treatment method will also cause the simulation results to deviate from the actual situation. Subsequent studies will further consider the influence of slope waterlogging and slope surface runoff during rainfall.

Author Contributions

Conceptualization, methodology, and formal analysis, writing review and editing, X.G.; program implementation, simulation, and writing original draft, L.S.; conceptualization, writing review, and editing and data curation, X.X.; funding acquisition, C.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Water Resources Science and Technology Program of Hunan Province (No. XSKJ2023059-16) and the National Natural Science Foundation of China (No. 12002118, 11932006, 12172121).

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interests.

References

  1. Zong, J.; Zhang, C.; Liu, L.; Liu, L. Modeling rainfall impact on slope stability: Computational insights into displacement and stress dynamics. Water 2024, 16, 554. [Google Scholar] [CrossRef]
  2. Song, K.; Han, L.; Ruan, D.; Li, H.; Ma, B. Stability prediction of rainfall-induced shallow landslides: A case study of mountainous area in China. Water 2023, 15, 2938. [Google Scholar] [CrossRef]
  3. Li, D.; Zhou, N.; Wu, X.; Yin, J. Seepage–stress coupling response of cofferdam under storm surge attack in Yangtze estuary. Mar. Georesour Geotechnol. 2021, 39, 515–526. [Google Scholar] [CrossRef]
  4. Xu, X.; Xu, W.; Xie, C.; Ali Khan, M.Y. Prediction of the long-term performance based on the seepage-stress-damage coupling theory: A case in south-to-north water diversion project in China. Appl. Sci. 2021, 11, 11413. [Google Scholar] [CrossRef]
  5. Zhou, D.; Zhang, Z.; Li, J.; Wang, X. Seepage-stress coupled modeling for rainfall induced loess landslide. Theor. Appl. Mech. Lett. 2019, 9, 7–13. [Google Scholar] [CrossRef]
  6. Liu, C.; Yan, Y.; Yang, H.-Q. Numerical modeling of small-scale unsaturated soil slope subjected to transient rainfall. Geosyst. Geoenviron. 2023, 2, 100193. [Google Scholar] [CrossRef]
  7. Hua, C.; Yao, L.; Song, C.; Ni, Q.; Chen, D. A new criterion for defining inhomogeneous slope failure using the strength reduction method. Comput. Model. Eng. Sci. 2022, 132, 413–434. [Google Scholar] [CrossRef]
  8. Wang, T.; Hu, S.; Yang, M.; Meng, D. Numerical method for predicting the nucleation and propagation of multiple cracks based on the modified extended finite element method. Theor. Appl. Fract. Mech. 2024, 133, 104568. [Google Scholar] [CrossRef]
  9. Deng, Y.; Xia, Y.; Wang, D.; Jin, Y. A study of Hydraulic fracture propagation in laminated shale using extended finite element. Comput. Geotech. 2024, 166, 105961. [Google Scholar] [CrossRef]
  10. Hai, T.; Jamila, R.; Hakim, N. Incremental alternating algorithm for damage and fracture modeling using phase-field method. J. Mech. Sci. Technol. 2024, 38, 1385–1392. [Google Scholar]
  11. Stefan, L.; Verena, K.; Lukas, M.; Munk, L. An enriched phase-field method for the efficient simulation of fracture processes. Comput. Mech. 2023, 71, 1015–1039. [Google Scholar]
  12. Silling, S.A. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids 2000, 48, 175–209. [Google Scholar] [CrossRef]
  13. Jabakhanji, R. Peridynamic Modeling of Coupled Mechanical Deformations and Transient Flow in Unsaturated Soils; Purdue University: West Lafayette, IN, USA, 2013. [Google Scholar]
  14. Gu, Q.; Lin, Z.; Wang, L.; Qiu, Z.; Huang, S.; Li, S. A novel peridynamic solution for modelling saturated soil-pore fluid interaction in liquefaction analysis. Comput. Geotech. 2023, 162, 105686. [Google Scholar] [CrossRef]
  15. Zhou, Y.; Zhang, M.; Pei, W.; Wang, Y. A non-local frost heave model based on peridynamics theory. Comput. Geotech. 2022, 145, 104675. [Google Scholar] [CrossRef]
  16. Menon, S.; Song, X. Updated lagrangian unsaturated periporomechanics for extreme large deformation in unsaturated porous media. Comput. Methods Appl. Mech. Eng. 2022, 400, 115511. [Google Scholar] [CrossRef]
  17. Song, X.; Silling, S.A. On the peridynamic effective force state and multiphase constitutive correspondence principle. J. Mech. Phys. Solids 2020, 145, 104161. [Google Scholar] [CrossRef]
  18. Song, X.; Hossein, P. Computational cosserat periporomechanics for strain localization and cracking in deformable porous media. Int. J. Solids Struct. 2024, 288, 112593. [Google Scholar] [CrossRef]
  19. Menon, S.; Song, X. Shear banding in unsaturated geomaterials through a strong nonlocal hydromechanical model. Eur. J. Environ. Civ. Eng. 2022, 26, 3357–3371. [Google Scholar] [CrossRef]
  20. Liu, P.; Gu, X.; Lu, Y.; Xia, X.; Madenci, E.; Zhang, Q. Peridynamics for mechanism analysis of soil desiccation cracking: Coupled hygro-mechanical model, staggered and monolithic solution. Comput. Methods Appl. Mech. Eng. 2023, 406, 115896. [Google Scholar] [CrossRef]
  21. Ren, B.; Fan, H.; Bergel, G.; Regueiro, R.A.; Lai, X.; Li, S. A peridynamics-SPH coupling approach to simulate soil fragmentation induced by shock waves. Comput. Mech. 2014, 55, 287–302. [Google Scholar] [CrossRef]
  22. Liu, P.; Gu, X.; Lu, Y.; Xia, X.; Madenci, E.; Zhang, Q. A coupled hygro-mechanical model for moisture diffusion and curling mechanism in saturated and unsaturated soil using ordinary state-based peridynamics. Comput. Geotech. 2024, 172, 106473. [Google Scholar] [CrossRef]
  23. Majid, S.; Yan, H.; Jivkov Andrey, P. Peridynamics modelling of clay erosion. Geotechnique 2022, 72, 510–521. [Google Scholar]
  24. Lai, X.; Liu, L.; Liu, Q.; Cao, D.; Wang, Z.; Zhai, P. Slope stability analysis by peridynamic theory. Appl. Mech. Mater. 2015, 744–746, 584–588. [Google Scholar] [CrossRef]
  25. Zhang, T.; Zhang, J.-Z. Numerical estimate of critical failure surface of slope by ordinary state-based peridynamic plastic model. Eng. Fail. Anal. 2022, 140, 106556. [Google Scholar] [CrossRef]
  26. Wang, R.; Li, S.; Liu, Y.; Hu, X.; Lai, X.; Beer, M. Peridynamics-based large-deformation simulations for near-fault landslides considering soil uncertainty. Comput. Geotech. 2024, 168, 106128. [Google Scholar] [CrossRef]
  27. Zhou, H.; Shen, F.; Gu, X.; Li, B. Regularization effect of peridynamic horizon on strain localization and soil slope instability analysis. Eng. Anal. Bound. Elem. 2024, 165, 105774. [Google Scholar] [CrossRef]
  28. Yang, Y.; Liu, Y. Modeling of cracks in two-dimensional elastic bodies by coupling the boundary element method with peridynamics. Int. J. Solids Struct. 2021, 217–218, 74–89. [Google Scholar] [CrossRef]
  29. Ni, T.; Pesavento, F.; Zaccariotto, M.; Galvanetto, U.; Zhu, Q.Z.; Schrefler, B.A. Hybrid FEM and peridynamic simulation of hydraulic fracture propagation in saturated porous media. Comput. Methods Appl. Mech. Eng. 2020, 366, 113101. [Google Scholar] [CrossRef]
  30. Liu, S.; Fang, G.; Fu, M.; Yan, X.; Meng, S.; Liang, J. A coupling model of element-based peridynamics and finite element method for elastic-plastic deformation and fracture analysis. Int. J. Mech. Sci. 2022, 220, 107170. [Google Scholar] [CrossRef]
  31. Sun, W.; Fish, J. Coupling of non-ordinary state-based peridynamics and finite element method for fracture propagation in saturated porous media. Int. J. Numer. Anal. Methods Geomech. 2021, 45, 1260–1281. [Google Scholar] [CrossRef]
  32. Jin, S.; Hwang, Y.K.; Hong, J.-W. Stabilized non-ordinary state-based peridynamics with irregular nodal distribution. Mech. Res. Commun. 2023, 130, 104130. [Google Scholar] [CrossRef]
  33. Dahaghi, A.K.; Gholami, V.; Moghadasi, J. A novel workflow to model permeability impairment through particle movement and deposition in porous media. Transp. Porous Media 2011, 86, 867–879. [Google Scholar] [CrossRef]
  34. Yang, J.; Yin, Z.; Laouafa, F.; Hicher, P.Y. Hydromechanical modeling of granular soils considering internal erosion. Can. Geotech. J. 2020, 57, 157–172. [Google Scholar] [CrossRef]
  35. Golay, F.; Bonelli, S. Numerical modeling of suffusion as an interfacial erosion process. Eur. J. Environ. Civ. Eng. 2011, 15, 1225–1241. [Google Scholar] [CrossRef]
  36. Chang, D. Internal Erosion and Overtopping Erosion of Earth Dams and Landslide Dams; The Hong Kong University of Science and Technology: Hong Kong, China, 2012. [Google Scholar]
  37. Donatella, S. Effects of the erosion and transport of fine particles due to seepage flow. Int. J. Geomecanics 2003, 3, 111–122. [Google Scholar]
  38. Zhang, L.; Wu, F.; Zhang, H.; Zhang, L.; Zhang, J. Influences of internal erosion on infiltration and slope stability. Bull. Eng. Geol. Environ. 2019, 78, 1818–1827. [Google Scholar] [CrossRef]
  39. Gu, X.; Zhang, Q.; Yu, Y. An effective way to control numerical instability of a nonordinary state-based peridynamic elastic model. Math. Probl. Eng. 2017, 1–7. [Google Scholar] [CrossRef]
  40. Fredlund, D.G.; Rahardjo, H.; Fredlund, M.D. Unsaturated Soil Mechanics in Engineering Practice; John Wiley& Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  41. Zhao, C.G.; Liu, Y.; Gao, F.P. Work and energy equations and the principle of generalized effective stress for unsaturated soils. Int. J. Numer. Anal. Methods Geomech. 2010, 34, 920–936. [Google Scholar] [CrossRef]
  42. Dean, E.T.R. Discussion: Editorial. Géotechnique 2005, 55, 415–417. [Google Scholar] [CrossRef]
  43. Lewis, R.W.; Schrefler, B.A. The Finite Element Method in the Static and Dynamic Deformation and Consoildation of Porous Media; John Wiley& Sons: Hoboken, NJ, USA, 1999. [Google Scholar]
  44. Van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef]
  45. Nikolic, M. Discrete element model for the failure analysis of partially saturated porous media with propagating cracks represented with embedded strong discontinuities. Comput. Methods Appl. Mech. Eng. 2022, 390, 114482. [Google Scholar] [CrossRef]
  46. Callari, C.; Armero, F.; Abati, A. Strong discontinuities in partially saturated poroplastic solids. Comput. Methods Appl. Mech. Eng. 2010, 199, 1513–1535. [Google Scholar] [CrossRef]
  47. Nikolic, M.; Ibrahimbegovic, A.; Miscevic, P. Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities. Comput. Methods Appl. Mech. Eng. 2016, 298, 407–427. [Google Scholar] [CrossRef]
  48. Zhang, Y.; Madenci, E.; Gu, X.; Zhang, Q. A coupled hydro-mechanical peridynamic model for unsaturated seepage and crack propagation in unsaturated expansive soils due to moisture change. Acta Geotech. 2023, 18, 6297–6313. [Google Scholar] [CrossRef]
  49. Warrick, A.W.; Lomen, D.O.; Yates, S.R. A generalized solution to infiltration. Soil Sci. Soc. Am. J. 1985, 49, 34–38. [Google Scholar] [CrossRef]
  50. Zhao, M. Soil Mechanics and Foundation Engineering; Wuhan University of Technology Press: Wuhan, China, 2014. (In Chinese) [Google Scholar]
  51. Cividini, A.; Gioda, G. Finite-element approach to the erosion and transport of fine particles in granular soils. Int. J. Geomech. 2004, 4, 191–198. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of interaction between material points before and after deformation.
Figure 1. Schematic diagram of interaction between material points before and after deformation.
Water 16 02210 g001
Figure 2. Schematic diagram of slope dimensions.
Figure 2. Schematic diagram of slope dimensions.
Water 16 02210 g002
Figure 3. Equivalent plastic strain distribution under stable slope condition.
Figure 3. Equivalent plastic strain distribution under stable slope condition.
Water 16 02210 g003
Figure 4. Top view of landslide damage.
Figure 4. Top view of landslide damage.
Water 16 02210 g004
Figure 5. Fluid–structure coupling scheme flow chart.
Figure 5. Fluid–structure coupling scheme flow chart.
Water 16 02210 g005
Figure 6. Geometric model of soil column.
Figure 6. Geometric model of soil column.
Water 16 02210 g006
Figure 7. Comparison diagram between finite element method solution and analytical solution of pressure head variation with soil column depth.
Figure 7. Comparison diagram between finite element method solution and analytical solution of pressure head variation with soil column depth.
Water 16 02210 g007
Figure 8. Comparison diagram between analytical and numerical solutions of pore water pressure distribution along the depth of soil column.
Figure 8. Comparison diagram between analytical and numerical solutions of pore water pressure distribution along the depth of soil column.
Water 16 02210 g008
Figure 9. Comparison between Abaqus solution and self-programmed solution for vertical displacement distribution along soil column height.
Figure 9. Comparison between Abaqus solution and self-programmed solution for vertical displacement distribution along soil column height.
Water 16 02210 g009
Figure 10. Initial saturation distribution of slope.
Figure 10. Initial saturation distribution of slope.
Water 16 02210 g010
Figure 11. The pressure head distribution of the slope changes with time when the rainfall intensity is 5 mm/h.
Figure 11. The pressure head distribution of the slope changes with time when the rainfall intensity is 5 mm/h.
Water 16 02210 g011
Figure 12. The variation of slope equivalent plastic strain distribution with time when the rainfall intensity is 5 mm/h.
Figure 12. The variation of slope equivalent plastic strain distribution with time when the rainfall intensity is 5 mm/h.
Water 16 02210 g012
Figure 13. Slope pressure head distribution under different rainfall intensities (m).
Figure 13. Slope pressure head distribution under different rainfall intensities (m).
Water 16 02210 g013
Figure 14. Distribution diagram of equivalent plastic strain of slope under different rainfall intensities.
Figure 14. Distribution diagram of equivalent plastic strain of slope under different rainfall intensities.
Water 16 02210 g014aWater 16 02210 g014b
Figure 15. Equivalent plastic strain at the slope foot over time under different rainfall intensities.
Figure 15. Equivalent plastic strain at the slope foot over time under different rainfall intensities.
Water 16 02210 g015
Figure 16. Schematic diagram of fine particle erosion from initiation to development.
Figure 16. Schematic diagram of fine particle erosion from initiation to development.
Water 16 02210 g016
Figure 17. Density distribution of liquefied particles in the slope under 10 mm/h rainfall (g/m3).
Figure 17. Density distribution of liquefied particles in the slope under 10 mm/h rainfall (g/m3).
Water 16 02210 g017
Figure 18. Distribution of liquefied particle density along the height on the left side of the slope under 10 mm/h rainfall.
Figure 18. Distribution of liquefied particle density along the height on the left side of the slope under 10 mm/h rainfall.
Water 16 02210 g018
Figure 19. Equivalent plastic strain distribution diagram in the slope under 10 mm/h rainfall.
Figure 19. Equivalent plastic strain distribution diagram in the slope under 10 mm/h rainfall.
Water 16 02210 g019
Figure 20. Density distribution diagram of liquefied particles in the slope under different initial saturated permeability coefficients (g/m3). (a) Initial saturated permeability 5 × 10 5 m/s; (b) initial saturated permeability 1 × 10 4 m/s; (c) initial saturated permeability 1.5 × 10 4 m/s.
Figure 20. Density distribution diagram of liquefied particles in the slope under different initial saturated permeability coefficients (g/m3). (a) Initial saturated permeability 5 × 10 5 m/s; (b) initial saturated permeability 1 × 10 4 m/s; (c) initial saturated permeability 1.5 × 10 4 m/s.
Water 16 02210 g020
Figure 21. Distribution of liquefied particle density along the height on the right side of the slope under different initial saturated permeability coefficients.
Figure 21. Distribution of liquefied particle density along the height on the right side of the slope under different initial saturated permeability coefficients.
Water 16 02210 g021
Figure 22. Equivalent plastic strain distribution diagram in slope under different initial saturated permeability coefficients. (a) Initial saturated permeability 5 × 10 5 m/s; (b) initial saturated permeability 1 × 10 4 m/s; (c) initial saturated permeability 1.5 × 10 4 m/s.
Figure 22. Equivalent plastic strain distribution diagram in slope under different initial saturated permeability coefficients. (a) Initial saturated permeability 5 × 10 5 m/s; (b) initial saturated permeability 1 × 10 4 m/s; (c) initial saturated permeability 1.5 × 10 4 m/s.
Water 16 02210 g022
Table 1. Calculation parameters of seepage–corrosion coupling model.
Table 1. Calculation parameters of seepage–corrosion coupling model.
ParameterMagnitudeParameterMagnitude
ρ f 0 0.4 g·cm−3 θ r 0
ρ f 0.386 g·cm−3 n 0 0.25
ρ s 2.368 g·cm−3 α w 0.08 kPa−1
v * 3.6 × 10 6 m·s−1 n w 2
α er 4.76 c 13.03 kPa
β er 6.95 × 10 3 φ 21
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Gu, X.; Song, L.; Xia, X.; Yu, C. Finite Element Method-Peridynamics Coupled Analysis of Slope Stability Affected by Rainfall Erosion. Water 2024, 16, 2210. https://doi.org/10.3390/w16152210

AMA Style

Gu X, Song L, Xia X, Yu C. Finite Element Method-Peridynamics Coupled Analysis of Slope Stability Affected by Rainfall Erosion. Water. 2024; 16(15):2210. https://doi.org/10.3390/w16152210

Chicago/Turabian Style

Gu, Xin, Laike Song, Xiaozhou Xia, and Cheng Yu. 2024. "Finite Element Method-Peridynamics Coupled Analysis of Slope Stability Affected by Rainfall Erosion" Water 16, no. 15: 2210. https://doi.org/10.3390/w16152210

APA Style

Gu, X., Song, L., Xia, X., & Yu, C. (2024). Finite Element Method-Peridynamics Coupled Analysis of Slope Stability Affected by Rainfall Erosion. Water, 16(15), 2210. https://doi.org/10.3390/w16152210

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop