Smoothed Particle Hydrodynamics (SPH) Analysis of Slope Soil–Retaining Wall Interaction and Retaining Wall Motion Response

: The occurrence of slope instability disasters seriously endangers the safety of people’s lives and property in China. Therefore, it is essential to study the slope instability process and the interaction between soil and retaining walls. In this paper, the smoothed-particle hydrodynamics (SPH) method, based on the elastoplastic constitutive model of rock and soil, was used to simulate the entire process of slope instability and the interaction between soil and retaining walls. The model, based on the classical elastic–plastic theory, includes linear elastic deformation and plastic deformation following the non-associated ﬂ ow rule under the Drucker–Prager (DP) yield criterion. By considering the plastic characteristics of geotechnical materials, this method can accurately simulate the slope movement process. The model was established, calculated, and compared with a slope example, thus verifying its feasibility. Furthermore, the motion response of the retaining wall under di ﬀ erent conditions was studied, which provides a new numerical simulation platform for the stability checking of the retaining wall and motion analysis after the interaction between the retaining wall and slope soil.


Introduction
Natural slopes are very common in complicated environments; however, many artificial slopes are produced by large-scale engineering construction and excavation [1][2][3].Landslides are a widespread geological disaster globally, causing thousands of deaths and billions in property loss [4].The study of slope instability and the interaction between slope soil and retaining structure plays an important role in analyzing the whole instability process.The problem of slope instability and its failure involves the process of large soil deformation and soil-structure interaction, which has been the research object of many scholars in recent years.The study of slope stability analysis, motion mechanism, and motion process is mainly conducted in two ways: first, numerical methods based on grids, such as the finite element method (FEM) and the finite difference method (FDM), are adopted to mainly solve the problems identified via slope stability analysis and deformation analysis; the other involves meshless methods, such as the discrete element method (DEM) and the material point method (MPM) [5].Geotechnical engineers usually use analytical and empirical methods to develop an appropriate model to effectively simulate the engineering site according to the design parameters and engineering characteristics of soil or rock materials [6].Kumar, V et al. [7] used the finite element method to analyze stability and added the shear strength reduction method to the finite element method to quantify the existing slope strength and infer the strain and displacement patterns.By analyzing the landslide activity mode, the shear strain and displacement of the slope were determined, and the harmfulness of the landslide zone was evaluated.N. Ko- take et al. [8] used the finite element method (FEM) to simulate a series of laboratory model test results of unreinforced sand slope and reinforced sand slope.The results of the FEM simulation of unreinforced sand slope and reinforced sand slope are consistent with the physical experiment results in terms of the load-settlement relationship, strain field showing failure mode, and tension distribution of steel bars in different reinforcement arrangements.Singh, SK et al. [9] used the finite difference method (FDM) and the limit equilibrium method (LEM) to conduct numerical analysis on typical slopes common in Indian railways to determine the stability of safety factors under static and seismic conditions.Maugeri, M et al. [10] put forward a mathematical model to predict the viscous deformation of landslides based on the FDM method, which has been verified in a landslide case in Gagliano Castelferrato, Italy.However, the numerical method based on the grid will appear as grid distortion, which will lead to large calculation errors or even nonconvergence [11].It is necessary to use a meshless method to simulate the large deformation of a slope in specific scenes.Li, SH et al. [12] studied the critical excavation depth of two ideal jointed rock slopes by using the DEM method.In addition, the influence of joints on failure modes in the DEM simulation and experimental observation was compared.It was found that the critical excavation depth predicted by DEM is lower than that predicted by LEM if the joint structure in the rock mass is not considered.Wang, HB et al. [13] put forward a discrete element modeling method (DSDM) based on displacement statistics.The stability of jointed rock slopes was analyzed by the coupled shear strength reduction method (SRM).Its innovation lies in using the global displacement statistics around the potential slip surface of jointed slopes instead of the local displacements of several characteristic points to define the criteria for terminating the iterative calculation process of DEM and determining FOS.Therefore, the DSDM method can effectively analyze the stability of jointed slopes with potential slip surfaces.Feng, ZK et al. [14] developed a high-performance MPM based on GPU and introduced the strength reduction method (SRM) for a quantitative analysis of slope failure risk assessment.By comparing the results of MPM with those of the limit equilibrium method (LEM), the results showed that the accuracy of MPM-SRM for landslide failure mechanics and deformation characteristics exceeds LEM, and based on this method, the stability and failure mechanism of east temple slope downstream of Xiaolangdi Dam in China were analyzed.Wu, FY et al. [15] established the foundation friction algorithm based on the MPM method and calculated the foundation friction considering pore water pressure.The effectiveness of MPM in simulating the dynamic process of soil landslide was verified by a large-scale debris flow test, and the depth of debris flow calculated by MPM is consistent with the experimental results.
SPH was originally proposed to solve astrophysical problems in three-dimensional open space and has been widely studied and expanded [16][17][18]; it has been applied to dynamic response problems with material strength and hydrodynamic problems with large deformation.SPH is mesh-free, immune to the mesh distortion problem in FEM/FDM, and able to directly simulate the large deformation of soils that occurs during landslides.A series of slope stability analyses were performed using an in-house SPH program [19].Yu, M et al. [20] used the SPH model to predict the mobility of unstable slopes after failure and extended the application of the SPH model from disaster simulation to the prediction and analysis of unstable landslides.The stability of the SPH method was proved by comparing the two groups.Mori, H et al. [21] modeled the shear strength parameters of soil as random variables and simulated the random field to explore the influence of soil variability on the runout distance.In addition, the uncertainty of rainfall characteristics was expressed by the Gumbel distribution, and the subsequent rainfall infiltration was simulated in a multiple seepage analysis to obtain the pore pressure profile of the slope, which was used as the initial condition of the SPH method.Combined with these different sources of uncertainty, according to the response uncertainty of landslide runout distance, the risk factors indicating the risk of nearby structures were quantified.
Lian, YJ et al. [22] put forward the first and fully validated smoothed-particle hydrodynamics (SPH) model to solve the flow-deformation coupling problem in unsaturated porous media that undergoes large deformation and post-failure behavior.The ability of the proposed SPH model was verified by a basic consolidation test and large-scale rainfallinduced slope failure test.The theoretical solution was in good agreement with the experimental results, which shows that the proposed SPH model can be easily extended to solve large-scale geotechnical engineering applications involving unsaturated seepage deformation coupling problems.Hu, YX et al. [23] established a coupled discrete element method (DEM) and smoothed-particle hydrodynamics (SPH) method to predict the whole process of landslide, wave generation, and wave propagation on the top slope of the Houziyan Reservoir in Sichuan.The DEM-SPH coupling simulation method proposed provides important enlightenment for reservoir landslide prevention and control.
The use of the SPH method has made good progress in the study of slope problems, but most of the current research focuses on the deformation of the soil after slope instability, and research on the impact of slope soil on the environment after instability is lacking.To further study the influence of slope instability, a soil-retaining wall interaction model was proposed by using the SPH method combined with the soil elastic-plastic constitutive model, and the model was verified.In addition, the motion law of retaining walls under different initial conditions was analyzed, which provides a new idea for the motion response analysis of retaining walls.

SPH Theory
The core idea of the SPH method is to describe a continuous fluid or solid by utilizing interacting particle groups, with each particle representing a specific region within the problem domain of SPH.Each particle carries various physical parameters including density, mass, velocity, and acceleration.By solving the control equation for the particle groups, the mechanical behavior of the entire system can be obtained.The SPH solution consists of two key steps: the first step involves the approximation of the kernel function, and the second step involves particle approximation.These steps are as follows: Kernel function approximation: The field functions (such as mass, velocity, density, etc.) are approximated by the following equation: (1) where f is a function of particle arbitrary space vector x; is Dirac function. ( can only approximate the space vector function at one point, so the SPH discrete model cannot be constructed.If the smoothing function w is used instead of , the formula above can be expressed as follows: ( where h is the smooth length, which indicates the influence range of the smooth function-that is, the support domain.The choice of smoothing function has a great influence on the accuracy, efficiency, and stability of the calculation results.The smoothing function in this paper adopts the cubic spline function [24]: where q is the relative distance between particles from to , ; is a fixed coefficient.In one-dimensional, two-dimensional, and three-dimensional problems, there are , , .The aim is to find the spatial derivative of Equation (3), and the derivative integral expression of the function is as follows: (5) Particle approximation: The partial differential equations to be solved in fluid force problems are transformed into a series of discrete ordinary differential equations by particle approximation.(6) where n is the number of particles in the support domain; is the particle volume, which is used to replace the infinite body element at the particle j; and represent the mass and density of j, respectively.
Similarly, the particle approximate form of the function space derivative is as follows: (7) where is the gradient of the nucleus, , is the distance between particles i and j.
According to Equations ( 6) and ( 7), the SPH form of the Navier-Stokes equation can be obtained.(8) where t is time, v is the velocity of particles, is the total stress tensor of the particle, x is the displacement of the particle, and f is the external force on the particle.is artificial viscosity, which is used to reduce the non-physical oscillation in the numerical results around the impact area [25].Equation ( 8) expresses the motion of SPH particles, following the laws of conservation of mass and momentum.(9) where , , , , c is the speed of sound, and α and β are constants and should be adjusted according to specific examples.When the SPH method simulates the mechanical behavior of solid materials, it will cause tension instability.In this paper, the artificial stress method is introduced, and its core idea is to introduce a kind of short-range repulsive force between a pair of adjacent particles to prevent tension when they come too close to each other [26]. ( where i and j represent a pair of adjacent particles; is the distance between two particles. For two-dimensional problems, the components of artificial stress tensor are as follows: (11) (12) (13) where is the angle. ( The diagonal component of the artificial stress tensor can be expressed as follows: (15) where is a constant coefficient between 0 and 1. is also calculated by a similar expression.The stress tensor of diagonal components, and , after rotating the coordinate system can be calculated by the following formula: (16)

Constitutive Model
In this paper, the constitutive model of soil adopts the elastic-plastic model, and its total stress tensor consists of elastic and plastic parts.(17) where is the total stress tensor, is the shear stress tensor, and p is the hydrostatic pressure.In this paper, according to its standard definition, hydrostatic pressure p is obtained from this constitutive equation: (18) where the compressive stress is negative.For ideal elastic-plastic materials, the strain-rate tensor is as follows: ( where consists of the elastic strain-rate tensor and fractional plastic strainrate tensor . ( According to Hooke's law, , where is deviatoric shear stress tensor, G is shear modulus ), and v is Poisson's ratio. ( The plastic strain-rate tensor can be calculated by the plastic flow law, as follows: (23) where is the plastic strain-rate operator, and g is the plastic buckling function.The total strain-rate tensor is expressed in the form of the stress-rate tensor, as follows: (24) In this paper, the Drucker-Prager yield criterion is adopted to determine the plastic flow state of the soil.The yield conditions of this criterion are as follows: (25) where and are the first and second invariants of the stress tensor, and , .
, and are Drucker-Prager constants, which can be obtained by mutual conversion with the parameters of the Mohr-Coulomb material, such as cohesion c and internal friction angle .In the plane-strain state, the Drucker-Prager constant is defined as , .
For the non-associated flow law, the plastic yield function is expressed as (26) where ψ is the dilatancy angle, where it is zero, indicating the plastic incompressibility of the material.The change rate of the plastic strain operator is obtained by solving (27) After various calculations, the stress-strain relationship of the Drucker-Prager ideal elastic-plastic geotechnical material model under the non-associated flow law can be expressed as follows: (28) Because this paper considers the problem of large deformation, it is necessary to calculate the differential strain rate of materials (the Jaumann rate): (29) where is the spin-rate tensor Finally, various calculated stress-strain relationships are as follows [27]: (31)

Boundary Conditions
As the difficulty of SPH, boundary conditions greatly affect the calculation accuracy..There are two kinds of boundary treatment methods: by distributing virtual particles on the boundary, the neighboring particles are repelled to prevent them from penetrating the boundary; virtual particles fill the boundary to reflect the symmetrical surface so that neighboring particles are supported and covered by particles.Morris et al. proposed the boundary condition of velocity without slip [28].When the distance between the virtual particle and the material particle is less than the threshold , the virtual particle velocity is taken as the inverse number of the material particle, and the pressure and stress are kept constant.
(32) where , , and represent soil velocity, pressure, and stress when the spacing is less than the threshold., , and represent the velocity, pressure, and stress of virtual particles.
(33) where is the velocity coefficient, and is the maximum velocity coefficient, which is used to prevent the material from being too fast when it is close to the boundary.

Morris's research shows that when
, the simulation result is better.According to the research of Bui et al., is between 1.5 and 2.0 [27].In this paper, when is 1.6, the effect is better.

Contact Model
In this paper, the soil-retaining wall interaction model for slopes is proposed, and the soil and retaining wall are represented by uniformly distributed particles.In the process of contact between the soil and retaining wall, the retaining wall is always regarded as a rigid body and does not deform.The external force types of soil particles are gravity and the contact force of the soil and the retaining wall.
As shown in Figure 1, a series of particles are generated on the surface of the retaining wall (k is adjacent to k + 1).When dp is less than the threshold d0, it is determined that contact occurs, and the vertical vector and horizontal vector of point P can be calculated by the following formula: where χ represents the penetration degree of particles (χ is 0, which means that penetration is not allowed), µ is the friction coefficient, and ▽t is the time step.If , then [29].Contact and stress analysis of soil-retaining wall is shown in For rigid body motion, its speed can be divided into translation speed and rotation speed, and the conservation equation (i.e., control equation) of rigid body motion can be expressed by the following formula: (36) where M0, I, v0, Ω0, and X0 are the mass, moment of inertia, translation speed, rotation speed, and centroid position of the rigid body, and Nc is the number of particles interacting with the rigid body.
According to the basic theoretical analysis of SPH and elastic-plastic constitutive model, Fortran language is used to simulate the whole process of slope and the interaction between slope soil and retaining wall.The initial conditions need to be set when the program runs-particle spacing, smooth length, density, mass, and coordinate information.The program runs based on the Visual Studio environment.The implementation process of the program is shown in Figure 2.

Experimental Control
The elastic-plastic constitutive model has been verified in previous research [26].In this paper, the soil-retaining wall interaction model of the slope was verified by the experiments of Taiping Mou et al. [30].The test model is shown in Figure 3.The test model consists of two parts: retaining wall and soil slope.The soil slope was made of silty clay, and the vertical retaining wall was made of cement mortar.In this experiment, geotechnical centrifuge was used to continuously agitate the model.Using centrifuge can increase the gravity acceleration of the model so as to simulate some special situations, such as earthquakes.By observing the test results, it could be found that the retaining wall only moved horizontally during the whole test process but did not rotate.Through the theoretical method described in the start of the article, the SPH particle model was constructed.Soil parameters, such as dilatancy angle, shear modulus, bulk modulus, etc., have a certain universality according to references in the literature [31].The main parameters of clay are as follows: internal friction angle = 35°, cohesive strength c = 23.93KN/m 2 , density = 1330 g/cm 3 , and dilation angle ψ = 13.95°.The density of the retaining wall refers to the density of concrete, which was 2400 Kg/m 3 here.When the soil type is silty soil, the friction coefficient µ between the soil and the bottom of the retaining wall is usually 0.25-0.35(The friction coefficients of different materials are shown in Table 1.), and in this model, it was 0.3.
The SPH model is established according to the above experimental standards, and the details of the SPH model are shown in Figure 4.
Table 1.Friction coefficient between rock and soil and the bottom of the retaining wall.

Types of Rock and Soil
Friction Factor Clay 0. The SPH simulation results for when the model was not centrifuged are shown in Figure 5.The arrow in the figure represents the velocity vector of particles.Figure 5 shows that the retaining wall with a normal design can bear the soil pressure effectively without translation.
In the SPH model, we used increasing gravity acceleration to simulate the process of the centrifuge used.By gradually increasing the acceleration of gravity, the retaining wall can be simulated to bear greater load and stress, which usually occurs in specific scenes, such as earthquakes and large-scale mechanical operations on the slope.When the acceleration of gravity increases to 3 g, the retaining wall will eventually be unable to resist such a load and translate.The simulation process and result analysis are shown in Figures 6 and 7.
The translation speed of the retaining wall under the action of earth pressure is positively correlated with the gravity acceleration of the model.After 0.84 s, the speed of the retaining wall gradually decreased to zero, and the displacement at this time was 1.43 mm.The simulation results are consistent with the experimental results, which verifies the accuracy of the model.
In order to study the influence of different soil environments on the translation of retaining walls, we used the friction coefficients of different types of soils to simulate the movement behavior of retaining walls.In this paper, we chose the friction coefficient µ = 0.2 and 0.1 for simulation and the simulation process and result analysis are shown in Figures 8-11   When the friction coefficient µ is 0.2, the results show that the retaining wall will be slightly displaced by soil force, and the final displacement is 0.21 mm.
When the friction coefficient µ is 0.1, the friction force on the retaining wall is not enough to support it to keep balance.In this case, the retaining wall will translate after being subjected to soil force, and the maximum translation distance is about 2.59 mm at 6.1 s.It is worth noting that the retaining wall shows a slight overturning trend at this moment.This result reminds us that, in practical engineering, if the friction force is not enough to keep the balance of the retaining wall, we need to use other methods to ensure the stability of the retaining wall.For example, in the event of an earthquake, it is difficult for the friction coefficient between the soil and the retaining wall to keep the balance of the retaining wall; therefore, the buried depth of the retaining wall should be increased in earthquake-prone areas.
The failure of the retaining wall has two forms: translation and overturning.Therefore, when designing retaining walls, it is necessary to carry out anti-translation checking and anti-overturning checking at the same time.Through a comparison between the model and the experiment above, the applicability of this SPH numerical model to translational failure is verified.This paper will further study the failure of overturning.The anti-overturning stability coefficient of the retaining wall needs to be calculated in the antioverturning checking calculation.The horizontal projection width b of the retaining wall bottom is the key parameter.This paper will analyze the overturning state of the retaining wall under different b values, and then verify the applicability of this model to the overturning failure of the retaining wall.The formula for checking the overturning resistance of retaining wall is as follows: (37) where is the overturning stability coefficient, G is the weight of the retaining wall, and is the horizontal distance from the center of the retaining wall to the toe of the wall. is the resultant force of dynamic earth pressure.
is the angle between the back of the wall and the horizontal projection of the wall bottom.
is the friction angle between the back of the wall and rock and soil.b is the horizontal projection length of the retaining wall bottom surface.z is the vertical distance from the point where the geotechnical pressure acts to the heel of the wall.
is the inclination angle of the bottom surface of the wall.In this paper, the bottom of the wall is placed horizontally; therefore is 0. Through the model above, we know that when µ is greater than 3, the retaining wall will not translate.The friction coefficient µ is controlled above 3, and other parameters remain unchanged.After changing b, the overturning state of the retaining wall is analyzed.When considering the original b value as 1, b is reduced to 0.9-0.7 times; the simulation process is shown in Figure 12.To verify whether the overturning results of the retaining wall in our model meet the design specifications, we further analyzed the included angle between the bottom edge of the retaining wall and the ground under different b values.For a more visual representation of the results, we imported the simulation results into CAD software 2019 for visualization.The result analysis is shown in Figure 13.
The analysis revealed that when the b value is 0.9 times the original value, the angle between the bottom edge of the retaining wall and the ground is approximately 2.70 degrees.When the b value is reduced to 0.8 times the original value, the angle increases to approximately 3.59 degrees.Further decreasing the b value to 0.7 times the original value results in a further increase in the angle to approximately 4.61 degrees.
The simulation results show that, with the decrease in the b value, the overturning degree of the retaining wall is also great, which proves that the model conforms to the anti-overturning checking code of the retaining wall.
This result serves as a reminder that, in order to prevent the overturning of the retaining wall, it is necessary to increase the length of the bottom edge in the design.If the conditions do not allow for an increase in the length of the bottom edge, other design aspects should be reinforced to ensure the stability of the retaining wall against overturning.One approach could be to increase the self-weight of the retaining wall.

Conclusions
In this paper, SPH modeling was carried out based on the elastic-plastic constitutive model, and an interaction model of slope soil and retaining wall was put forward.The simulation results were compared with the experiment, and the following conclusions were obtained: 1.By comparing the simulation results with the test results, the applicability of this model to the translation failure of retaining walls was verified.Displacement under different base friction coefficients was analyzed, which further verified that the model conformed to the anti-translation checking conditions in the code.2. SPH analysis was carried out on the motion response of retaining walls with different base lengths.The overturning situation and overturning angle of different retaining walls were compared.It was verified that the model conformed to the anti-overturning checking conditions in the code.3. Overall, the model accurately describes the motion response of the retaining wall under different conditions after soil impact.The model proved to be suitable for practical engineering cases and can provide a brand new numerical simulation platform for the design and checking calculation of retaining walls and the impact of soil on structures after landslides.

Figure 1 .
Figure 1.Soil-retaining wall contact and stress analysis diagram.In the figure, pink particles represent soil particles and black particles represent the retaining wall particles.The trapezoid in the figure represents the retaining wall.

Figure 2 .
Figure 2. Flow chart of SPH model algorithm.

Figure 3 .
Figure 3. Model diagram of the landslide experimental test of Mou et al [30].

Figure 4 .
Figure 4. Schematic diagram of the slope and retaining wall.

Figure 7 .
Figure 7. Time-displacement diagram of retaining wall with acceleration of 3 g.

Figure 13 .
Figure 13.Relationship diagram of the angle between the bottom edge of the retaining wall and ground.