# SPH Simulation of Structures Impacted by Tailing Debris Flow and Its Application to the Buffering Effect Analysis of Debris Checking Dams

Authors: Qing-Yun Zeng; Jian-Ping Pan; Han-Zheng Sun
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Since the tailing dam fails catastrophically with substantial instantaneous deformation, it is difficult to measure the migration of debris flow caused by the failure of the tailings dam. A simulation model of tailing debris flow based on Smoothed Particle Hydrodynamics (SPH) theory of elastic-plastic constitutive equation has been established by considering the viscoplasticity of mud and the elastic-plastic characteristics of tailing sand to investigate the impact effect of tailing flow on the downstream structures. By comparing the experimental and two different simulation results obtained, it can be concluded that SPH elastic-plastic constitutive model can effectively simulate the accumulation and migration processes of the tailing debris flow, which indicates that the SPH model has good applicability to solve geotechnical large deformation problems of similar tailings flow slide. Then, the verified simulation model developed based on a series of simulations of tailing debris flow propagations was used to determine the momentum reduction on the downstream structure resulting from the presence of a simple checking dam perpendicular to the direction of propagation and to determine the characteristics of stresses applied to this structure in terms of peak impact force and evolution over time to the main flow direction.


---

## Body

## 1. Introduction

Tailings dam is one of the most complex geotechnical structures. As the rainfall increases, the tailing infiltration line in the reservoir increases, and once the saturation is reached, dam break occurs in the form of high-speed moving sand-bearing mudflow [1, 2]. These flows commonly occur without a prior warning in mountainous regions around the world and therefore cause massive loss of both human lives and properties. Debris flow caused due to the collapse of saturated tailing soil is usually difficult to solve in fluid dynamics. The relevant parameters are not only difficult to be measured in the field, but the calculation is also pretty complicated. Therefore, the numerical simulation calculation has become the first choice to study the tailing debris flow. The literature shows that most of the existing simulation methods applied to viscous debris flow are primarily grid-based finite element method (FEM) [3] and Boltzmann algorithm [4]. FEM is bounded by the element mesh to deal with large soil deformation. FEM has cell grid constraints that either interrupt the calculation or provide the inaccurate result of the calculation [5]. Similarly, the Boltzmann algorithm often requires a perfect grid to solve the higher-order differential equations, which greatly increases the computational cost [6].Nowadays, the meshless method has been developed to overcome the deficiencies of the above-stated calculation method applied in the mechanical behavior of geomaterial analysis. The whole solution region in the meshless method is discretized into independent nodes and has no dependence on the grid, and the calculation methods under the different framework of the meshless theory are established, such as DEM (discrete element method), EFG (element-free galerkin method, improved DEM), and SPH (smoothed particle hydrodynamics method), especially the DEM and EFG methods show a highly computational accuracy in the application of rock displacement and deformation in rock mechanics or geotechnical mechanics [7–10]. However, the applicability of methods such as DEM and EFG in large deformation of soil flow problems of similar debris flow or landslide have some defects, because these methods are not easy to determine some nonphysical parameters such as the friction and elastic parameters of the contact zone which are to be calibrated as they are not directly obtained from experiments and they are difficult to directly describe the stress-strain relationship of soil, which limited their application in this research direction [11–13]. The saturated tailing soil mainly experiences the stage of soil stiffness degradation and the stage of tailing debris flow movement in the process of large deformation of soil flow. When the saturated soil has been transformed from solid to liquid, it will be difficult reasonably to calculate and predict the complex dynamic behavior of soil by using traditional solid mechanics [13]. Therefore, it is necessary to introduce elastic-plastic mechanics and fluid dynamics to analyze the different phase states in the process of saturated tailing soil deformation.Recently, a meshless particle method known as SPH has attracted more attention. The SPH method can get rid of the finite element method's dependence on the grid and shows obvious advantages in the problems related to grid distortion and grid movement [14, 15]. On the other hand, in dealing with the problems related to dynamic response which occurs during large deformation of debris flow, SPH can easily introduce constitutive equations describing the properties of various materials including solid mechanics and hydrodynamics, and the advantages of Lagrangian particularization will describe the moving state of a solid-fluid mixture in the form of particles, which is beneficial to accurately grasp the deformation characteristics of different positions in the process of large deformation of debris flow [16, 17]. Then, several researchers had applied SPH to the migration impact simulation of viscous debris flow. Antonio et al. [18] have successfully established the simulation calculation model of viscous debris flow by referring to the test parameters of indoor debris flow sliding experiments conducted by Laigle and Coussot [19]. Huang et al. [20] utilized SPH to carry out the migration rule analysis of landslide and debris flow disaster under the impact of the Wenchuan earthquake. Dai et al. [21] established the coupled SPH model to simulate the mountain debris flow-structure interaction. However, the SPH method is not widely used in the case of non-Newtonian flow such as tailing debris flow. It can be concluded from the related literature about the SPH control equation that the characteristics of the stress-strain relation of tailings debris flow are rarely reported. As a result, many researchers pay attention to the research and application of debris flow constitutive model based on SPH theory.The present paper focused on the application of the SPH method to the simulations of structures impacted by tailing debris flow that involves large deformations of soil flow. Firstly, the main features of the SPH model were presented with specific attention considered to the treatment of boundary conditions and rheological behavior and SPH fluid-solid coupled model based on the constitution of elastic-plastic constitutive mechanics for tailings debris flow in Section2 was established. Then, in Section 3, the model is validated by simulating the impact effect of tailing flow on the downstream structures with the indoor model test and two different numerical simulations (SPH and FEM) to check the applicability of SPH in simulating the tailing debris flow. In Section 4, the verified SPH simulation model has been applied to derive the respective relationship between the height of checking dam and height of mud surface, impact force, and momentum reduction of the tailing debris flow that occurs at the downstream structure. The results obtained from the calculations were then introduced and discussed in detail. Finally, several conclusions were summarized in Section 5.


## 2. SPH Formulations of Structures Impacted by Tailing Debris Flow

To get rid of the grid constraint caused by the mesh division in the FEM model, the SPH method uses the free distribution of discrete particles that carry field variables instead of its problem domain. The aforementioned method also introduces the kernel function to simplify the complex interactive relationship between particles [22]. Therefore, tailing debris flow problems can be solved by tracking the movements of the particles and the changes of the carried field variables.


### 2.1. SPH Integral Interpolation Function Approximation

For continuous fluid motion, the SPH integral interpolation functionf (x) can be accurately simulated by a large number of particles containing mass and momentum [23]. The function f (x) at position x (x, y, z) is approximated by the following expression:(1)fx=∫Ωfx′Wx−x′,hdx′,where Ω represents the computational domain at position x, dx′ is the element of volume over the domain, W is the weight equation or smooth kernel function, and h is smooth length, which determines the size of the computational domain of the kernel function W.Since the distribution of the initial interpolation points is not random, it changes constantly with the motion of the fluid after the initial value is taken. To describe the fluid motion more accurately, (1) is discretized into particle form, which is given by the following equation:(2)fxi=∑j=1NmjρjfxjWij,where i represents the concerning particle, j is a neighboring particle in Ω, and mj and ρj represent the mass and density of particle j, respectively. Smooth length h is 1.05–1.5 times of particle spacing l0, which determines the number of influenced particles in the computational domain (Figure 1).Figure 1
The computational domain Ω of particlei.The three-dimensional form of cubic spline kernel functionWij=WR,h is adopted in the following form [16], where relative displacement vector rij=xi−xj, rij=xi−xj=l0 and R=rij/h:(3)WR,h=32πh3×23−R2+12R3,0≤R<1,162−R3,1≤R<2,0,R≥2.The smooth lengthh should satisfy two requirements [24]. First, the value of h should be smaller than the influence radius of the individual particles, so that less number of iterations is required for interpolation calculation. Second, h must be larger than the particle spacing l0 to ensure the existence of correlation particles in the range of influence; otherwise, the lack of correlation between particles might lead to the particle penetration of the finite element mesh in the coupling calculation of SPH and FEM. The purpose of applying the smooth length is to ensure that the total mass of the nearest neighbor particles is constant or maintained within an allowable range. According to the discussion on the relationship between smooth length and particle spacing, the rational range of smooth length h should be l0<h≤1.5l0 [25]. Based on these requirements and the cost of calculation, the economical and reasonable value of h=1.05l0 has been selected for the present research work.


### 2.2. Governing Equations of Viscous Debris Flow

When the tailing dam runs at a high water level, the tailings are usually close to saturation. Since the shear strength of the tailings dam is very low, the sliding surface appears easily, which in turn forms the dam breaking debris flow. The tailings debris flow satisfies the characteristic description of high sediment flow in terms of mechanical properties and movement mechanism, showing obvious incompressible viscous non-Newtonian flow characteristics [26]. The most widely used non-Newtonian models for debris flow simulations are the Bingham [27, 28] and Herschel-Bulkley [18, 29] models.In the present work, referring to the study of rheological control of mountain debris flow [20, 30, 31], the viscosity between tailing particles in high-speed motion is very small in comparison to the gravity of the tailing fluid and the viscosity resistance between slurry; therefore, it can be neglected. At the same time, the mixture of tailing particles and rainwater has been considered as a uniform continuous incompressible Bingham fluid. The governing equations of viscous debris flow for the Bingham model are the continuity equation and momentum equation. The following equation expresses the application of the aforementioned equations (continuity and momentum equation) in the present research work:(4)dρdt=−ρ∇⋅u,ρdudt+u⋅∇u=−∇P+∇⋅μ∇u+∇uT+ρg+σ=κ∇a,where u is the velocity vector, ∇ represents a divergence symbol, ρ is the density of debris flow, P is hydrostatic pressure, σ= is the total stress tensor, κ is the interface curvature, μ represents the viscosity coefficient, g is the gravitational acceleration vector, and a is the distance function of the contact surface between air and mud.To maintain a good correlation between the governing equations and the SPH kernel function, governing equations are written in the following Lagrange form [23, 32]:(5)dρidt=∑j=1Nmjuijα∂Wij∂xiα,(6)duiαdt=−∑j=1NmjPiρi2+Pjρj2+Πij∂Wij∂xiα+gα+∑j=1Nmjσiαβ+τiαβρi2+σjαβ+τjαβρj2∂Wij∂xiα,where gα is the gravitational acceleration, u is the velocity of debris flow, σ is the principal stress tensor, τ is the shear stress tensor, Pi and Pj represent the pressure on particles i and j, respectively, and α and β denote the Cartesian components x, y, and z with the Einstein convention applied to repeated indices.To improve the stability of numerical calculations and avoid the shock of solid-liquid coupling contact of dam break fluid, an artificial viscosityΠij is added to the momentum equation which denotes dissipations of the kinetic energy of the shocked region into heat. Monaghan and Gingold [33] adopted this approach and introduced it in the following form:(7)Πij=−q1c¯ijφij+q2φij2ρ¯ij,uij⋅rij<0,0,uij⋅rij≥0,ϕij=huijrijrij2+0.01h2,where rij=ri−rj is the relative distance, uij=ui−uj is the relative velocity, and c¯ij=ci+cj/2 and ρ¯ij=ρi+ρj/2 represent the mean sound velocity and average density, respectively. To avoid the overflow of the calculated results, the parameters ϕ2=0.01h2 are satisfied under the condition that the smooth length h=h¯ij=hi+hj/2 is dynamic. For the treatment of dam, break debris flows under viscoelastic characteristics, q1 and q2 are considered to be constant. In the current work, we have assumed q2 = 0.01–0.1 and q2 = 0 [34].


### 2.3. Constitutive Equations of Tailing Debris Flow

Debris flow, although possesses a complex composition, can be broadly considered as a mixture of water, clay, and granular matter. To simplify the calculation, the stress-strain relation of viscoplastic fluids such as the tailing debris flow is assumed to follow the Bingham model [21, 27, 28].According to the characteristics of Bingham flow, the liquid-phase stress-strain relation of tailing debris flow is described by the following expression:(8)ταβ=μ+τyεΠεαβ,where ταβ is the shear stress tensor, μ is the viscosity coefficient, τy represents the yield strength, εαβ is the shear strain tensor, ε∏ is the second invariant of the tensor of shear strain, and α and β denote the Cartesian components x, y, and z with the Einstein convention applied to repeated indices.For ideal viscoplastic materials, the shear strain rateε˙αβ is given by the following equation:(9)ε˙αβ=∂εαβ∂t=12∂uα∂xβ+∂uβ∂xα,uα=ux,uy,uzT,uβ=uy,uz,uxT,xα=x,y,zT,xβ=y,z,xT,where ε˙iαβ of the particle i in the SPH form can be expressed as the following equation [35]:(10)ε˙iαβ=∂εiαβ∂t=12∑j=1Nmjρjujα−uiα∂Wij∂xiβ+∑j=1Nmjρjujβ−uiβ∂Wij∂xiα.Considering the solid-liquid cohesive effect of saturated tailings, the molar Coulomb yield criterion is introduced to solve the shear stress:(11)τy=Ptanφ+C,where P is the hydrostatic pressure, ϕ is the internal friction angle, and C is the cohesive strength.For the solid-phase part of debris flow, the stress-strain relation of elastic-plastic material in solid mechanics principle is given by the following equation:(12)σαβ=−Pδαβ+sαβ,where the rate of partial shear stress tensor sαβ can be modified according to Jaumann criterion [12] which is as follows:(13)∂siαβ∂t=2Gεiαβ−13ε˙iγγ−siαγω˙iβγ+ω˙iαγsiβγ,where G is the shear modulus of tailings and ω˙iαγω˙iβγ represents the spin rate tensor. For ideal elastoplastic materials, the spin rate tensor must satisfy the following condition:(14)ω˙αβ=∂ωαβ∂t=12∂uα∂xβ−∂uβ∂xα=12∇⋅u,ω˙αγ=12∂uα∂xγ−∂uγ∂xα,ω˙βγ=12∂uβ∂xγ−∂uγ∂xβ,where ω˙iαγ and ω˙iβγ of the particle i in SPH form can be written as follows [35]:(15)ω˙iαγ=∂ωiαγ∂t=12∑j=1Nmjρjujα−uiα∂Wij∂xiγ−∑j=1Nmjρjujγ−uiγ∂Wij∂xiα,ω˙iβγ=∂ωiβγ∂t=12∑j=1Nmjρjujβ−uiβ∂Wij∂xiγ−∑j=1Nmjρjujγ−uiγ∂Wij∂xiβ.According to the variation of particle unit velocity in fluid mechanics [36], the relative velocity between two particles can be expressed as a matrix:(16)duα=duxduyduz=∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂zdxα=∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂z=dxdydz,where the square matrix in the upper formula can be decomposed into(17)∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂z=012∂ux∂y−∂uy∂x12∂ux∂z−∂uz∂x12∂uy∂x−∂ux∂y012∂uy∂z−∂uz∂y12∂uz∂x−∂ux∂z12∂uz∂y−∂uy∂z0+∂ux∂x12∂ux∂y+∂uy∂x12∂ux∂z+∂uz∂x12∂ux∂y+∂uy∂x∂uy∂y12∂uy∂z+∂uz∂y12∂ux∂z+∂uz∂x12∂uy∂z+∂uz∂y∂uz∂z=0ω˙xyω˙xz−ω˙xy0ω˙yz−ω˙xz−ω˙yz0+ε˙xxε˙xyε˙xzε˙xyε˙yyε˙yzε˙xzε˙yzε˙zz=0ω˙xyω˙xz−ω˙xy0ω˙yz−ω˙xz−ω˙yz0+ε˙xx000ε˙yy000ε˙zz+0ε˙xyε˙xzε˙xy0ε˙yzε˙xzε˙yz0=ω˙αβ+ε˙γγ+ε˙αβ.The first matrix in the upper formula (17) is ω˙αβ (antisymmetric), and ε˙αβ and ε˙γγ matrices are symmetric. These three matrices are also called second-order tensors in hydrodynamics [36].Among the nine components in the antisymmetric matrixω˙αβ, there are only three independent components which are as follows:(18)ω˙xy=12∂ux∂y−∂uy∂x,ω˙xz=12∂ux∂z−∂uz∂x,ω˙yz=12∂uy∂z−∂uz∂y.These three components are exactly the three components of the rotational angular velocity vector of the fluid particle, soω˙αβ is called the rotational velocity tensor. At the same time, ω˙αβ=ω˙xyi+ω˙yzj+ω˙xzk, that is, half the curl of the velocity vector is satisfied (equation (14)).Among the nine components of symmetric matricesε˙αβ and ε˙γγ, there are only six independent components as follows:(19)ε˙xx=∂ux∂x,ε˙yy=∂uy∂y,ε˙zz=∂uz∂z,ε˙xy=ε˙yx=12∂ux∂y+∂uy∂x,ε˙yz=ε˙zy=12∂uz∂y+∂uy∂z,ε˙xz=ε˙zx=12∂ux∂z+∂uz∂x,then ε˙γγ=ε˙xxi+ε˙yyj+ε˙zzk and ε˙αβ=ε˙xyi+ε˙yzj+ε˙xzk:(20)ε˙γγ=∂uα∂xα,ε˙αβ=∂εαβ∂t=12∂uα∂xβ+∂uβ∂xα,where ε˙γγ is the linear strain rate tensor of fluid particle in fluid mechanics in the direction of three axes and ε˙αβ is the fluid particle’s angular deformation rate tensor. Therefore, the tensor ε˙αβ+ε˙γγ is called the strain rate tensor in fluid mechanics, and ω˙αβ+ε˙αβ+ε˙γγ is called the velocity gradient tensor.Moreover, to ensure the weak compressibility and continuity of tailing debris flow, changes in hydrostatic pressurePi within the particle motion domain can be described directly by the equation of state (EOS) in the form of particle density. The SPH equation of state in the form of weak compression was proposed by Monaghan [37], which is as follows:(21)Pi=c02ρ0λρiρ0λ−1,where constant λ=7, ρ0 is the initial average density of tailing debris flow, and c0 is the sound velocity in the corresponding fluid medium. In the current work, c0 = 340 m/s. In the numerical model, it is assumed that the tailing is elastoplastic materials and the mud is weakly compressible viscoplastic fluid. The elastic-plastic constitutive model of the tailings flow is finally established by using the equation of state (21) coupled with (8) and (12) to determine the contact correlation between the mud and the tailings (Figure 2).Figure 2
Modeling of tailing-mud interaction in SPH simulation.


### 2.4. Boundary Condition Used in Current SPH Model

While applying SPH, to save computation time, the boundary was assumed to be pressure-free or constituted of weak-pressure interfaces. However, the idea of “node-to-surface” unidirectional contact control based on compensation by introducing the gradient of momentum and mass into the interpolation kernel function was chosen as it becomes the core idea of dealing with boundary problems [38]. Meanwhile, to avoid particle penetration at the boundary, an improved virtual particle compensation method was introduced. The simulation accuracy was improved substantially by setting virtual particles with no mass and no volume at the boundary to avoid the interference of real boundary particles and the moving fluid particles [39, 40].While applying the virtual particle compensation method, a discrete layer of particles is set at the boundary where the displacement and mass of the boundary particles are assumed to be zero. According to Newton's third law, when the boundary particles come in contact with the fluid particles, repulsive forces are produced in the vertical direction. Consequently, the velocity vector and the displacement vector of the fluid particles change, and the fluid particles are separated and they do not penetrate, as presented in Figure3. This force contribution from the boundary particles is added to the right-hand side of the momentum equation ((6)).Figure 3
Boundary particle distribution in the fictitious particle method.


### 2.5. Operational Framework of SPH Model

The SPH governing equations of the final elastoplastic constitutive model are obtained by substituting the SPH form variable ((8) and ((12)) from Section 2.3 into the momentum (6) and combining with continuity (5). According to the SPH governing equations, the flow chart of the SPH program is shown in Figure 4.Figure 4
Flow chart of SPH program.


## 2.1. SPH Integral Interpolation Function Approximation

For continuous fluid motion, the SPH integral interpolation functionf (x) can be accurately simulated by a large number of particles containing mass and momentum [23]. The function f (x) at position x (x, y, z) is approximated by the following expression:(1)fx=∫Ωfx′Wx−x′,hdx′,where Ω represents the computational domain at position x, dx′ is the element of volume over the domain, W is the weight equation or smooth kernel function, and h is smooth length, which determines the size of the computational domain of the kernel function W.Since the distribution of the initial interpolation points is not random, it changes constantly with the motion of the fluid after the initial value is taken. To describe the fluid motion more accurately, (1) is discretized into particle form, which is given by the following equation:(2)fxi=∑j=1NmjρjfxjWij,where i represents the concerning particle, j is a neighboring particle in Ω, and mj and ρj represent the mass and density of particle j, respectively. Smooth length h is 1.05–1.5 times of particle spacing l0, which determines the number of influenced particles in the computational domain (Figure 1).Figure 1
The computational domain Ω of particlei.The three-dimensional form of cubic spline kernel functionWij=WR,h is adopted in the following form [16], where relative displacement vector rij=xi−xj, rij=xi−xj=l0 and R=rij/h:(3)WR,h=32πh3×23−R2+12R3,0≤R<1,162−R3,1≤R<2,0,R≥2.The smooth lengthh should satisfy two requirements [24]. First, the value of h should be smaller than the influence radius of the individual particles, so that less number of iterations is required for interpolation calculation. Second, h must be larger than the particle spacing l0 to ensure the existence of correlation particles in the range of influence; otherwise, the lack of correlation between particles might lead to the particle penetration of the finite element mesh in the coupling calculation of SPH and FEM. The purpose of applying the smooth length is to ensure that the total mass of the nearest neighbor particles is constant or maintained within an allowable range. According to the discussion on the relationship between smooth length and particle spacing, the rational range of smooth length h should be l0<h≤1.5l0 [25]. Based on these requirements and the cost of calculation, the economical and reasonable value of h=1.05l0 has been selected for the present research work.


## 2.2. Governing Equations of Viscous Debris Flow

When the tailing dam runs at a high water level, the tailings are usually close to saturation. Since the shear strength of the tailings dam is very low, the sliding surface appears easily, which in turn forms the dam breaking debris flow. The tailings debris flow satisfies the characteristic description of high sediment flow in terms of mechanical properties and movement mechanism, showing obvious incompressible viscous non-Newtonian flow characteristics [26]. The most widely used non-Newtonian models for debris flow simulations are the Bingham [27, 28] and Herschel-Bulkley [18, 29] models.In the present work, referring to the study of rheological control of mountain debris flow [20, 30, 31], the viscosity between tailing particles in high-speed motion is very small in comparison to the gravity of the tailing fluid and the viscosity resistance between slurry; therefore, it can be neglected. At the same time, the mixture of tailing particles and rainwater has been considered as a uniform continuous incompressible Bingham fluid. The governing equations of viscous debris flow for the Bingham model are the continuity equation and momentum equation. The following equation expresses the application of the aforementioned equations (continuity and momentum equation) in the present research work:(4)dρdt=−ρ∇⋅u,ρdudt+u⋅∇u=−∇P+∇⋅μ∇u+∇uT+ρg+σ=κ∇a,where u is the velocity vector, ∇ represents a divergence symbol, ρ is the density of debris flow, P is hydrostatic pressure, σ= is the total stress tensor, κ is the interface curvature, μ represents the viscosity coefficient, g is the gravitational acceleration vector, and a is the distance function of the contact surface between air and mud.To maintain a good correlation between the governing equations and the SPH kernel function, governing equations are written in the following Lagrange form [23, 32]:(5)dρidt=∑j=1Nmjuijα∂Wij∂xiα,(6)duiαdt=−∑j=1NmjPiρi2+Pjρj2+Πij∂Wij∂xiα+gα+∑j=1Nmjσiαβ+τiαβρi2+σjαβ+τjαβρj2∂Wij∂xiα,where gα is the gravitational acceleration, u is the velocity of debris flow, σ is the principal stress tensor, τ is the shear stress tensor, Pi and Pj represent the pressure on particles i and j, respectively, and α and β denote the Cartesian components x, y, and z with the Einstein convention applied to repeated indices.To improve the stability of numerical calculations and avoid the shock of solid-liquid coupling contact of dam break fluid, an artificial viscosityΠij is added to the momentum equation which denotes dissipations of the kinetic energy of the shocked region into heat. Monaghan and Gingold [33] adopted this approach and introduced it in the following form:(7)Πij=−q1c¯ijφij+q2φij2ρ¯ij,uij⋅rij<0,0,uij⋅rij≥0,ϕij=huijrijrij2+0.01h2,where rij=ri−rj is the relative distance, uij=ui−uj is the relative velocity, and c¯ij=ci+cj/2 and ρ¯ij=ρi+ρj/2 represent the mean sound velocity and average density, respectively. To avoid the overflow of the calculated results, the parameters ϕ2=0.01h2 are satisfied under the condition that the smooth length h=h¯ij=hi+hj/2 is dynamic. For the treatment of dam, break debris flows under viscoelastic characteristics, q1 and q2 are considered to be constant. In the current work, we have assumed q2 = 0.01–0.1 and q2 = 0 [34].


## 2.3. Constitutive Equations of Tailing Debris Flow

Debris flow, although possesses a complex composition, can be broadly considered as a mixture of water, clay, and granular matter. To simplify the calculation, the stress-strain relation of viscoplastic fluids such as the tailing debris flow is assumed to follow the Bingham model [21, 27, 28].According to the characteristics of Bingham flow, the liquid-phase stress-strain relation of tailing debris flow is described by the following expression:(8)ταβ=μ+τyεΠεαβ,where ταβ is the shear stress tensor, μ is the viscosity coefficient, τy represents the yield strength, εαβ is the shear strain tensor, ε∏ is the second invariant of the tensor of shear strain, and α and β denote the Cartesian components x, y, and z with the Einstein convention applied to repeated indices.For ideal viscoplastic materials, the shear strain rateε˙αβ is given by the following equation:(9)ε˙αβ=∂εαβ∂t=12∂uα∂xβ+∂uβ∂xα,uα=ux,uy,uzT,uβ=uy,uz,uxT,xα=x,y,zT,xβ=y,z,xT,where ε˙iαβ of the particle i in the SPH form can be expressed as the following equation [35]:(10)ε˙iαβ=∂εiαβ∂t=12∑j=1Nmjρjujα−uiα∂Wij∂xiβ+∑j=1Nmjρjujβ−uiβ∂Wij∂xiα.Considering the solid-liquid cohesive effect of saturated tailings, the molar Coulomb yield criterion is introduced to solve the shear stress:(11)τy=Ptanφ+C,where P is the hydrostatic pressure, ϕ is the internal friction angle, and C is the cohesive strength.For the solid-phase part of debris flow, the stress-strain relation of elastic-plastic material in solid mechanics principle is given by the following equation:(12)σαβ=−Pδαβ+sαβ,where the rate of partial shear stress tensor sαβ can be modified according to Jaumann criterion [12] which is as follows:(13)∂siαβ∂t=2Gεiαβ−13ε˙iγγ−siαγω˙iβγ+ω˙iαγsiβγ,where G is the shear modulus of tailings and ω˙iαγω˙iβγ represents the spin rate tensor. For ideal elastoplastic materials, the spin rate tensor must satisfy the following condition:(14)ω˙αβ=∂ωαβ∂t=12∂uα∂xβ−∂uβ∂xα=12∇⋅u,ω˙αγ=12∂uα∂xγ−∂uγ∂xα,ω˙βγ=12∂uβ∂xγ−∂uγ∂xβ,where ω˙iαγ and ω˙iβγ of the particle i in SPH form can be written as follows [35]:(15)ω˙iαγ=∂ωiαγ∂t=12∑j=1Nmjρjujα−uiα∂Wij∂xiγ−∑j=1Nmjρjujγ−uiγ∂Wij∂xiα,ω˙iβγ=∂ωiβγ∂t=12∑j=1Nmjρjujβ−uiβ∂Wij∂xiγ−∑j=1Nmjρjujγ−uiγ∂Wij∂xiβ.According to the variation of particle unit velocity in fluid mechanics [36], the relative velocity between two particles can be expressed as a matrix:(16)duα=duxduyduz=∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂zdxα=∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂z=dxdydz,where the square matrix in the upper formula can be decomposed into(17)∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂z=012∂ux∂y−∂uy∂x12∂ux∂z−∂uz∂x12∂uy∂x−∂ux∂y012∂uy∂z−∂uz∂y12∂uz∂x−∂ux∂z12∂uz∂y−∂uy∂z0+∂ux∂x12∂ux∂y+∂uy∂x12∂ux∂z+∂uz∂x12∂ux∂y+∂uy∂x∂uy∂y12∂uy∂z+∂uz∂y12∂ux∂z+∂uz∂x12∂uy∂z+∂uz∂y∂uz∂z=0ω˙xyω˙xz−ω˙xy0ω˙yz−ω˙xz−ω˙yz0+ε˙xxε˙xyε˙xzε˙xyε˙yyε˙yzε˙xzε˙yzε˙zz=0ω˙xyω˙xz−ω˙xy0ω˙yz−ω˙xz−ω˙yz0+ε˙xx000ε˙yy000ε˙zz+0ε˙xyε˙xzε˙xy0ε˙yzε˙xzε˙yz0=ω˙αβ+ε˙γγ+ε˙αβ.The first matrix in the upper formula (17) is ω˙αβ (antisymmetric), and ε˙αβ and ε˙γγ matrices are symmetric. These three matrices are also called second-order tensors in hydrodynamics [36].Among the nine components in the antisymmetric matrixω˙αβ, there are only three independent components which are as follows:(18)ω˙xy=12∂ux∂y−∂uy∂x,ω˙xz=12∂ux∂z−∂uz∂x,ω˙yz=12∂uy∂z−∂uz∂y.These three components are exactly the three components of the rotational angular velocity vector of the fluid particle, soω˙αβ is called the rotational velocity tensor. At the same time, ω˙αβ=ω˙xyi+ω˙yzj+ω˙xzk, that is, half the curl of the velocity vector is satisfied (equation (14)).Among the nine components of symmetric matricesε˙αβ and ε˙γγ, there are only six independent components as follows:(19)ε˙xx=∂ux∂x,ε˙yy=∂uy∂y,ε˙zz=∂uz∂z,ε˙xy=ε˙yx=12∂ux∂y+∂uy∂x,ε˙yz=ε˙zy=12∂uz∂y+∂uy∂z,ε˙xz=ε˙zx=12∂ux∂z+∂uz∂x,then ε˙γγ=ε˙xxi+ε˙yyj+ε˙zzk and ε˙αβ=ε˙xyi+ε˙yzj+ε˙xzk:(20)ε˙γγ=∂uα∂xα,ε˙αβ=∂εαβ∂t=12∂uα∂xβ+∂uβ∂xα,where ε˙γγ is the linear strain rate tensor of fluid particle in fluid mechanics in the direction of three axes and ε˙αβ is the fluid particle’s angular deformation rate tensor. Therefore, the tensor ε˙αβ+ε˙γγ is called the strain rate tensor in fluid mechanics, and ω˙αβ+ε˙αβ+ε˙γγ is called the velocity gradient tensor.Moreover, to ensure the weak compressibility and continuity of tailing debris flow, changes in hydrostatic pressurePi within the particle motion domain can be described directly by the equation of state (EOS) in the form of particle density. The SPH equation of state in the form of weak compression was proposed by Monaghan [37], which is as follows:(21)Pi=c02ρ0λρiρ0λ−1,where constant λ=7, ρ0 is the initial average density of tailing debris flow, and c0 is the sound velocity in the corresponding fluid medium. In the current work, c0 = 340 m/s. In the numerical model, it is assumed that the tailing is elastoplastic materials and the mud is weakly compressible viscoplastic fluid. The elastic-plastic constitutive model of the tailings flow is finally established by using the equation of state (21) coupled with (8) and (12) to determine the contact correlation between the mud and the tailings (Figure 2).Figure 2
Modeling of tailing-mud interaction in SPH simulation.


## 2.4. Boundary Condition Used in Current SPH Model

While applying SPH, to save computation time, the boundary was assumed to be pressure-free or constituted of weak-pressure interfaces. However, the idea of “node-to-surface” unidirectional contact control based on compensation by introducing the gradient of momentum and mass into the interpolation kernel function was chosen as it becomes the core idea of dealing with boundary problems [38]. Meanwhile, to avoid particle penetration at the boundary, an improved virtual particle compensation method was introduced. The simulation accuracy was improved substantially by setting virtual particles with no mass and no volume at the boundary to avoid the interference of real boundary particles and the moving fluid particles [39, 40].While applying the virtual particle compensation method, a discrete layer of particles is set at the boundary where the displacement and mass of the boundary particles are assumed to be zero. According to Newton's third law, when the boundary particles come in contact with the fluid particles, repulsive forces are produced in the vertical direction. Consequently, the velocity vector and the displacement vector of the fluid particles change, and the fluid particles are separated and they do not penetrate, as presented in Figure3. This force contribution from the boundary particles is added to the right-hand side of the momentum equation ((6)).Figure 3
Boundary particle distribution in the fictitious particle method.


## 2.5. Operational Framework of SPH Model

The SPH governing equations of the final elastoplastic constitutive model are obtained by substituting the SPH form variable ((8) and ((12)) from Section 2.3 into the momentum (6) and combining with continuity (5). According to the SPH governing equations, the flow chart of the SPH program is shown in Figure 4.Figure 4
Flow chart of SPH program.


## 3. Verification of the Developed SPH Model

### 3.1. Experimental Setup for Tailing Debris Flow Model Test

Since the failure of tailing dams take place catastrophically, the relevant parameters are difficult to measure in the field, and therefore, the laboratory model test is one of the most widespread methods to research tailing debris flow. The tailing debris flow model test has been simulated using the fluid-structure coupled SPH model to verify its numerical accuracy. The schematic of the tailing flow modeling test is presented in Figure5. The experiments were performed in a 4.3 m long × 0.35 m wide × 0.8 m high-discharge flume, in which the horizontal section length was 2.3 m and the slope horizontal projection length was 2.0 m. The stainless steel bottom plate of the flume was 0.35 m wide and two sidewalls were made of 0.25 m high clear acrylic plate to facilitate camera recording. Rigid checking dam with a height of 3.5 cm at 65 cm downstream of the liftgate and concrete structure of length, width, and height of 0.15 m, 0.15 m, and 0.55 m, respectively, were placed at Z = 2.2 m (note that for all the figures in which z-axis is displayed, the origin has been considered corresponding to the liftgate of tailing debris flow reservoir).Figure 5
Schematic diagram illustrating the test device for tailing flow slide and impact.The weight of tailings and mud mixture of 30 kg with a volume concentration ratio of 1 : 3 were placed at the tailing debris flow reservoir with an initial slope towards the outflow boundary on the right-hand side. The average particle size of tailing,d50 = 0.23 m, was obtained from a mine tailing pond in Jiangxi Province, China, and the physical and mechanical properties were estimated in the laboratory. Dynamic load cells and displacement sensors were set up at the dam and downstream structures, and the data were collected simultaneously under multichannels. The camera and slurry height measuring the scale were set up at the relevant places to record the slide and impact behavior of tailing flow.


### 3.2. SPH Simulation Model for Tailing Debris Flow Impact

In this section, the SPH model has been used to simulate the tailing flow impact downstream structure under the action of check dam blocking. The SPH model is consistent with the size of the indoor experimental model and the initial calculation domain of SPH shown in Figure6. The bottom plate of the flume, checking dam, and structure are divided into FEM grids, and SPH particles assignment is used to describe the tailing flow. The debris flow around the tailing pond in the field consists of silt and heavy sand clay, and the behaviors of tailing flow are viscous matrix flow. According to the indoor rheological experiments and real-time debris flow observation in the field, the coefficient of viscosity varies from 0.5 to 2.0 Pa·s [41, 42]. So, the viscosity coefficient of tailing flow in the current work was chosen, μ = 1.0 Pa·s. Internal friction angle between tailing particles ϕ = 26.1 and cohesive strength C = 0.8 kPa. However, the internal friction angle between the mud and the tailing particles was ignored, and only the cohesive strength was considered. Accordingly, the spin parameters, displacement parameters, and other relevant initial parameters were set ((8) and (12)) to obtain the solution of momentum equation (6) and realize the position tracking of the tailing flow particles. To treat the coupling contact between FEM mesh and SPH particle, the former was set as a master segment and the latter was assigned to be the slave segment. Also, a dynamic friction coefficient of 0.12 between each other was assumed, and once the calculation begins, the tailing flow under the action of gravity will start to move from rest.Figure 6
Initial configuration of the numerical computational domain.To carry out the FEM analysis, “solid” rigid material with density, Young’s modulus, and Poisson’s ratio of 7860 kg/m3, 100 GPa, and 0.25, respectively, was selected. For the SPH validation model, the “viscoplastic” fluid materials were assumed to be comprised of about 9000 discretized SPH particles, and the particles were assumed to be uniformly distributed in the whole reservoir. There were 3000 particles, which represent tailing and are made up of “elastic-plastic” material. After defining the initial density and particle arrangement of the material, the continuity equation (5) was solved to realize the mass tracking of tailing flow particles. It was observed from the results of SPH simulation of debris flow [18] that smaller the “MACH number” of controlling fluid compressibility, more accurate are the results of numerical calculation; therefore, MACH number less than 0.3 was recommended. To satisfy the weak compressibility of tailing flow, a pressure reduction coefficient (PC) = 0 was defined in the SPH model. The boundary constraint setting was carried out by applying a layer of virtual particles (about 2000) at X = 0.175 m plane and X = −0.175 m plane, respectively. The distance between the SPH particles was assumed to be 0.01 m. The parameters of the relevant SPH model are calculated and presented in Table 1.Table 1
Parameters used in the SPH simulation of tailing flow.
Mud densityρ (kg/m3)Fine tailing densityρ1 (kg/m3)Coarse tailing densityρ2 (kg/m3)Young’s modulus of tailingsE (GPa)The angle of internal friction of tailingsϕ°Cohesion strengthC (kPa)Total number of particlesNUnit time stepDt (s)Constant applied to the smoothing length of the particlesTerminal timeDT (s)1690169021001626.10.8160004 × 10−41.0514.0


### 3.3. Numerical Comparison with Physical Results

The comparison of the flow patterns of the dam impacted by tailings flow at 0.65 m downstream of the liftgate is presented in Figure7. It was observed that the tested and the simulated results are in good agreement with each other. Further, it was observed from SPH simulation that the tailings with elastic-plastic characteristics carry ascertain the degree of deposition, while the mud particles with viscoplastic characteristics carry some suspended tailings across the blocking device. Array flow characteristics such as “debris flow head” and “debris flow body” of mountain debris flow appeared in the model experiment as well as in the SPH and FEM simulations (Figures 7(b) and 7(c)). As the liftgate was opened, the tailing flow was rapidly flushed forward in the form of a lump of the water-sand mixture, which accumulated strong kinetic energy of impact and formed “debris flow head” (Time = 2.5 s). When time = 2.8 s, the subsequent tailing flow movement became smooth, which is known as “debris flow body.” Besides, the above characteristics also confirmed that it is easy for the saturated tailing reservoir under long-term high water level operation to develop into viscous debris flow after unstable failure [43].Figure 7
Visual comparison of experimental (a) and numerical (SPH (b) and FEM (c)) slurry surface profiles of the dam under impact.
(a)(b)(c)Meanwhile, a comparison between the CFX simulation results based on the FEM grid division and the SPH simulation results is shown in Figure8, and the accumulation and migration processes of the tailing debris flow were consistent. It could be seen that the SPH method was feasible for simulating large deformation of fluid-solid coupling. Besides, in the view of the apparent particle spatter and the discriminating degree of fluid accumulation at the boundary of the solid wall, the SPH method showed better adaptability (Figures 7 and 8). It was noticed that, in the next frame where the time, t = 1.5 s, the “debris flow head” touches the checking dam, which results in a splash of the tailing flow particles and when the time, t = 2.5 s, the peak slurry height reaches 12 cm. The damage caused by the instantaneous impact of the “debris flow head” on the structure is significant, as evident in Figures 8 and 9. Also, the peak slurry height reaches 20 cm at a time, t = 4.5 s. However, increasing the time further to 5.25 s results in a reduction in the height of the slurry present in front of the structure and the particles gather and subsequently flatten.Figure 8
Numerical time history of the density distribution of tailing debris flow based on SPH simulation (upper right corner are the changes of the corresponding flow state under the FEM method).
(a)(b)(c)(d)(e)(f)Figure 9
Visual comparison of experimental and numerical (SPH) slurry surface profiles of the structure under impact.The graph presented in Figure10 depicts a comparison of the numerical and experimental time history of slurry surface elevation acting on checking dam and structure. It can be seen that the SPH and FEM numerical calculation results are consistent with the experimental data, especially at the time range of start to peak appearance. Although the liquid level entered the gentle stage, at time t = 4.5 s and time t = 7 s, the liquid level fluctuation occurred at the checking dam and structure, respectively. This behavior is also evident in the SPH simulation presented in Figures 8(d) and 8(f), respectively. Nevertheless, after the peak value appearance, FEM simulation results were larger than those of SPH and experimental results, and the liquid level fluctuation in FEM simulation was found to be insignificant. The SPH and experimental results have shown better consistency after the peak value appearance. Dai et al. [21] while studying the behavior of debris flow observed that the sudden liquid level fluctuations in the stationary phase are caused by the solid-liquid coupling cohesive action between tailings and mud, which eventually results in large specific gravity tailing transport lag. So, it becomes imperative to complete the material composite setting of viscoplastic mud and elastic-plastic tailings instead of considering the debris flow to be a homogeneous liquid substance.Figure 10
Comparison of numerical and experimental time history of slurry surface elevation acting on the checking dam and structure.By controlling the ratio of particle density to initial density, (21) was utilized to plot the pressure distribution curve. Meanwhile, the experimental pressure was measured by the load cells (electric resistance strain gauge) at the checking dam and the structure. It was noticed from the presented Figure 11 that the numerically calculated pressure was in good agreement with the experimental data, especially the SPH results showed better consistency of the peak values and more obvious characteristic of pressure fluctuation caused by liquid level fluctuation. Also, the instantaneous increase of the impact pressure in a short time is due to the strong destructive energy released at the “debris flow head.” Therefore, the tailing dam and the downstream structure should be separated by a safe distance so that the kinetic energy of the “debris flow head” does not exert a detrimental impact. From this section, it can be concluded that the solid-liquid coupling SPH model can simulate the impact process of tailings flow slide which is similar to mudflow mountain migration characteristics having intermittent continuous flow characteristics and estimate the impact force on checking dam and downstream structure.Figure 11
Comparison of numerical and experimental time history of pressure acting on the checking dam and structure.


## 3.1. Experimental Setup for Tailing Debris Flow Model Test

Since the failure of tailing dams take place catastrophically, the relevant parameters are difficult to measure in the field, and therefore, the laboratory model test is one of the most widespread methods to research tailing debris flow. The tailing debris flow model test has been simulated using the fluid-structure coupled SPH model to verify its numerical accuracy. The schematic of the tailing flow modeling test is presented in Figure5. The experiments were performed in a 4.3 m long × 0.35 m wide × 0.8 m high-discharge flume, in which the horizontal section length was 2.3 m and the slope horizontal projection length was 2.0 m. The stainless steel bottom plate of the flume was 0.35 m wide and two sidewalls were made of 0.25 m high clear acrylic plate to facilitate camera recording. Rigid checking dam with a height of 3.5 cm at 65 cm downstream of the liftgate and concrete structure of length, width, and height of 0.15 m, 0.15 m, and 0.55 m, respectively, were placed at Z = 2.2 m (note that for all the figures in which z-axis is displayed, the origin has been considered corresponding to the liftgate of tailing debris flow reservoir).Figure 5
Schematic diagram illustrating the test device for tailing flow slide and impact.The weight of tailings and mud mixture of 30 kg with a volume concentration ratio of 1 : 3 were placed at the tailing debris flow reservoir with an initial slope towards the outflow boundary on the right-hand side. The average particle size of tailing,d50 = 0.23 m, was obtained from a mine tailing pond in Jiangxi Province, China, and the physical and mechanical properties were estimated in the laboratory. Dynamic load cells and displacement sensors were set up at the dam and downstream structures, and the data were collected simultaneously under multichannels. The camera and slurry height measuring the scale were set up at the relevant places to record the slide and impact behavior of tailing flow.


## 3.2. SPH Simulation Model for Tailing Debris Flow Impact

In this section, the SPH model has been used to simulate the tailing flow impact downstream structure under the action of check dam blocking. The SPH model is consistent with the size of the indoor experimental model and the initial calculation domain of SPH shown in Figure6. The bottom plate of the flume, checking dam, and structure are divided into FEM grids, and SPH particles assignment is used to describe the tailing flow. The debris flow around the tailing pond in the field consists of silt and heavy sand clay, and the behaviors of tailing flow are viscous matrix flow. According to the indoor rheological experiments and real-time debris flow observation in the field, the coefficient of viscosity varies from 0.5 to 2.0 Pa·s [41, 42]. So, the viscosity coefficient of tailing flow in the current work was chosen, μ = 1.0 Pa·s. Internal friction angle between tailing particles ϕ = 26.1 and cohesive strength C = 0.8 kPa. However, the internal friction angle between the mud and the tailing particles was ignored, and only the cohesive strength was considered. Accordingly, the spin parameters, displacement parameters, and other relevant initial parameters were set ((8) and (12)) to obtain the solution of momentum equation (6) and realize the position tracking of the tailing flow particles. To treat the coupling contact between FEM mesh and SPH particle, the former was set as a master segment and the latter was assigned to be the slave segment. Also, a dynamic friction coefficient of 0.12 between each other was assumed, and once the calculation begins, the tailing flow under the action of gravity will start to move from rest.Figure 6
Initial configuration of the numerical computational domain.To carry out the FEM analysis, “solid” rigid material with density, Young’s modulus, and Poisson’s ratio of 7860 kg/m3, 100 GPa, and 0.25, respectively, was selected. For the SPH validation model, the “viscoplastic” fluid materials were assumed to be comprised of about 9000 discretized SPH particles, and the particles were assumed to be uniformly distributed in the whole reservoir. There were 3000 particles, which represent tailing and are made up of “elastic-plastic” material. After defining the initial density and particle arrangement of the material, the continuity equation (5) was solved to realize the mass tracking of tailing flow particles. It was observed from the results of SPH simulation of debris flow [18] that smaller the “MACH number” of controlling fluid compressibility, more accurate are the results of numerical calculation; therefore, MACH number less than 0.3 was recommended. To satisfy the weak compressibility of tailing flow, a pressure reduction coefficient (PC) = 0 was defined in the SPH model. The boundary constraint setting was carried out by applying a layer of virtual particles (about 2000) at X = 0.175 m plane and X = −0.175 m plane, respectively. The distance between the SPH particles was assumed to be 0.01 m. The parameters of the relevant SPH model are calculated and presented in Table 1.Table 1
Parameters used in the SPH simulation of tailing flow.
Mud densityρ (kg/m3)Fine tailing densityρ1 (kg/m3)Coarse tailing densityρ2 (kg/m3)Young’s modulus of tailingsE (GPa)The angle of internal friction of tailingsϕ°Cohesion strengthC (kPa)Total number of particlesNUnit time stepDt (s)Constant applied to the smoothing length of the particlesTerminal timeDT (s)1690169021001626.10.8160004 × 10−41.0514.0


## 3.3. Numerical Comparison with Physical Results

The comparison of the flow patterns of the dam impacted by tailings flow at 0.65 m downstream of the liftgate is presented in Figure7. It was observed that the tested and the simulated results are in good agreement with each other. Further, it was observed from SPH simulation that the tailings with elastic-plastic characteristics carry ascertain the degree of deposition, while the mud particles with viscoplastic characteristics carry some suspended tailings across the blocking device. Array flow characteristics such as “debris flow head” and “debris flow body” of mountain debris flow appeared in the model experiment as well as in the SPH and FEM simulations (Figures 7(b) and 7(c)). As the liftgate was opened, the tailing flow was rapidly flushed forward in the form of a lump of the water-sand mixture, which accumulated strong kinetic energy of impact and formed “debris flow head” (Time = 2.5 s). When time = 2.8 s, the subsequent tailing flow movement became smooth, which is known as “debris flow body.” Besides, the above characteristics also confirmed that it is easy for the saturated tailing reservoir under long-term high water level operation to develop into viscous debris flow after unstable failure [43].Figure 7
Visual comparison of experimental (a) and numerical (SPH (b) and FEM (c)) slurry surface profiles of the dam under impact.
(a)(b)(c)Meanwhile, a comparison between the CFX simulation results based on the FEM grid division and the SPH simulation results is shown in Figure8, and the accumulation and migration processes of the tailing debris flow were consistent. It could be seen that the SPH method was feasible for simulating large deformation of fluid-solid coupling. Besides, in the view of the apparent particle spatter and the discriminating degree of fluid accumulation at the boundary of the solid wall, the SPH method showed better adaptability (Figures 7 and 8). It was noticed that, in the next frame where the time, t = 1.5 s, the “debris flow head” touches the checking dam, which results in a splash of the tailing flow particles and when the time, t = 2.5 s, the peak slurry height reaches 12 cm. The damage caused by the instantaneous impact of the “debris flow head” on the structure is significant, as evident in Figures 8 and 9. Also, the peak slurry height reaches 20 cm at a time, t = 4.5 s. However, increasing the time further to 5.25 s results in a reduction in the height of the slurry present in front of the structure and the particles gather and subsequently flatten.Figure 8
Numerical time history of the density distribution of tailing debris flow based on SPH simulation (upper right corner are the changes of the corresponding flow state under the FEM method).
(a)(b)(c)(d)(e)(f)Figure 9
Visual comparison of experimental and numerical (SPH) slurry surface profiles of the structure under impact.The graph presented in Figure10 depicts a comparison of the numerical and experimental time history of slurry surface elevation acting on checking dam and structure. It can be seen that the SPH and FEM numerical calculation results are consistent with the experimental data, especially at the time range of start to peak appearance. Although the liquid level entered the gentle stage, at time t = 4.5 s and time t = 7 s, the liquid level fluctuation occurred at the checking dam and structure, respectively. This behavior is also evident in the SPH simulation presented in Figures 8(d) and 8(f), respectively. Nevertheless, after the peak value appearance, FEM simulation results were larger than those of SPH and experimental results, and the liquid level fluctuation in FEM simulation was found to be insignificant. The SPH and experimental results have shown better consistency after the peak value appearance. Dai et al. [21] while studying the behavior of debris flow observed that the sudden liquid level fluctuations in the stationary phase are caused by the solid-liquid coupling cohesive action between tailings and mud, which eventually results in large specific gravity tailing transport lag. So, it becomes imperative to complete the material composite setting of viscoplastic mud and elastic-plastic tailings instead of considering the debris flow to be a homogeneous liquid substance.Figure 10
Comparison of numerical and experimental time history of slurry surface elevation acting on the checking dam and structure.By controlling the ratio of particle density to initial density, (21) was utilized to plot the pressure distribution curve. Meanwhile, the experimental pressure was measured by the load cells (electric resistance strain gauge) at the checking dam and the structure. It was noticed from the presented Figure 11 that the numerically calculated pressure was in good agreement with the experimental data, especially the SPH results showed better consistency of the peak values and more obvious characteristic of pressure fluctuation caused by liquid level fluctuation. Also, the instantaneous increase of the impact pressure in a short time is due to the strong destructive energy released at the “debris flow head.” Therefore, the tailing dam and the downstream structure should be separated by a safe distance so that the kinetic energy of the “debris flow head” does not exert a detrimental impact. From this section, it can be concluded that the solid-liquid coupling SPH model can simulate the impact process of tailings flow slide which is similar to mudflow mountain migration characteristics having intermittent continuous flow characteristics and estimate the impact force on checking dam and downstream structure.Figure 11
Comparison of numerical and experimental time history of pressure acting on the checking dam and structure.


## 4. Numerical Simulation of Checking Dam Height Effect

When the tailing flow slide across the checking dam is stopped, then the intercepted sediment approximates a triangular wedge as presented in Figure8(f). The tailing flow interception volume under the action of checking dam which is 3.5 cm in height accounts for about 30% of the tailing debris flow reservoir volume. Therefore, to reduce the impact of tailing debris flow on downstream structures, it is necessary to study the effect of checking dam height on debris flow disaster. In this context, the verified SPH solid-liquid coupling simulation model can be used as presented in Figure 6 to modify the height of checking dam h from 0 cm, 2 cm, 3.5 cm, and 5 cm to 6.5 cm at 65 cm downstream of the liftgate. The calculation parameters are consistent with Table 1.


### 4.1. Effect of Checking Dam Height on Debris Flow

The simulation results of the impact pattern of tailing flow under different checking dams heights are presented in Figure12. At the time, t = 1.75 s, four groups of numerical simulation barring simulation 1 (h = 0) produce particles spatter of varying degrees under the action of checking the dam. With the increase in checking dam height, the impact of “debris flow head” is enhanced, and it also takes more time for tailing flow to pass through the dam. At the time, t = 4.9 s, under the action of checking dam and tailing-mud cohesion, some of the tailings were accumulated near the dam.Figure 12
Comparison of impact form of tailing flow under different checking dam heights.


### 4.2. Run-Up and Impact Force

The influence of checking dam height exerted on both slurry run-up and the impact force acting on the downstream structure have been presented in Figures13 and 14. It was noted that the density field variable of tailing flow was obtained by solving the continuity equation (5). To assess the pressure distribution, the equation of state (21) was solved. Meanwhile, the area integral of pressure was performed to estimate the magnitude of impact force acting on the structure. Due to the buffer effect, while checking the dam, a slightly different spatial distribution of these particles could perhaps result in a different pressure determination at the front end of the downstream structure. The slurry run-up peak on the structure, shown in Figure 13, is greater for h = 0 cm than it is for the other values of h. Consequently, the associated impact force is also greater for the condition of checking dam existent. With the increase in the checking dam height, the rate of the slurry run-up in front of the downstream structure decreased, and a lag in the peak value was observed. Nevertheless, further investigation is required to properly elucidate this variability.Figure 13
Numerical time history of slurry-surface elevation acting on the structure under different checking dam conditions.Figure 14
Numerical time history of impact force acting on the structure under different checking dam conditions.So, the influence of checking dam height on both slurry run-up peak and the impact force peak acting on the checking dam and downstream structure have been investigated and presented in Figures15 and 16. The blocking height has a significant buffering effect on the breaking tailing flow, and it varies linearly with dam height. Increasing the dam height by 1.5 cm, the peak value of slurry run-up at the downstream structure is reduced by 15.95%, and the associated average reduction of the peak value of the impact force is 30.75%. On the other hand, the effective blocking efficiency at checking the dam also shows a linear distribution and the absolute slope value is consistent with the associated linear distribution at the downstream structure.Figure 15
Slurry depth peak vs checking dam’s height acting on the structure and checking dam.Figure 16
Impact force peak vs checking dam’s height acting on the structure and checking the dam.


### 4.3. Velocity Distribution around the Structure

As the impact forces, presented in Figure14 at the time, t = 4.55 s, for different checking dam heights correspond to key moments in the tailing flow-structure interaction, the velocity vectors in the x-direction are presented sequentially in Figure 17 to investigate the velocity distribution around the structure. When tailing debris flow impinges on the structure, the forward flow particles rapidly extrude and get divided into either side and the velocity of flow in the x-direction increases rapidly, which results in forming a nearly symmetrical “velocity-dense zone” (shown in yellow and blue color) at ±45° corners in front of the structure. However, at the rear end of the structure, there is a slight presence of the “velocity-dense zone” in the opposite direction, and an intermediate confluence is formed as the wake downstream from the structure develops. This occurs because the reverse flow rate component is produced by the reaction force of the boundary particles according to Newton’s third law [44]. The divided tailing flow is then reoriented towards the centerline of the flume, and since the X component of velocity is small as compared with the Z component, the divided tailing flow eventually combines further while moving downstream.Figure 17
X-velocity field in a plane close to structure after gate opening for time, t = 4.55 s, for the following check dam heights: (a) 2 cm; (b) 3.5 cm; (c) 5 cm; (d) 6.5 cm.
(a)(b)(c)(d)Figure18 shows the comparison of numerically calculated pressure cross-sectional view of tailing flow-structure interaction 0.01 m upstream of structure for time t = 4.55 s under the action of different checking dam heights. As the checking dam height increases, at the time, t = 4.55 s, the range of “velocity-dense zone” decreases substantially. The divided tailing flow flowing around the structure reconnect later in the simulations is lagged (Figures 17(b)–17(d)) and less impact via kinetic energy is exerted on the side pillars of the structure for dam height h = 2 cm (Figures 17(a) and 18(a)). Although the pressure distribution before the structure is weakened by the checking dam, the concentrated distribution of pressure at side pillars is obvious (Figure 18)). From the point view of disaster prevention and control, the appearance of “velocity-dense zone” and cross-sectional pressure provides a certain reference basis for increasing the stiffness of both sides of the pillar in the structural design. Combining with the blocking dam set up upstream, the system can slow down the erosion of the side pillar foundation in the front of the structure by tailing flow and therefore avoid cracking and failure of the walls near the pillars.Figure 18
Numerical cross-sectional view of tailings flow-structure interaction 0.01 m upstream of the structure after gate opening for time,t = 4.55 s, for the following check dam heights. (a) 2 cm; (b) 3.5 cm; (c) 5 cm; (d) 6.5 cm.
(a)(b)(c)(d)


## 4.1. Effect of Checking Dam Height on Debris Flow

The simulation results of the impact pattern of tailing flow under different checking dams heights are presented in Figure12. At the time, t = 1.75 s, four groups of numerical simulation barring simulation 1 (h = 0) produce particles spatter of varying degrees under the action of checking the dam. With the increase in checking dam height, the impact of “debris flow head” is enhanced, and it also takes more time for tailing flow to pass through the dam. At the time, t = 4.9 s, under the action of checking dam and tailing-mud cohesion, some of the tailings were accumulated near the dam.Figure 12
Comparison of impact form of tailing flow under different checking dam heights.


## 4.2. Run-Up and Impact Force

The influence of checking dam height exerted on both slurry run-up and the impact force acting on the downstream structure have been presented in Figures13 and 14. It was noted that the density field variable of tailing flow was obtained by solving the continuity equation (5). To assess the pressure distribution, the equation of state (21) was solved. Meanwhile, the area integral of pressure was performed to estimate the magnitude of impact force acting on the structure. Due to the buffer effect, while checking the dam, a slightly different spatial distribution of these particles could perhaps result in a different pressure determination at the front end of the downstream structure. The slurry run-up peak on the structure, shown in Figure 13, is greater for h = 0 cm than it is for the other values of h. Consequently, the associated impact force is also greater for the condition of checking dam existent. With the increase in the checking dam height, the rate of the slurry run-up in front of the downstream structure decreased, and a lag in the peak value was observed. Nevertheless, further investigation is required to properly elucidate this variability.Figure 13
Numerical time history of slurry-surface elevation acting on the structure under different checking dam conditions.Figure 14
Numerical time history of impact force acting on the structure under different checking dam conditions.So, the influence of checking dam height on both slurry run-up peak and the impact force peak acting on the checking dam and downstream structure have been investigated and presented in Figures15 and 16. The blocking height has a significant buffering effect on the breaking tailing flow, and it varies linearly with dam height. Increasing the dam height by 1.5 cm, the peak value of slurry run-up at the downstream structure is reduced by 15.95%, and the associated average reduction of the peak value of the impact force is 30.75%. On the other hand, the effective blocking efficiency at checking the dam also shows a linear distribution and the absolute slope value is consistent with the associated linear distribution at the downstream structure.Figure 15
Slurry depth peak vs checking dam’s height acting on the structure and checking dam.Figure 16
Impact force peak vs checking dam’s height acting on the structure and checking the dam.


## 4.3. Velocity Distribution around the Structure

As the impact forces, presented in Figure14 at the time, t = 4.55 s, for different checking dam heights correspond to key moments in the tailing flow-structure interaction, the velocity vectors in the x-direction are presented sequentially in Figure 17 to investigate the velocity distribution around the structure. When tailing debris flow impinges on the structure, the forward flow particles rapidly extrude and get divided into either side and the velocity of flow in the x-direction increases rapidly, which results in forming a nearly symmetrical “velocity-dense zone” (shown in yellow and blue color) at ±45° corners in front of the structure. However, at the rear end of the structure, there is a slight presence of the “velocity-dense zone” in the opposite direction, and an intermediate confluence is formed as the wake downstream from the structure develops. This occurs because the reverse flow rate component is produced by the reaction force of the boundary particles according to Newton’s third law [44]. The divided tailing flow is then reoriented towards the centerline of the flume, and since the X component of velocity is small as compared with the Z component, the divided tailing flow eventually combines further while moving downstream.Figure 17
X-velocity field in a plane close to structure after gate opening for time, t = 4.55 s, for the following check dam heights: (a) 2 cm; (b) 3.5 cm; (c) 5 cm; (d) 6.5 cm.
(a)(b)(c)(d)Figure18 shows the comparison of numerically calculated pressure cross-sectional view of tailing flow-structure interaction 0.01 m upstream of structure for time t = 4.55 s under the action of different checking dam heights. As the checking dam height increases, at the time, t = 4.55 s, the range of “velocity-dense zone” decreases substantially. The divided tailing flow flowing around the structure reconnect later in the simulations is lagged (Figures 17(b)–17(d)) and less impact via kinetic energy is exerted on the side pillars of the structure for dam height h = 2 cm (Figures 17(a) and 18(a)). Although the pressure distribution before the structure is weakened by the checking dam, the concentrated distribution of pressure at side pillars is obvious (Figure 18)). From the point view of disaster prevention and control, the appearance of “velocity-dense zone” and cross-sectional pressure provides a certain reference basis for increasing the stiffness of both sides of the pillar in the structural design. Combining with the blocking dam set up upstream, the system can slow down the erosion of the side pillar foundation in the front of the structure by tailing flow and therefore avoid cracking and failure of the walls near the pillars.Figure 18
Numerical cross-sectional view of tailings flow-structure interaction 0.01 m upstream of the structure after gate opening for time,t = 4.55 s, for the following check dam heights. (a) 2 cm; (b) 3.5 cm; (c) 5 cm; (d) 6.5 cm.
(a)(b)(c)(d)


## 5. Conclusions

In the present research work, the Bingham flow equation of motion and mud-tailing constitutive model was modified, and smooth particle hydrodynamics (SPH) discrete formula of N-S equation is established to evaluate the impact simulation of tailing flow on the downstream structure under the action of checking dam. The following conclusions can be drawn from the aforementioned study:(1)
The impact process of the tailing flow showed the characteristics of “debris flows head” and “debris flows body” which are similar to the migration characteristics of mountain debris flow. The mud properties showed the formation of continuous flow, which has typical Bingham flow characteristics. The “debris flows head” exerted a substantial instantaneous impact on the structure; consequently, impact damage was easily caused to the downstream structure. To avoid the aforementioned problem, the tailing dam and the downstream structure should maintain a safe distance and meet the requirements of safety regulations.(2)
The numerical calculation results, based on the SPH method, of impact process of the tailing flow, were in good agreement with the experimental results, especially the depth of mud inundation and the impact strength before the structure. The simulation program by using an SPH method based on the elastic-plastic constitutive model had good applicability in solving geotechnical large deformation problems of similar tailing flow slides. The SPH algorithm can be used to dynamically reflect the disaster-forming process, predict and evaluate the impact damage of the affected bodies under different disaster intensities, and lay a foundation for early warning and scientific guidance for disaster prevention and mitigation.(3)
The SPH simulation model was used to analyze the buffer effect of checking dam height on the dam break tailing flow. The results showed that the blocking height had a linear distribution on the buffer effect of the tailing flow, and changes in the “velocity-dense zone” in the range of ±45° on the front part of the structure directly affect the disaster-bearing capacity of the side pillar foundation of the structure.However, the SPH numerical model can further be developed better to employ the pore water pressure in tailing flow-like regime since the present SPH model is still controlled by hydrostatic pressure which ignores the influence of pore water pressure in tailing-mud interaction. Besides, the detailed analysis of the viscosity coefficient and its effects on the tailings flow process is to be conducted in the further works.



---

*Source: 9304921-2020-07-07.xml*