Numerical simulation of impact and entrainment behaviors of debris flow by using SPH–DEM–FEM coupling method

2.1 Governing equations

For the special solid–liquid flow, such as debris flows, the slurry velocities before and after the debris particles may differ because of the large difference between the sizes of the debris particles. The difference in the slurry velocity, in turn, forms a shape resistance and causes a certain difference in the contact resistance within the debris particles. In contrast, the set of rheological equations for solid–liquid flows, which is usually a combination of the continuous equations based on computational fluid dynamics (CFD) and Navier–Stokes (N–S) equations, does not consider the effect of gravel particle size differentiation and cannot adequately determine the mass and momentum conservation characterization of the debris flow. Therefore, based on the original CFD rheological control equation, this article introduces contact resistance such as drag force and buoyancy force to achieve momentum exchange and dynamic equilibrium between continuous slurry and discontinuous debris particles [50,51,52].

The continuity and momentum equations for fluids are as follows [50,51]:

where ρ, v, P, τ, and g represent density, velocity vector, pressure, shear stress tensors, and gravitational acceleration vector, respectively. f f d and f f b represent drag force and buoyancy per unit volume in the liquid phase. ∇ is a divergence symbol, and the subscript f represents “fluid.” D D t = ∂ ∂ t + v ⋅ ∇ is the substantial derivative.

The motion of solid particles is governed by the following equation [51,52]:

where f s d , f s c , and f s m represent drag force per unit volume in the solid phase, net force per unit volume of the interaction forces between the solid particles, and other miscellaneous effects (e.g., lubrication force or buoyancy on solid particle), respectively. The subscript s refers to “solid particle.”

As shown in Figure 1, both slurry and debris particles carried in the debris flow exist in the form of particle flow; the discrete distribution, instead of the continuum problem domain, is apparently more appropriate for analyzing the impact of debris flow on the structure. Therefore, we further improve equations (1)–(3), and we discretely obtain the rheological control equations in the form of SPH and DEM. The treatment of contact relations for solid–liquid two-phase flows is introduced in Section 3.

Figure 1

Modeling of debris–slurry–structure interaction in SPH–DEM–FEM simulation.

2.2 Slurry governed by SPH

In the coupled SPH–DEM method, fluids are represented by SPH particles, solid particles entrained by fluid are described by DEM particles, and the free distribution of discrete particles is used instead of its problem domain. Each particle is assigned its own variable information (abscissa and ordinate, material density, velocity, mass, stress tensor, etc.). A discrete modification of the continuity equation and the N–S equation are introduced to realize particle flow control. Based on the spline interpolation kernel function, the random interaction between particles is simplified to obtain the SPH equation [53,54]:

where x i α and v i j α = v i α − v j α represent the position and velocity vector, respectively; g ^α is gravitational acceleration; σ i α β is the stress tensor of SPH particles; ρ is the density of fluid; Fd _s→f,i , Fb _s→f,i represent the drag force and buoyancy exerted on SPH particle i by DEM particles, and details of SPH–DEM interaction are introduced in Section 3; α and β represent the Cartesian components x, y, z; subscripts i and j represent arbitrary SPH (fluid) particles.

The kernel function W _ij selects the different dimensional forms of the kernel function of cubic spline interpolation with the change in α _d [15]:

where r _ij = x _i − x _j represents the distance vector; | r _ij | = |x _i − x _j | = L ₀; R = |r_ij|/h at position x i α (x, y, z); α _d represents a constant defined for one-, two-, and three-dimensional spaces as in α d = 1 h , 15 7 π h 2 , and 3 2 π h 3 , respectively.

The smooth length h and influence radius of the single particle are smaller, reducing the number of interpolation calculation iteration while saving computer storage. However, the smooth length h must be addressed for particle spacing L ₀(L ₀ < h ≤ 1.5L ₀) to ensure that related particles exist within the influence range and to prevent particles from penetrating into the finite element mesh in the coupling calculation [30]. Considering these requirements and cost of calculation, an economical and reasonable value of h = 1.05L ₀ was selected in this study.

In investigating large deformations of the solid–liquid flow based on the SPH–DEM model (i.e., debris flow), the solid and liquid phases are often connected by fluid–solid interaction forces that contain the drag and buoyancy of the particles. However, the internal flow resistance is expressed by viscous shear stress and needs to be implemented through a material model [55,56]. In this study, we selected two material models to define the internal flow resistance.

The null material (*MAT_NULL in LSDYNA) is generally defined for water, whose stress–strain relationship can be represented by the Bingham flow constitutive control equation [57]:

where μ _eff is effective viscosity and ε ^αβ is shear strain tensor.

where μ is the apparent viscosity coefficient; τ _y represents the yield strength, which controls the flow state before fluid instability failure and satisfies the Mohr–Coulomb yield criterion; ε _∏ denotes the second invariant of the shear strain tensor; ε ̇ is the rate of shear strain; m and n are the constant and power law index controlling the growth of stress at low and high shear rates, respectively. However, the null material model ignores yield strength (τ _y = 0) and behaves like a fluid. Therefore, the null material model can be considered as a special case of the Bingham model.

As slurry is typically considered as a homogeneous continuous mixture of water and clay and exhibits elastic–plastic behavior patterns [58,59,60], slurry can be characterized by the elastic–plastic–hydrodynamic material (*MAT_ELASTIC_PLASTIC_HYDRO [MEPH]),which exhibits continuous behaviors of fluid and elastic–plastic solid. The stress–strain relationship of slurry is based on the elastic–plastic constitutive model, as follows:

where P represents hydrodynamic pressure, δ ^αβ is Dirac delta function, and s ^αβ is partial stress tensor. The rate of change of partial stress tensor is based on the Jaumann criterion [58]:

where G, w ̇ , and ε ̇ represent shear modulus, spin rate tensor, and shear strain rate tensor, respectively. To simplify the numerical calculation, the internal fluid resistance of slurry in the MEPH material model can be characterized by the yield stress; equation (10) can be modified for a time step Δt as follows [29,61]:

where R ^αβ is the rotation matrix containing basis vectors. The effective value of the deviatoric shear stress s _eff satisfies the expression of the deviatoric shear stress tensor ( s ^αβ ) _t+Δt as follows:

In the MEPH material model, the s _eff and the yield stress τ _y are compared to determine whether the material undergoes plastic deformation, which affects the output of the deviatoric shear stress tensor.

where Δε _p and E _h represent plastic strain increment and plastic hardening modulus, respectively. Finally, the stress tensor σ t + Δ t α β is calculated as follows:

In addition, Peng et al. [62] applied the equation of state (*EOS_MURNAGHAN) to implement the SPH method for the simulation of large-scale landslide accidents and verified the EOS of fluid in SPH is feasible for describing large deformation problems of soil (e.g., landslides and debris flows). To calculate the pressure, the following equation of state (*EOS_MURNAGHAN) is used to model weakly compressible fluid flow with SPH elements [61,62]. Small changes in relative density yield large pressure; thus, simulations reproduce the dynamic pressure changes in motion.

where the constant λ = 7; ρ ₀, ρ _i , and c ₀ = c(ρ ₀) represent the initial density of the fluid, density of the fluid in the calculation domain of SPH particles, and sound velocity related to initial density, respectively.

Based on the coupled equation of state (equation (16)), stress–strain equation (equation (7) or equation (9)) of the liquid phase, and related equation improvement, the solution of the SPH governing equation (equations (4) and (5)) can be obtained to establish the rheological model of the slurry.

2.3 Debris particle governed by DEM

DEM is a meshless method that was derived from the development of molecular dynamics theory proposed by Cundall and Strack [17] and first used in the field of rock mechanics. The basic principle is to determine the force and displacement of particles at each time step based on the interaction force and Newton’s law of motion. The aim is to obtain the macroscopic motion law of particle flow by tracking the microscopic motion of single particles.

Based on the governing equation of solid particles in the solid–liquid flows, the dominant forces exerted on the solid phase are fluid–solid interaction forces and those between the solid particles [50]. In this study, we used spherical DEM particles to simulate debris particles and weighted various physico-mechanical metrics of debris particles according to the diameter size of the DEM particles. As actual debris particles are not strictly spheres, we adopted contact torque to replace the effect of irregular debris particles, and we approximately simulated particle rolling by updating angular changes. Thus, the translation and rotation governing equations of arbitrary solid particle are written as follows:

where Fc _kl represents the interaction forces between solid particles and m _k refers to mass. Fd _f→s,k and Fb _f→s,k are drag force and buoyancy, respectively, which are exerted on solid particle k by fluid. The quantities are defined in Section 3. I k = m k d k 2 / 10 is rotary inertia, d _k is particle diameter, ω _k denotes rotational velocity, and Tc _kl is torque of interaction forces. The subscript k and l denotes arbitrary DEM (solid) particles.

To describe the displacement variation of solid particles, a first-order Taylor expansion of the velocity is used, and the specific derivation procedure is shown as follows:

The equation of the combination is then obtained after substituting equation (20) into equation (19):

where Δt is the time step and Δt → 0. Subsequently, the effect of the quadratic part of Δt is negligible, while the velocity at time t + Δ t 2 is resolved as follows:

Similarly, the iterative expression of the angular velocity is obtained as follows:

Finally, the position and angle of the DEM particle are updated by integrating time t accordingly:

In DEM, solid particles typically make contact through a joint operation between two series damping springs, as shown in Figure 2 [63,64]. Thus, by adding the forces generated by the spring damper in the normal and tangent directions, the interparticle contact force Fc _kl in DEM can be obtained.

where F c kl n represents the normal contact forces and F c kl s represents tangential contact forces. The interparticle damping reflects the energy dissipation during the collision, and considering the particle force under damping, it can be expressed as follows [63]:

where K, Δx, η, and Δv are the coefficient of elasticity, displacement increment, damping coefficient of the dashpot, and velocity increment, respectively. The superscripts n and s represent the normal and tangential directions, respectively. In this study, we used the solid material model (i.e., *MAT_FHWA_SOIL or *MAT_ELASTIC) to describe the stress of debris particle, and the damping and elasticity coefficients were defined by element control (*CONTROL_ DISCRETE_ELEMENT).

Figure 2

Illustration of DEM contact model for solid–solid particles.

2.4 Structures governed by FEM

Figure 3 illustrates the structural damage that occurs when the debris flow contacts the check dam. In the SPH–DEM–FEM coupling model, the structure is modeled by the FEM elements, while the failure process of the structure is simulated based on the element erosion algorithm. We adopted the material model (*MAT_PLASTIC_KINEMATIC) to effectively describe the isotropic and kinematic hardening of the plastic strain model; the material model is suitable for the impact simulation of beams, shells, and solid elements. Several failure criteria are available in the platform, including effective plastic strain, critical principal stress, and critical principal strain [45]. In this study, we selected the effective plastic strain to determine structural failure:

where ε eff p is the effective plastic strain; ε f p represents effective plastic strain at failure, which can be determined in the initial material definition of the model system; and ε ̇ p is the plastic part of the strain rate tensor.

Figure 3

Illustration of structural damage caused by coupling impact between debris flow and check dam.

Satisfying equation (29) demonstrates that the structure has failed. Meanwhile, the failed FEM element no longer bears any load and is automatically deleted in the later node calculation, as illustrated in Figure 3. This failure criterion is a simplification of the destruction process, and the failure element is deleted instead of being cracked or separated; however, the failure criterion can also ensure the continuous operation of the numerical calculation while sufficiently reflecting the contact damage of the debris flow and the structure. To minimize the calculation cost, the rigid material model (*MAT_RIGID) is selected to simulate some weakly deformed or unimportant buildings. Although different material models have different requirements for deformation control, the impact forces need to be calculated, and the corresponding calculation parameters should be inputted into the simulation model.

3.1 Coupled SPH–DEM

During debris flow motion, the coupling contact of the solid–liquid phase is highly complex and involves impact destruction of debris flow, structure buffering, and postdiffusion accumulation distribution. To investigate the slippery shift of debris flows, we developed a numerical model of solid–liquid coupling using SPH–DEM and adopted the drag force and buoyancy as contact forces (as shown in Figure 4), which connected the slurry particle (SPH) and debris particle (DEM). The node-to-node contact (*DEFINE_SPH_DE_COUPLING in LS-DYNA) between the solid and fluid was selected to couple the SPH and the DEM solvers.

Figure 4

Contact force and torque of the slurry particles acting on the debris particles.

The drag force is related to the relative flow velocity among the solid and liquid phases, the local density of neighboring DEM particles in the calculation domain of the SPH particle, and the effective cross-section area of the debris particles. Based on the current theory, the debris particles in mudslides are approximately spherical. Referring to the two-fluid calculation model [35,52,65], the drag force impact on the DEM particle by the SPH particle is expressed as follows:

where v ¯ SPH,k α denotes the average velocity of SPH particles near the DEM particle k, v DEM,k α represents the velocity of DEM particle k, ε _k is the porosity of the DEM particle determined by the specific gravity of the DEM particles within the unit volume, and A _k represents the effective cross-section area of the DEM particle. The fluid resistance coefficient ζ _k is given by equation (32). The upper equation [66] is used in regions with ε _k ≤ 0.8, and the lower equation [67] is applied in regions with ε _k > 0.8. ν _f, d _k, C _d, and ρ _f denote the fluid dynamic viscosity, diameter of the DEM particle, drag coefficient of the DEM particle, and the reference density of slurry in debris particles for determining saturation, respectively.

An approximate expression of C _d was obtained by linearizing the N–S equation, but the equation is only applicable in cases of low Reynolds numbers (Re < 0.4). For high Reynolds numbers occurring in the debris flow, the relationship between C _d and Re can be determined only through multiple tests. Therefore, the C _d values in this study were characterized using Brown and Lawler empirical equations [68] that are suitable in turbulent, laminar, and transition regions.

When Re < 3.5, the friction velocity ( ∣ v ¯ SPH,k α − v DEM,k α ∣ ) is considerably small, the square term is negligible, and the value of the drag force is dominated by the friction of the particle surface. When Re ≥ 3.5, the friction velocity gradually increases, and the top streamline of the particle is separated, resulting in a pressure difference between the front and back and generating a shape resistance [69]. The drag force is generally the combined force of the surface friction F ₁ and the shape resistance F ₂, as shown in Figure 4. Equation (31) is compatible with the aforementioned two scenarios.

Buoyancy, another component of the contact force in the SPH–DEM coupling model, satisfies the following expression in hydrostatic equilibrium:

Hydrostatic pressure gradient ∇P _SPH,i does not apply to the equation of state (equation (16)). By introducing the viscosity of the hydrodynamics and the relative flow velocity, the upper equation suitable for SPH–DEM coupling model can be written as follows [64]:

Since Newton’s Third Law governs the interaction force, the drag force Fd _s→f,I and buoyancy Fb _s→f,I exerted by the DEM particle on the SPH particle and the reacting force on the SPH fluid by the DEM object are of equal values but opposite directions.

3.2 Coupled SPH–FEM

To implement the contact analysis between slurry and structure, a penalty method based on a master-slave algorithm was adopted. In this coupling contact model, SPH particles (slurry) represent the slave segment, FEM elements (structure) represent the master segment, and the relevant contact control parameters (i.e., friction coefficient, scale factor on default slave penalty stiffness, and scale factor on default master penalty stiffness) are applied to evaluate the interactions. By automatically searching for the contact displacement and relative velocity on contact surfaces combined with contact stiffness and damping coefficient (as shown in Figure 5(a)), the expression of contact force is established [29,45,70]:

where Fc _SPH→FEM; K, η, Δ x , and Δ v represent the contact force between the SPH particle and FEM element, contact stiffness, damping coefficient, contact displacement vector, and relative velocity, respectively. For solid elements in FEM, the normal contact force F c sf n and tangential contact force F c sf s can be written as follows:

Figure 5

Illustration of interaction schemes of contact elements: (a) SPH particle and FEM element contact and (b) DEM particle and FEM element contact.

For shell elements in FEM,

where α = 0.1 is a default scale factor on contact stiffness; μ, G, V, and A represent sliding friction coefficient, volume modulus of FEM element, volume of FEM element, and contact area of master segment surface, respectively. To avoid the stiffness difference between the master and slave segments interrupting the calculation, we adopted the soft constraint penalty formulation to expand the definition of contact stiffness. The program automatically searches for the upper limit of the contact stiffness in calculating the contact force. SOFSCL is the scale factor for the soft constraint penalty formulation and is defaulted to 0.1. m _i is the mass of the SPH particle on the FEM surface, and Δt is the time step.

3.3 Coupled DEM–FEM

In the coupled SPH–FEM model, some SPH particles penetrate the contact surface, and the related contact force is mainly determined by the displacement generated when the penetration occurs. For the contact between the DEM particles and FEM elements, the DEM particle size is large and cannot fully penetrate the structure; however, a certain contact overlap exists between the DEM particle k and contact surface (as shown in Figure 5(b)). Hence, the contact between the DEM particles (debris) and FEM elements (structure) can also be expressed by the master–slave algorithm. The relevant contact force is calculated as in the SPH–FEM contact described in Section 3.2 [45].

where Fc _DEM→FEM represent the contact force between the DEM particle and FEM element. Δ x _df = Δ x _k denotes the overlap distance vector between the initial contact point x _k,t and the change contact point x _k,t+Δt after time step Δt. When positive, this overlap is used to calculate a contact force Fc _DEM→FEM using the following equations. For solid elements in FEM, the normal contact force F c sf n and tangential contact force F c sf s can be expressed as follows:

For shell elements in FEM,

where m _k, η _df, and Δ v _df represent the mass of the DEM particle, contact damping coefficient of DEM–FEM, and relative velocity of the DEM particle on the FEM surface, respectively. The superscripts n and s represent normal and tangential directions, respectively.

3.4 Operating process of the coupling model

SPH–DEM–FEM simulation is a method based on the Euler–Lagrangian reference frame, which is a new attempt at discrete simulation. The flow chart of the SPH program is shown in Figure 6 based on the crucial issues considered earlier in the SPH–DEM–FEM numerical simulation. The computation procedure is summarized in the following steps:

Figure 6

Flow chart of SPH–DEM–FEM program.

Step 1. The initial stress values of the SPH particles at time step t + Δt are obtained from the stress–strain relationship equations of the elastic–plasticity intrinsic model established in Section 2.2 (see equations (14)–(16)). The equations of motion control (equations (4) and (5)) are solved to generate the motion potential of the SPH particles.

Step 2. The contact reaction forces between SPH–DEM and SPH–FEM are obtained by calculating the coupled contact equations (31), (37), (38), and (43), which produce erosive entrainment of DEM particles and cause a deformation of FEM elements.

Step 3. Time-step iteration is used to continuously update the motion trajectories of discrete particles (SPH and DEM particles) and the deformation failure of FEM elements.

Step 4. The moving discrete particles are filled or bypassed according to the boundary shape constructed by the FEM elements; the discrete particles stop moving as the energy is gradually dissipated by the contact friction of the FEM elements and the contact damping between the particles. Finally, the distribution range of the discrete particles and the accumulated erosion depth generated by the SPH particles to the DEM particles are obtained.

4.1 Test 1: The solid–liquid phase dam break test

To verify the numerical accuracy of the coupling analysis method in describing the solid–fluid–structure interaction and erosion entrainment behavior, the solid–liquid two-phase flow dam break test implemented by Sun et al. [50] was simulated using the established SPH–DEM–FEM coupling model (Figure 7). Since the analytical model conditions with low porosity can be satisfied by the solid DEM particles filled with fluid SPH particles, the effectiveness of the interaction forces with a series of dense particles can be verified here. We note that the overlap between the DEM particle and the SPH particle was caused by the visualization of the SPH particle and did not influence the simulation [35].

Figure 7

2D schematic of the SPH–DEM–FEM coupling model.

As shown in Figure 7, fluid SPH particles with density of 1,000 kg/m³ and dynamic viscosity of 0.001 Pa·s were orderly filled in the tank of width of 0.05 m, length in the paper direction of 0.1 m, and height of 0.1 m. Solid DEM particles were randomly distributed and filled the bottom sediment layer based on the distribution range of particle size, particle density, and required volume. In the DEM calculation, the larger the distribution range of the particle size, the greater the computational cost of searching for particle pairs. To reduce the computational costs, reference values of solid particles with a density of 2,500 kg/m³, an identical radius of 0.002 m, and a total mass of 200 g were adopted in this study. The boundary outside the fluid were all fixed by the FEM elements as a nonslip condition in addition to the moving water gate. Based on the experimental material characteristics of the two-phase flow dam break, the NULL material model combined with the Murnaghan equation of state (EOS) in LS-DYNA describing an incompressible fluid was selected for defining water, and the ELASTIC material model was selected for characterizing solid particles. The FEM section was associated with the RIGID material model with density, Young’s modulus, and Poisson ratio of 7,860 kg/m³, 100 GPa, and 0.25, respectively. Furthermore, defining the solid–liquid phase contact in numerical simulation is essential. Node-to-node contact option between the solid and fluid was selected for coupling the SPH and the DEM solvers. Node-to-surface contact option can be used for defining the contact between discrete (SPH or DEM) particles and boundary FEM element. The relevant numerical calculation parameters are presented in Table 1. When the experiment started, the water gate of the tank was lifted at a speed of 0.68 m/s, and the SPH particle drove the DEM particles to flow forward, forming the two-phase flow dam break.

Figure 8 displays the experimental and numerical results of the solid–fluid flow dam break test. The profile of the simulated fluid surface for a period of 0–0.2 s after the dam breached is demonstrated, which is almost consistent with that of the study by Sun et al. [50], especially when T = 0.2 s. It can be observed that the front flow positions of the SPH–DEM–FEM coupling model is more acceptable than the computed result of the SPH–DEM coupling model implemented by Wu et al. [35]. Based on the output of the velocity distribution cloud map (Figure 9), the present coupling model can adequately capture the peak of the dam break flow velocity and the corresponding position. As shown by the refined area in Figure 9(b) and (d), the erosion entrainment behavior of the SPH particles on the DEM particles is evident. The differential variability of velocity between the solid particles in the bottom precipitation layer and the fluid is also adequately reflected in the velocity profile because of the damping contact influence between the DEM particle and the SPH particle. Nonphysical splashing of SPH particles appears in the present simulation (see Figures 8, 9(c) and (d)), and a similar water particle splash scenario is observed from the coupled SPH–DEM calculation [41,62]. Peng et al. [62] reported that the smooth length must be addressed for particle spacing L ₀(L ₀ < h ≤ 1.5L ₀) and that when the initial smooth length is too small, problems such as numerical fracture will occur, inducing nonphysical splashing of SPH particles. The smaller the smoothing length, the smaller the constraint between SPH particles. In addition, the boundary constraints make the density variation of SPH particles more obvious near the boundary than in other regions, which considerably exaggerates the density discontinuity of boundary particles and their computational pressure [41]. The exaggerated pressure exerted by the boundary particles causes the nonphysical gap between the fluid and boundary in the SPH model. The occurrence of nonphysical particle splash is related to the small smooth length (h = 1.05L ₀) taken in the current numerical model. Figure 9(e) and (f) illustrate that the subsequent motion of the particles is mainly smooth, while the nonphysical particle splash disappears. Thus, the effect of a slight SPH particle splash on the overall motion of the fluid is acceptable.

Figure 8

Experimental and numerical results of solid–fluid flow dam break test.

Figure 9

Demonstration of the velocity distribution cloud map at each time.

In addition, the soft constraint and contact stiffness were optimized by introducing the scale factor SOFSCL and α, respectively; thus, the present coupling model eliminates the defect of extremely fast flow velocity caused by unstable motion in the SPH–DEM model [35]. The optimization result is reflected in the dimensionless relationship of the displacement in the Z-direction with time, as shown in Figure 10. Considering that the material in the present dam break model test only changes the numerical height in the Y-direction, the corresponding dynamic solid particle volume concentration ( V solid / V fluid ) can be obtained by capturing the average accumulation height of solid particles and fluid particles in the storage tank. Figure 11 illustrates the dimensionless expression of solid particle volume concentration with time. It can be seen that the dynamic change equation of the volume concentration of solid particles has a high consistency with the distribution law of the test results, which can reflect the change of the entrainment volume in the storage tank in real time.

Figure 10

Dimensionless expression of displacement in the Z-direction with time.

Figure 11

Dimensionless expression of solid particle volume concentration with time.

4.2 Test 2: liquid-solid flow in a rotating cylindrical tank

To evaluate the performance of the coupled SPH–DEM–FEM model in describing the fluid–solid interaction, the rotation model experiment of a liquid–solid flow in a cylindrical tank [71] was simulated. Both the inner diameter of the cylinder and its depth in the Z-direction were 0.1 m, and the tank rotated at 102 rpm. In the tank, the water depth was 0.05 m, and the particle bed was formed by the accumulation of solid particles with a particle size of 0.0027 m. Two tests were implemented in a previous study, in which the initial height of the particle bed was set as 0.022 m (200 solid particles in cross-section, case A) and 0.028 m (300 solid particles in cross-section, case B), respectively. To reproduce the rotation model experiments, the same physical parameters were adopted in the current simulations. Figure 12 shows the schematic of the numerical computational domain. The tank was simulated by FEM elements and defined by a RIGID material model with density, Young’s modulus, and Poisson’s ratio of 7,930 kg/m³, 200 GPa, and 0.3, respectively [72]. A combination of the null material model and Murnaghan EOS was used for the water definition. Herein, a sensitivity analysis of the contact damping coefficient of DEM–FEM η _df was performed with values of 0.05, 0.1, 0.2, 0.3, and 0.4. The initial distance between adjacent particles was 0.002 m, and the water was cross-sectionally discretized into 988 SPH particles. In addition, an equal number of DEM particles was adopted to simulate the solid particles, while the elastic material model was used to describe its physical properties. The associated elasticity coefficient and damping coefficients are presented in Table 1, and additional input parameters are summarized in Table 2.

Figure 12

Schematic diagram of the numerical computational domain.

Figure 13 illustrates the flow pattern changes of the solid–fluid flow in the rotating tank when the solid–liquid flow reaches a quasi-steady state. The SPH–DEM–FEM model intuitively reflects the transport of solid particles in water, instead of assuming the solid–liquid two-phase flow as a homogeneous mixed fluid. In the numerical cases, the particle bed profiles are consistent with those of the experimental model, both showing bilinear features. The solid–liquid flow velocity in the Y-direction at a certain time is shown in Figure 13(b). The change in the direction of the solid particles velocity occurs at both ends of the particle bed, and the flow velocity reaches the maximum value near the moving boundary. When the solid particles move with the rotating tank, the water particles drag the solid particles, causing the exfoliation and accumulation of solid particles on the surface of the particle bed. Figure 13(b) shows that fluid particles cluster at the moving boundary, which is similar to the case reported by Sun et al. [43]. The numerical accuracy near the boundary is reduced as the boundary weight function at the moving boundary depends on the contact particle size, while multiple particle sizes are in the solid–liquid two-phase flow [43]. With an increase in the wall particle to fluid particle size ratio near the moving boundary, fluid particles tend to leak or pressure tends to fluctuate at the boundary; thus, particles accumulate near the boundary. A suitable motion boundary control equation is helpful in addressing the FSI problem, which is also an important scope for future improvements of the current numerical model.

Figure 13

Experimental and numerical results of solid–fluid flow motion in the rotating tank (η _df = 0.05).

In addition, Table 3 summarizes a comparison between the simulated and experimental results of the height and width of the particle bed. In each numerical model, the simulation results deviate from the experimental results, and the relative error of the numerical results of Case A is more prominent than that of Case B; however, the degree of deviation is within an acceptable range. In general, the simulation results in this study are consistent with the experimental results, indicating that the coupled SPH–DEM–FEM model can effectively describe the interaction of fluid–solid flows with moving boundaries.

In this study, the contact force of the coupled DEM–FEM is the main driving force of the two-phase flow motion. According to equation (43), the influence of the contact damping coefficient of DEM–FEM η _df on the fluid–solid interaction cannot be ignored. Figure 14 illustrates a visual comparison of the motion form profiles of two-phase flows under different contact damping coefficients occurring at the same time as shown in Figure 13. As η _df increased, the bilinear characteristic of the particle bed gradually weakened, but the effect of exfoliation and accumulation caused by the action of the water particles on the solid particles gradually increased. In addition, the solid particle flow velocity near the moving boundary gradually increased with increasing η _df, which also caused the particle bed geometry to increase in height and decrease in width. Figure 15 illustrates the change in geometric dimensions of the particle bed as a function of the contact damping coefficient. By fitting the scatter values under a combination of height/width and η _df, a relation curve satisfying the exponential distribution law was obtained. Furthermore, by importing the experimentally measured height/width values into the fitting equation, the obtained contact damping coefficients between the solid particles and moving boundaries were 0.036 (Case A) and 0.018 (Case B), respectively. The inversion results were different from those of η _df = 0.05 (currently used in numerical calculations), which was likely the main source of the discrepancy between the numerical simulation and experimental results.

Figure 14

Effect of the contact damping coefficient of DEM–FEM on the macroscopic behavior of two-phase flows.

Figure 15

Dimensionless expression of the particle bed size with the contact damping coefficient of DEM–FEM.

As demonstrated in this section, the SPH–DEM–FEM coupling model can numerically simulate solid–liquid two-phase flow dynamic interaction, especially for capturing and illustrating the solid–fluid coupling flow erosion entrainment behavior with good applicability. Therefore, the present coupling model can provide theoretical support and technical guarantee for the application of debris flow run-out assessment in a 3D geological environment.

5.1 Yohutagawa debris flow event

For several days in October 2010, the city of Amamioshima in Kyushu, Japan, was hit by heavy rainfall due to Typhoon Megi. The rainfall considerably raised the flood level in the Yohutagawa mountain area (28°24′N, 129°32′E). Consequently, the highly weathered andesite bearing limit was broken, while a substantial amount of slurry and debris particles were drained, resulting in the Yohutagawa debris flow (Figure 16(a) and (b)). Han et al.[73] and Wu et al. [9] reported that the debris flow event occurred in a hilly mountainous area with an elevation range of 20–230 m, affecting an area of about 0.24 km² (500 m × 480 m × 210 m, as shown in Figure 16(c)). In the source zone with elevation distribution of 170–230 m, the flood water (5,843 m³) impacted loose rocks (2,854 m³) downstream along the river channel, which gradually developed into a debris flow downstream along the 750 m long river valley. At the gully outlet (50 m above sea level), the debris flow was buffered by a check dam with an average height of 2 m; 5,621 m³ debris flow material accumulated in the check dam, while the remaining 3,076 m³ material washed into the downstream village to form an alluvial fan. Based on the demonstration of disaster impact regions (as shown in Figure 17), a large portion of the debris flow energy was consumed by the buffer effect of check dam; thus, only two structures were slightly damaged, while the kindergarten and the elementary schools were not impacted by the debris flow.

Figure 16

(a) Overview of the study area: Yohutagawa mountain area located in the city of Amamioshima in Kyushu, Japan (28°24′N, 129°32′E) © Google Earth. (b) Details of the deposition area [8]. (c) 3D topographic map of the Yohutagawa debris flow.

Figure 17

Topographic contour map of the Yohutagawa debris flow distribution [73].

5.2 Parameter determination and model setup

The Amamioshima island digital terrain data with a spatial resolution of 10 m were downloaded via Google Earth and GIS. Since the exported data were in the form of irregular scattered points, the SURFER software was imported for interpolation fitting to automatically identify and capture the terrain data, draw a smooth topographic contour map, and output ordered data. By compiling the processed coordinate data, we generated point-line-surface and output the shell element. We also completed the FEM grid section of the shell element, selected the RIGID material model, and established the finite element model of 3D geological terrain. In addition, finite element modeling was performed for the check dam and the two buildings affected by the debris flow. As the structures at the site were only slightly damaged and did not fail completely, the plastic-kinematics material model was considered for definition [29].

Based on the investigations described by Wu et al. [74] and Han et al. [73], the debris flow material in the Source Zone was mainly fine-grained clay minerals with homogeneous mud formed by flooding; the erosion entrainment of slurry on loose andesite particles occurred after slurry initiation into the transfer zone. In the simulation modeling, 11,300 SPH particles with MEPH material model with an initial spatial spacing of 1 m were selected in the source zone for simulating 5,800 m³ of slurry. Meanwhile, based on the entrainment volume and erosion depth distribution [9], 1,550 DEM spherical particles with elastic material model with initial spatial spacing of 1–2 m were continuously installed in the transfer zone along the gully to simulate the dynamic entrainment process of 2800 m³ debris particles. When the slurry entered the transfer zone, the SPH–DEM coupling contact occurred, thereby simulating the slurry erosion on debris particles. As mentioned in Section 3, the contact friction coefficient is an important parameter for evaluating debris flow resistance. Since the infiltration capacity between debris particles is larger, the superporous water pressure is smaller, and the contact friction between gravel particles is larger. Thus, the sliding friction coefficient 1.0 and rolling friction coefficient 0.1 were selected based on a similar mountain debris flow case [45]. The physical and mechanical properties of the strongly weathered andesite particles in this mudslide accident are insufficiently described in the relevant literature; hence, we referred to similar andesite material parameters from Japan in the simulation modeling [75]. To simplify the calculation model, the plastic hardening modulus of slurry in the yield stress control equation (equation (10)) was set to 0, the damping effect in the contact between SPH–FEM and DEM–FEM was neglected, and only the damping effect in the SPH–DEM contact was considered. As the simulation results are not sensitive to the shear modulus and yield strength [76,77], reasonable values were taken in this study based on similar research. The definition of the coupling contact mode among SPH, DEM, and FEM as well as the selection of related contact parameters were made according to the description in Section 4. The adopted parameters of Murnaghan EOS and the material parameters are summarized in Table 4. Figure 18 illustrates the initial configuration for Yohutagawa debris flow event simulation.

Figure 18

3D debris flow initiation and particle generation (time = 0 s).

5.3 Model results

The simulated velocity variations in the Yohutagawa debris flow are shown in Figure 19. After the mud particles in the source zone were activated, they made SPH–DEM coupling contact with the loose andesite particles in the transfer zone. This phenomenon resulted in erosion entrainment behavior, which gradually developed into a discharged debris flow. As shown in Figure 19(b), the established coupled calculation model simulated the erosion entrainment process between the gully bed-sediment and debris flow. Owing to the large slope of the source zone, the debris flow accelerated. As the erosion entrainment developed, the energy dissipation gradually increased and the velocity of debris flow decreased in the transfer area. However, as the potential energy was converted to kinetic energy, the flow velocity gradually increased again. Owing to the relatively steep topography at the gully throat, the flow velocity at the throat was significantly greater than that at other areas (see Figure 19(c) and (d)); none exceeded the peak values. As shown in Figure 19(d), the debris flow started contacting the check dam, while the flow velocity decreased substantially. As the subsequent debris flow arrived, the debris flow materials were partly intercepted in the check dam, while the remaining materials crossed the check dam to impact downstream buildings. The crossing materials formed an alluvial fan (see Figures 19(e) and (f), 16(b)).

Figure 19

Velocity–time history of the simulated Yohutagawa debris flow.

Figure 20 illustrates the influence of slurry entrainment on debris particles in the transfer zone. By capturing the elevation data of debris particles in the transfer zone, the elevation difference obtained after a comparison with the initial elevation was the depth of debris particles eroded. Figure 20(a) and (b) illustrate the development of entrainment from 30 s to 60 s. As shown in Figure 17, the monitoring point (MP) at 125 m.a.s.l. was in the transfer zone. Figure 20(c) illustrates the time history of debris flow depth and accumulative erosion depth at the MP, obtained by numerical calculation. The cross-section erosion at the MP stabilized after time = 90 s and reached the maximum erosion depth of 1.05 m, which satisfied the range of the maximum erosion depth (1–1.5 m) obtained from field observations. However, the flow depth fluctuated slightly at this time, but the erosion effect on the solid particles was negligible.

Figure 20

Influence of slurry entrainment on debris particles in the transfer zone.

The entire debris flow movement lasted about 550 s, and the final debris flow material accumulation is shown in Figure 21. The maximum accumulation depth of the debris flow material was approximately 4.17 m, which is nearly consistent with postsite observations [57]. The shape of the simulated debris flow accumulation largely coincides with the shape of the measured one. In the 3D terrain, the meshing method commonly used in engineering earthwork calculations coupling with the GIS technology can be used to calculate the fluid volume based on the fluid distribution on the contour map and the fluid depth observation results. Referring to the change law between the debris volume concentration and the slope of the trench bed (i.e., Takahashi model [79]), the volume of the debris particles entrained by slurry erosion is calculated inversely. Importing the information shown in Figure 21 into the GIS system, the volumes of the debris flow in intercept area and deposition area were calculated to be 4,907 and 3,579 m³, respectively. In addition, with the help of Takahashi [79] mountain debris flow volume concentration test derivation formula (see equation (48)), the corresponding debris particle volume concentration were obtained according to the slope of intercept area (tan θ = 35/150, C _m = 0.39) and deposition area (tan θ = 20/118, C _m = 0.23), and the volume of debris particles (2736.9 m³) entrained in the slurry was obtained by inversion with the debris flow volume output from the GIS system. Compared with the calculation results, the debris particle volume of the eroded entrainment was largely consistent with the onsite investigation (2,854 m³).

Figure 21

Comparison of the simulated and measured debris flow accumulation at time = 550 s.

Figure 22 shows the time history curve of the impact force when the debris flow made contact with the structure. As shown in Figure 22(a), a comparison of the peak values of the slurry impact force at three different locations shows that the check dam has a significant effect on the energy dissipation of the debris flow. As shown in Figure 22(b), the shock and vibration of the debris particles in contact with check dam were more obvious than those of slurry. The maximum impact force generated by the debris particles at the check dam reached 21 × 10³ kN. The part of the debris particles after crossing the check dam gradually accumulated with slurry, and the peak impact force was 2.3 × 10³ kN. This value was approximately ten times different from that of the barrier dam.

Figure 22

Impact force–time history curve: (a) by slurry and (b) by debris particles.

A comparison of the resultant displacement of debris flow at time = 550 s (see Figure 23) shows that the alluvial fan scale and migration distance of debris flow without the check dam were much greater than those with the check dam. Based on the continuous distribution of particle displacement, most of the debris particles in contact with the slurry in the early stage were entrained to the middle and lower reaches of the alluvial fan. Meanwhile, the debris in the middle and later sections of the transfer zone stayed in the upper parts of the alluvial fan because of the small overall potential energy difference. Figure 24 illustrates the velocity propagation of the debris flow heads. In the case of a check dam, when the debris flow reached the check dam at time = 60 s, the flow velocity reached the peak value of 23.75 m/s, which is mostly consistent with the calculation result (23.08 m/s) obtained by Wu et al. [9] based on the GIS model. The debris flow was affected by the dam buffer, and the flow velocity dropped rapidly until the flow rate increased slightly after the debris flow crossed the dam. For the case without a check dam, when time = 65 s, the debris flow velocity peaked at 25.59 m/s.

Figure 23

Resultant displacement of the debris flow particles at time = 550 s: (a) with a check dam and (b) without a check dam.

Figure 24

Velocity propagation of the debris flow heads.