Dynamic Response Analysis of Retaining Dam under the Impact of Solid-Liquid Two-Phase Debris Flow Based on the Coupled SPH- DEM-FEM Method

Based on the coupled SPH-DEM-FEM numerical method, this paper analyzes the dynamic interaction of solid debris flow particleliquid debris flow slurry-retaining dam in order to explore the dynamic response of retaining dam under the impact of the solidliquid two-phase debris flow and delves into the process of the debris flow impact on the dam, the impact force of debris flow, and the elastic-plastic time-history characteristics of the dam under different slopes of trapezoidal grooves. The calculation results show that the coupled SPH-DEM-FEMmethod can vividly simulate the impact behavior of the solid-liquid two-phase debris flow on the dam, reproduce the impact, climbing, and siltation in the process of the debris flow impact; the dynamic time-history curve of the retaining dam is consistent with the law of the literature, and the result of the debris flow impact force obtained is close to that of the empirical formula. Moreover, this paper studies the impact force distribution of the debris flow impact process. The results have a certain reference value for the study of the dynamic response of the retaining dam under the impact of the solid-liquid two-phase debris flow and the engineering design of the debris flow-retaining dam.


Introduction
As a very common natural disaster in mountainous areas, debris flow is a special flood with a large number of solid particles and liquid slurry produced by the provenance body under the action of rainstorm and earthquake [1]. Debris flow occurs in various mountain areas of the world, and the enormous impact force it produces can cause structural damage to all kinds of buildings and structures [2]. According to statistics, from 1950 to 2011, a total of 213 major debris flow disasters were recorded, resulting in a total of 77,779 deaths [3]. In particular, the August 7 Debris Flow Disaster occurred in Zhouqu of Gansu Province in China, which is one of the worst geological disasters in the world, left 1,557 people dead and 284 missing and caused a direct economic losses as high as 400 million. Data show that the fast-moving solid-liquid two-phase debris flow can unleash devastating power to destroy a large number of structures [4].
In order to mitigate the damage caused by debris flow, a blocking structure is usually set up in the flowing area of debris flow, which can effectively block the fluid, lower the flow velocity, and reduce the damage downstream [5][6][7]. The impact of the solid-liquid two-phase debris flow should be taken into consideration in the design of debris flowretaining dam; otherwise, damage can be caused to the dam, and the safety of life and property downstream can be seriously threatened. Therefore, the dynamic coupling of solid debris flow particles and liquid phase debris flow slurry and retaining dam should be considered in the design of debris flow-retaining dam.
According to the existing literatures, the motion analysis of debris flow and the dynamic response analysis of retaining dam are conducted mostly through field investigation, laboratory test, theoretical analysis, and numerical simulation. First-hand data obtained from field investigation is the basis for analysis. The advance of Beidou, remote sensing, UAV, and other technologies have greatly facilitated field investigation [8][9][10]. The laboratory test is easier to carry out than the field test, but it is difficult to achieve due to the similarity of the physical simulation test, and the coupling effect of solid debris flow particles is often overlooked [11,12]. In 1961, the Chinese Academy of Sciences carried out a large number of on-site debris flow observations in Jiangjia Ravine, Yunnan Province, and obtained a large amount of precious data [13]. Chen [14] et al. simplified the debris flow to a twodimensional fluid model and put forward a formula for calculating the impact force. The calculation of the impact force of the large rock mainly includes simplified calculation based on cantilever beam, calculation based on the energy method, and calculation based on Hertz elastic theory [15]. In recent years, numerical analysis methods have been widely used. In general, there are two kinds of numerical analysis methods: grid-based and particle-based, and there are two description methods: Euler description and Lagrangian description. The finite element method is a common continuous medium calculation method based on grid, but in the case of irregular bottom bed, complex mesh deformation is common in debris flow simulation, which is easy to cause mesh distortion, resulting in inaccurate calculation and even failure. Smooth particle hydrodynamics (SPH) is one of the earliest Lagrangian meshless continuum methods. It was initially used by Lucy and Gingold to solve celestial problems, but has been widely used in fluid calculation due to its good adaptability to solve large deformation problems. The discrete element method (DEM) was first proposed by Cundall for rock mechanics and later gradually applied to the field of soil mechanics. Wang [15] et al. used the SPH method to simulate the dynamic behavior of debris flow liquid, but did not carry out fluid-solid coupling analysis. Zhao [16] simulated the movement process of viscous debris flow based on the SPH method.
The hybrid numerical method is often used to analyze the coupling between debris flow and retaining dam. Liu et al. [17,18] used the SPH-FEM method and DEM-FEM method to analyze the dynamic behavior of debris flow and retaining dam and the coupling effect of debris flow and flexible retaining net. It reflects the applicability and accuracy of the SPH-DEM-FEM coupling numerical method to simulate solid-liquid two-phase debris flowretaining dam coupling analysis. The schematic diagram of SPH-DEM simulation of the solid-liquid two-phase debris flow is shown in Figure 1.
To sum up, the dynamic response analysis of the retaining dam under the impact of the solid-liquid two-phase debris flow is rarely studied. The existing literatures focus only on large rock impact dam, debris flow slurry impact dam, or debris flow slurry impact dam, while the coupling analysis of the solid-liquid two-phase flow debris flow-retaining dam is largely overlooked. In this paper, the numerical analysis method is used to study the dynamic response of the retaining dam under the impact of the solid-liquid two-phase debris flow, which provides theoretical and technical support for the design of the debris flow prevention and control system.

Numerical Computation Method
As computing technology advances in recent years, a variety of numerical analysis software has been used to solve many complex physical problems. Fluid-solid coupling analysis is a case in point. The three-dimensional finite element software ABAQUS can be used to analyze the coupling between debris flow and structure.

DEM Governing Equation.
The solid debris flow particles are modeled by the discrete element method (DEM). It is assumed that the solid debris flow particles are represented by elastic spheres of different sizes, and the irregularly shaped particles can be formed by several spheres. Such setting not only simplifies the complexity of debris flow particles but also simulates the interaction with debris flow slurry and retaining dam to the maximum extent [19].
The collision and energy dissipation of particles will occur in the process of movement. The force expression of particles considering damping is as shown in formula (1): In formula (1), k n , k t is the normal and tangential stiffness, respectively; λ, δ is the normal and tangential overlap, respectively, and nis the unit vector.
The resultant of all the external forces and the resultant moment of all the external forces on the particles are shown in formula (2): In formula (2), k is the total amount of particle j in contact with particle i, and f s is the friction force between particles.

SPH Governing Equation.
SPH is a meshless particle method based on Lagrangian frame, so it overcomes the mesh distortion of the traditional grid method and has good adaptability to fluid motion problems under complex boundary conditions. The approximate calculation of the fluid equation between particles through the kernel function, in essence, is the use of the kernel density estimation method to search and calculate the attributes of the neighboring particles in the particle i support domain by the summation method.
The kernel approximation of the SPH method primarily refers to the use of smoothing function and its derivative, and the kernel function is shown in formula (3): In formula (3), δðx − x ′ Þ is a Dirac function; f ðx ′ Þ is a function of position vector x.
The particle approximate form of the SPH method is shown in formula (4): In formula (4), l is the smooth length, ∇ w is the gradient of particles, and Nis particles in the support domain.

SPH-DEM-FEM Coupling Control Equation.
The DEM-FEM coupling control equation is shown in formula (5): In formula (5), f n,ij and f t,ij are normal contact force and tangential contact force;f a is the external load matrix of finite element; f b is the contact force matrix between finite element and discrete element.
The SPH-FEM coupling control equation is shown in formula (6): In formula (6), F is the external load matrix, and F c is the contact force matrix between SPH and FEM.
The SPH-DEM coupling control equation [20] is shown in formula (7): In formula (7), s and f represent DEM particles and SPH particles, respectively;

Optimization of Boundary Problems Based on the Virtual
Particle Method. In response to the SPH boundary penetration problem, this paper uses the virtual particle method proposed by Chen et al. [20] to optimize boundary problems. This method comprises the following steps: three or four layers of virtual particles are arranged along the boundary curve on the outside of the fluid according to the length of the smooth particles of the fluid. Virtual particles can be regarded as the expansion of fluid particles and conditionally involve in the calculation of continuity equations. When fluid particles and boundary virtual particles are close to or away from the boundary normal, the density of fluid particles increases or decreases accordingly. Then, the pressure of the virtual particles is obtained by interpolating the fluid pressure field. When the fluid is relatively close to the boundary, the virtual particles exert repulsion on the fluid particles through the pressure gradient in order to prevent fluid particles from crossing the   4 Geofluids boundary. The boundary equation of the virtual particle method is shown in formula (8): In formula (8), ðdρ i /dtÞ b , ðdv i /dtÞ b are the density and velocity increment, respectively, of the fluid particle i caused by the boundary virtual particle b, and n b represents the unit normal vector of the virtual particle b.

Establishment of the Calculation Model
As the impact force of debris flow comprises the impact force of the debris flow solid-liquid two-phase flow and the interaction force between them, this paper constructs the interaction model between the two-phase flow and retaining dam by taking Nanjiaogou in the western mountain area of Fangshan District of Beijing as a case study. Fangshan District is located in the Taihang mountain system, and the mountain trend is roughly Northeast. The image map and floor plan of the Nanjiaogou small watershed are shown in Figure 2. The height difference of the Nanjiaogou small watershed is 1.12 km, the watershed area is 24.3 km2, and the slope of the main gully is 19°-47°, which is a large debris flow gully. The longitudinal profile of the Nanjiaogou main gully is shown in Figure 3. The average debris flow source particle density in the gully watershed is 2650 kg/m 3 , the average friction angle is 23°, and the equivalent diameter of the source particles in the gully is 0.2 cm-0.4 m.
In order to ensure the appearance of "dragon head, dragon body and dragon tail", the retaining dam is set at a certain distance. The debris flow moves along the debris flow channel under the action of gravity. The average equivalent diameter of provenance particles is 0.3 m, and the slope of the formation area of the main gully is 30°.

Geometric
Model. The geometric model of debris flow and retaining dam is shown in Figure 4, which is composed of debris flow solid particles, debris flow slurry body, debris flow channel, and retaining dam. The specific geometric model parameters are shown in Table 1. Since the topographic conditions of the site are complex and this paper aims to analyze the dynamic response of the retaining dam, the numerical model is simplified as follows: The bottom and side of the debris flow gully are simplified to a plane, and the average specific drop of the debris flow gully is set to 30%; the debris flow is composed of solid particles and slurry and is given the initial velocity at the same time as the large rock to impact the retaining dam.

Material Model.
For debris flow gullies, the use of rigid shell materials for simulation can significantly shorten the simulation time. A total of 2916 shell elements are used in this paper. The density of debris flow gully is 3000 kg, Poisson's ratio is 0.24, and elasticity modulus is 30GPa.
The measured average density, friction angle, and Poisson's ratio of solid particles in debris flow provenance are 2650 kg/m 3 , 23°, and 0.3, respectively. Discrete element particle simulation is used to simulate a total of 2172 particles. The slurry density of debris flow is 1000 kg/m 3 , and the elastoplastic hydrodynamic model is used as the constitutive model. SPH element is used to simulate the constitutive model, and Us-Up equation of state is used to simulate the     In view of the surface roughness of debris flow solid particles, the friction coefficient between debris flow solid particles is 0.3, the friction coefficient between debris flow solid particles and gully and retaining dam is 0.3, and the friction factor between debris flow liquid and gully and retaining dam is 0.12.
The debris flow slurry and solid particles are given the gravity acceleration of 9.8 m/s 2 , and the solid-liquid phase of debris flow moves along the debris flow gully under the action of gravity.
The calculation time is 10s, the calculation has a total of 1171895 incremental steps, and the CPU time is 28 h 15 min 50s. The process simulates the complete impact process of debris flow. particles-liquid debris flow slurry jointly involved in the interaction and began to impact the retaining dam. At the same time, the solid-liquid two-phase debris flow "dragon body" impacted the retaining dam and began to climb; at t = 2:9 s, the "dragon body" impacted the retaining dam, and most of the debris flow was blocked, and the debris flow continued to impact the retaining dam. The splashing and diffusion of fluid particles can be prevented using the SPH method, showing a good visual effect of detailed simulation of solid-liquid interaction; at t = 5:2 s, the solid-liquid twophase debris flow silted back and finally formed a static load at the bottom of the retaining dam.

Climbing and Silting
Process. The climbing process of the solid-liquid two-phase debris flow impacting the barrier dam is shown in Figure 6. At t = 2:1 s, the debris flow reached the bottom of the retaining dam, impacted the retaining dam and began to climb, and the climbing height was 0.6 m at this point; at t = 2:3 s, the debris flow continued to climb, and local siltation occurs. At this point, the splashing and diffusion of SPH fluid particles can be observed, and the climbing height was 2.2 m; at t = 2:9 s, the "dragon body" impacted the retaining dam, and the debris flow climbed to the highest, resulting in obvious siltation and particle splash, and the highest climbing height was 4.5 m; at t = 5:2 s, the solidliquid two-phase debris flow silted and finally formed a static load at the bottom of the retaining dam, and the siltation height was the height of the retaining dam.  Figure 7. The results showed that when the slope was 50°, at t = 2:6 s, the displacement in the middle position of the retaining dam crest reached the maximum value of 5.61 mm and decreased rapidly to 3.04 mm after reaching the peak, and the displacement value was the peak value of the displacement in the middle position of the V, Magnitude +1.158e+01 +1.062e+01 +9.650e+00 +8.685e+00 +7.720e+00 +6.755e+00 +5.790e+00 +4.825e+00 +3.860e+00 +2.895e+00 +1.930e+00 +9.650e-01 +0.000e+00    8 Geofluids retaining dam crest when the slope was about 20°. When the slope was 30°and 40°, the peak displacement of the middle position of the retaining dam crest was 2.75 mm and 3.55 mm, respectively, and then decreased to the static displacement, which was about the peak displacement when the slope was 20°. The simulation results are consistent with the test results of Hungr [21], which shows that the method can be applied to simulate the impact of debris flowretaining dam.

Time-History
Analysis of the Impact Force. Figure 8 shows the impact time history of the dam with different slopes (20°, 30°, 40°, and 50°). As shown in this figure, when the slope was 50°, at t = 2:6 s, the impact force of the solid-liquid two-phase debris flow on the retaining dam reached its peak value, and the peak impact force was 9896 kN and then decreased rapidly to 3032 kN.
When the slope was about 20°, the impact force of the retaining dam reached its peak value. When the slope was 30°and 40°, the impact force peak value of the retaining dam was 5626 kN and 6917 kN, respectively, and then decreased to the static load, and the static load value was about the displacement peak value when the slope was 20°.
In order to verify the rationality of the results, the numerical results are compared according to Kwan's [22,23] semiempirical formula for calculating the peak impact force of debris flow based on hydrodynamic theory.
where the dynamic pressure coefficient is κ = aFr b ; ρ is the density of the solid-liquid two-phase flow; v is the impact V, Magnitude +1.158e+01 +1.062e+01 +9.650e+00 +8.685e+00 +7.720e+00 +6.755e+00 +5.790e+00 +4.825e+00 +3.860e+00 +2.895e+00 +1.930e+00 +9.650e-01 +0.000e+00  The numerical calculation methods and existing empirical formula calculation results in this paper are shown in Table 2, with a calculation error between 7.8% and 16.2%. This is because the collision between solid particles and the coupling between solid particles and debris flow slurry dissipate a certain amount of energy. The impact force calculated by the numerical calculation method in this paper is close to that of the empirical value, which further illustrates the rationality of the numerical calculation method in this paper. 4.5. Impact Force Distribution. The distribution of the impact force of the solid-liquid two-phase debris flow impact dam is shown in Figure 9. At t = 2:1 s, when the debris flow reached the bottom of the retaining dam and began to impact the retaining dam, the impact force was mainly distributed in the middle position at the bottom of the retaining dam; at t = 2:5 s, the debris flow continued to climb, and the impact force of the debris flow was mainly distributed in the middle and lower position of the retaining dam and spread rapidly to the upper part of the retaining dam; at t = 2:9 s, the "dragon body" impacted the retaining dam, and the impact force of debris flow was distributed in the whole retaining dam, which was consistent with what we observed during the impact process; at t = 5:2 s, the solid-liquid two-phase debris flow silted back and finally formed a static load at the bottom of the retaining dam. At this point, the impact force was distributed in the whole retaining dam, and due to the action of particles, the impact force distribution showed a point-like distribution, which could not be shown if the method of fluid mechanics is used.
The results show that the coupled SPH-DEM-FEM method can vividly simulate the impact behavior of the solid-liquid two-phase debris flow on the retaining dam, and reproduce impact, climbing, and siltation in the process of debris flow impact; the displacement time-history analysis of the retaining dam indicates that the simulation results are consistent with the existing research results. The impact force time-history analysis of the retaining dam is in the same order of magnitude as that of the existing empirical formula, and the results are basically credible. The climbing and impact force distribution of debris flow is consistent with the observed impact process, which shows that this method can be used to simulate the solid-liquid two-phase debris flow, and the analysis results have a certain reference value for the engineering design of the debris flow impactretaining dam.

Conclusion
(1) Based on the coupled SPH-DEM-FEM numerical method, this paper analyzes the dynamic interaction of the solid debris flow particle-liquid phase debris flow slurry-retaining dam and reproduces the whole impact, climbing, and siltation process of the debris flow impact-retaining dam  10 Geofluids (2) The dynamic response analysis of the debris flow impact-retaining dam under different slopes is compared. The dynamic response law can provide some reference for the design of the debris flow-retaining dam (3) The variation law of the impact force of debris flow under different slopes is compared, and it is found by comparison that the impact force is close to that of the empirical formula. The numerical calculation method can provide some reference for the calculation of the debris flow impact force and the design of the further retaining dam (4) The impact force distribution of the whole process of the solid-liquid two-phase debris flow impactretaining dam is studied. The impact force distribution has a certain reference value for the engineering design of the debris flow-retaining dam (d) t = 5:2 s Figure 9: Impact force distribution.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.