Semi-analytical mpm-dem coupled numerical simulation method for interaction between debris flow and structures
The MPM-DEM semi-analytical coupling method was used to achieve efficient simulation of the collaborative impact and blockage process of boulders and driftwood in debris flows. This method solves the problems of low computational efficiency and insufficient accuracy in existing technologies, and provides a scientific analysis of the stress on structures and the probability of blockage. It is applicable to debris flow disaster prevention engineering design.
Patent Information
- Application Number
- CN202511485010.3
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-10-17
- Publication Date
- 2025-12-26
- Estimated Expiration
- 2045-10-17
AI Technical Summary
Existing technologies struggle to simultaneously simulate the synergistic impact effects and blockage processes of boulders and driftwood at an acceptable computational cost. Furthermore, traditional models suffer from low computational efficiency and poor stability under complex conditions, failing to accurately describe the interactions between boulders, driftwood, and structures in debris flows.
The MPM-DEM semi-analytical coupling method is adopted, which treats mud fluid as a continuous medium and boulders and driftwood as discrete solids. Through grid-particle information exchange and bidirectional transmission of particle-grid coupling force, the synchronous solution of fluid-solid multiphase coupling is realized. The contact force is accurately calculated by combining the Hertz-Mindlin model. The semi-analytical method is used to handle fluid-solid contact force to ensure numerical stability and efficiency.
It achieves high-fidelity modeling of the synergistic impact mechanism of boulders and driftwood, with stable and efficient calculations. It can graphically output the mechanical properties and blockage probability of structures, providing a scientific basis for the design of protective structures and is suitable for large-scale complex working condition simulation.
Smart Images

Figure CN120951729B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application relates to the field of geotechnical engineering disaster simulation, in particular to a semi-analytical MPM-DEM coupling numerical simulation method for interaction between debris flow and structure, which adopts coupling of a material point method (MPM) and a discrete element method (DEM) (referred to as MPM-DEM coupling) to realize efficient simulation analysis on the interaction process between a block stone-containing debris flow (solid-liquid two-phase medium) and a structure. BACKGROUND
[0002] Debris flow is a destructive mountain disaster in which a mixture of loose deposits and water moves at high speed along a gully or slope under the action of gravity, and the fluid usually contains block stones (boulders and pebbles) and tree trunks (driftwood) with a large particle size span. The high-speed moving block stones will exert an impulse impact load on protective structures such as retaining dams, row piles and bridge piers, with small local contact area and high peak force, which can easily trigger secondary damage such as concrete erosion and steel bar buckling. The greater the geometric shape coefficient of the block stone, the more significant the impact amplification effect, which has become a key control factor in impact-resistant design of structures. Standing trees in mountain forests play a dual role of "blocker" and "participant" in the debris flow disaster chain: unbroken standing trees can increase the roughness of the slope and dissipate part of the kinetic energy, while broken driftwood under the action of strong flow is converted into a rigid discrete solid moving with the flow and is easily jammed with block stones at narrow gullies, fence dams or bridge orifices to form "bridge arch blockage". The blockage body rapidly raises the potential energy of the upstream water-mud-stone mixture and may suddenly break, significantly amplifying the downstream flood peak and the disaster chain impact range. A large number of field investigations and physical model tests have shown that the cooperative blockage of block stones and driftwood is one of the key trigger factors for the expansion of secondary debris flow disasters.
[0003] Traditional empirical formula and continuum / particle flow model which equivalent debris flow as single-phase viscous flow or pure granular flow, are difficult to describe the coupling behavior of block stone multi-scale shape effect, driftwood posture-locked dynamics and mud-water phase non-Newtonian yield flow. The continuum method (such as classical CFD) is limited to Euler grid and is difficult to accurately depict the discrete solid impact, blockage and large displacement contact; pure DEM can capture rigid body collision, but it is difficult to efficiently handle the viscous-plastic-shear thinning properties of viscous slurry fluid. In order to make up for their respective defects, researchers began to try to couple particle method and continuum method, among which the coupling framework based on material point method (MPM) and discrete element method (DEM) is concerned because it can track fluid continuous field and large particle discrete body at the same time. However, the fully analytical MPM-DEM coupling still has obvious disadvantages under complex real working conditions: it needs to accurately integrate all interactions (fluid-solid, solid-solid, solid-structure) at each time step and perform global contact search, resulting in limited calculation step, huge storage and CPU / GPU resource consumption; when facing millions of particles, large-scale domain decomposition and high-frequency impact scenarios, parallel efficiency often decreases, convergence is difficult and even calculation is unstable. The existing public work focuses on the single coupling scene of "block stone-mud" or "driftwood-water flow", and has not yet been able to consider the double solid phase of block stone and driftwood, mud fluid and its complex interaction with engineering blocking structure at an acceptable computational cost. In view of this technical gap, the present application proposes a semi-analytical simulation method based on MPM-DEM, which effectively avoids the high computational overhead of the fully analytical framework while ensuring coupling accuracy, to finely reproduce the block stone and driftwood collaborative impact-blocking-breach process in debris flow, and provide efficient and reliable numerical support for protective structure design and disaster chain risk assessment. SUMMARY
[0004] The main purpose of the present application is to overcome the limitations of the prior art and provide a numerical simulation method that can realistically reproduce the block stone and driftwood collaborative impact effect in block stone-containing debris flow. This method can simultaneously simulate the complex interaction between mud flow and two types of discrete solids (block stone, driftwood) and protective structures, and quantitatively analyze the impact load and blocking process of the structure, providing a scientific basis for debris flow disaster prevention engineering.
[0005] The specific technical solutions of the present application are as follows:
[0006] A semi-analytical MPM-DEM coupled numerical simulation method for the interaction between debris flows and structures is proposed. Fine-grained debris flow is treated as a continuous medium and discretized using the material point method (MPM). Spherical particles (rock slabs) and clustered particles (driftwood) are used as the solid phase, modeled using the discrete element method (DEM). Structures and boundaries are embedded in the computational domain as rigid walls. Simultaneous solution of the fluid-solid multiphase coupling is achieved through mesh-particle (P2G / G2P) information exchange and bidirectional transfer of particle-mesh coupling forces. The fluid-solid contact forces are calculated semi-analytically, ensuring numerical stability and efficiency. The interactions between DEM particles and the particle-wall interactions are accurately calculated based on the Hertz-Mindlin model.
[0007] Specifically, the following steps are included:
[0008] Step 1: Determine the model size and calculation parameters through on-site investigation, measurement and evaluation, including the tank inclination, fluid volume, stone content and driftwood size;
[0009] Step 2: Establish the model;
[0010] Debris flow: constructed using MPM;
[0011] Stones: Constructed using SPHERE particles from the DEM;
[0012] Driftwood: Constructed using CLUMP particles from the DEM;
[0013] Structure: Constructed using wall elements from the DEM;
[0014] Step 3: Treat the fine-grained mud and water phase in the debris flow as a continuous medium and represent them discretely using MPM;
[0015] The continuous phase (mud or water) is governed by the mass conservation equation and the momentum conservation equation:
[0016] mass conservation equation:
[0017] ,
[0018] Momentum conservation equation:
[0019] ,
[0020] In the formula, The density of the material; For material velocity; For stress tensor; The coupling force between DEM and MPM; stress tensor Determined by the constitutive relations of the fluid; For time; This is the acceleration due to gravity.
[0021] Step 4: Treat large boulders and driftwood in the debris flow as discrete solids and represent them using a Discrete Image Model (DEM). The motion of DEM particles (stones, driftwood) follows Newton's second law:
[0022] DEM particle translation equation:
[0023] ,
[0024] DEM particle rotation equation:
[0025] ,
[0026] In the formula, DEM particle mass; The acceleration of DEM particles; For all and The neighboring objects (particles + walls); Contact force of DEM particles; The coupling force between DEM and MPM; The moment of inertia of the particle; This refers to the particle's angular velocity. This refers to the contact torque between particles; This is the coupling torque;
[0027] Damping coefficient Using Abraham's piecewise empirical formula:
[0028] ,
[0029] The first section is the Stokes region, where viscosity is dominant and particle Reynolds number is high. Less than 1; the middle section follows the Abraham empirical formula, with a smooth transition and a particle Reynolds number. The range is between 1 and 1000; the third segment has a Reynolds number greater than 1000. It tends to a constant of 0.44.
[0030] Second law of damping (dragging force):
[0031] ,
[0032] in, It is the fluid density. It is the particle radius. It is porosity. It is relative velocity. It is the porosity-dependency index:
[0033] ,
[0034] Step 6, update the position information and velocity information of the solid phase and liquid phase particles;
[0035] Step 7, visualize the current simulation effect, output physical information, enter the next calculation or end the current calculation.
[0036] Compared with the prior art, the present application has the following beneficial effects:
[0037] (1) Realize high-fidelity modeling of block stone-wood impact mechanism: the present application successfully solves the problem that the traditional single model cannot simultaneously simulate the collision of discrete large particles and viscous flow through the MPM-DEM coupling model of double solid phase + mud fluid, so as to truly reproduce the physical process of block stone and wood jointly blocking and transmitting impact in front of the structure, and improve the cognitive accuracy of the multi-phase action mechanism of debris flow.
[0038] (2) Stable and efficient calculation, suitable for complex large-scale working conditions: compared with pure DEM or traditional CFD method, the present application can complete three-dimensional large-scale simulation containing millions of fluid particles and tens of thousands of particles in a reasonable time, and the numerical stability and accuracy are also verified to be excellent.
[0039] (3) Rich result parameters, which can directly serve engineering design: the present application can not only obtain mechanical indexes such as stress time history and peak value of the structure, but also output functional indexes such as wood-block stone blocking probability, blocking duration and discharge capacity, so as to evaluate the performance of the protective structure from multiple angles. These data and indexes provide a scientific basis for quantitative design and optimization of debris flow protection engineering, and make up for the shortcomings of relying on experience in the past.
[0040] (4) Strong universality, large expansion space: the modeling and implementation of the present method have universality, and can be extended and applied to different types of debris flow protection research as needed. For example, a non-Newtonian mud constitutive model can be further introduced to simulate the specific mud flow slurry characteristics, or the structure can be set as elastic-plastic to evaluate the deformation response of the structure itself. This reflects the broad application prospect and innovation of the present method in the field of mountain disaster numerical simulation. BRIEF DESCRIPTION OF DRAWINGS
[0041] Figure 1 The flow chart of the MPM-DEM coupling simulation numerical method of the present application.
[0042] Figure 2 The modeling flow chart of the present application;
[0043] Figure 3 The simulation diagram of the embodiment. DETAILED DESCRIPTION
[0044] The specific technical solutions of the present application are illustrated by combining with the embodiments.
[0045] The present application provides a semi-analytical MPM-DEM coupled numerical simulation method for the interaction between debris flow and structures, as shown in the flowchart. It includes: Figure 1
[0046] I. Modeling, as shown in Figure 2
[0047] Debris flow: constructed by MPM;
[0048] Stone: constructed by SPHERE particles in DEM;
[0049] Driftwood: constructed by CLUMP particles in DEM;
[0050] Structure: constructed by wall units in DEM.
[0051] II. Numerical implementation:
[0052] 1. Discrete representation of fine-grained mud and water in debris flow as a continuous medium by MPM;
[0053] (1) The continuous phase (mud or water) is controlled by the mass conservation equation and the momentum conservation equation:
[0054] Mass conservation equation:
[0055] ,
[0056] Momentum conservation equation:
[0057] ,
[0058] (2) Numerical integration scheme:
[0059] MPM explicit integration:
[0060] The explicit integration process of MPM is usually divided into three stages: "particle to grid (P2G)", "grid solving" and "grid to particle (G2P)".
[0061] a. P2G: particle mass and momentum converge to grid nodes
[0062] For each material point , its mass and momentum are distributed to each node according to the interpolation function :
[0063] ,
[0064] In the formula, is the material pointmass; mass; for a material point at time velocity; for a material point at node a quadratic B-spline interpolation function.
[0065] b, mesh solving: node force and momentum update
[0066] Calculate the force on each node
[0067] ,
[0068] Then update the momentum:
[0069] ,
[0070] Then update the momentum:
[0071] c, G2P: mesh particle velocity and position update
[0072] Velocity update (PIC / FLIP)
[0073] ,
[0074] where control the weight of PIC and FLIP
[0075] Position update
[0076] ,
[0077] Strain rate tensor:
[0078] Strain increment:
[0079] Stress update:
[0080] 2, the large particles in the debris flow and driftwood are regarded as discrete solids, and are discretely represented by DEM.
[0081] Among them, the driftwood is fine and complex in shape, and is modeled by a clump (CLUMP) particle composed of multiple spherical particles to reflect the true geometric shape and inertia characteristics of the driftwood; the block stone is simplified as an equivalent spherical particle for modeling. The protective structures (such as stone dams, row piles, etc.) and the channel boundaries are embedded in the simulation domain in the form of rigid wall elements, and do not participate in the movement but can interact with the fluid and particles.
[0082] (1) DEM particle (stone, driftwood) motion satisfies Newton's second law:
[0083] DEM particle translational equation:
[0084] ,
[0085] DEM particle rotation equation:
[0086] ,
[0087] (2) Hertz-Mindin contact:
[0088] Normal force
[0089]
[0090] where , ,
[0091] Damping: hunt-crossley:
[0092] ,
[0093] where e is the normal restitution coefficient;
[0094] Tangential force :
[0095] ,
[0096] Tangential stiffness:
[0097] Contact radius:
[0098] Mindlin model that two elastic balls on the contact surface will form a circular contact area with a radius of a. At this time, the tangential stiffness is proportional to the size of the circular contact area, which is derived from the elastic theory . The greater the stiffness, the higher the tangential spring modulus, and the stronger the resistance to sliding.
[0099] Cumulative tangential displacement (update of ):
[0100] ,
[0101] : Cumulative tangential displacement, which is the history of the tangential spring's deformation variable retained during the contact between two particles;
[0102] : Tangential relative velocity, calculated as:
[0103] ,
[0104] where is the relative velocity of the two particles at the contact point, is the normal unit vector;
[0105] Trial tangential force:
[0106] ,
[0107] ,
[0108] Moment of force:
[0109] Tangential moment of force : , is the particle radius
[0110] Rolling moment of force :
[0111] ,
[0112] ,
[0113] : Cumulative roll angle;
[0114] ,
[0115] : Roll stiffness;
[0116] ,
[0117] : Rolling resistance coefficient, same as friction coefficient;
[0118] DEM velocity Verlet algorithm:
[0119] For each DEM particle , the classical velocity-Verlet scheme is used for explicit integration to update the position and velocity with second order accuracy:
[0120] a. Position half-step update:
[0121]
[0122] b. Intermediate velocity update:
[0123]
[0124] c. Force recalculation:
[0125]
[0126] d. Velocity full-step update:
[0127]
[0128] where, is the total force at the two time steps; is the same time step for MPM, ensuring coupling synchronization:
[0129]
[0130] where, is the grid size; is the sound speed; is the particle mass; is the contact stiffness.
[0131] 3. MPM-DEM semi-analytical coupling algorithm:
[0132] To balance the large-scale three-dimensional calculation efficiency of the multi-phase system of mud-block stone-driftwood, the invention adopts a semi-analytical coupling strategy: the fluid-solid (mud ⇄ block stone / driftwood) contact force is approximated by resistance-pressure gradient, and the solid-solid and solid-structure still maintain the Hertz-Mindlin contact precision; the coupling process is completed synchronously in a single explicit time step. The core steps are as follows:
[0133] Fluid-solid resistance: for any DEM particle p (radius , velocity ), its relative velocity with the fluid is:
[0134]
[0135] Particle Reynolds number :
[0136]
[0137] where, is the particle diameter, is the voidage correction, considering the increase of effective fluid velocity after group effect;
[0138] Damping coefficient Abraham's piecewise empirical formula is adopted:
[0139]
[0140] The first segment is the stokes region (viscosity dominated); the middle segment is the Abraham empirical formula, smoothly transitioning to a constant of 0.44 at high Reynolds numbers.
[0141] Quadratic drag law (drag force):
[0142]
[0143] where, is the fluid density, is the particle radius, is the porosity, is the relative velocity, and the exponent is the porosity-dependent exponent:
[0144]
[0145] Gaussian kernel weight:
[0146]
[0147] Volume fraction mapping, enabling conservative interpolation between particle-grid and grid-particle;
[0148] Particle-to-grid solid volume mapping:
[0149]
[0150] Updating porosity:
[0151] Particle local porosity:
[0152] Force on particles:
[0153] Reaction on fluid grid nodes:
[0154] Efficient approximate calculation of mud-boulder two-phase interaction (buoyancy, drag, pressure field) through empirical resistance + fluid stress divergence term.
[0155] III. Results output and structural response analysis:
[0156] The simulation results of the embodiment are shown in Figure 3 .
[0157] The method of the present application can capture the full-time response of the impact load of the protective structure during simulation. For example, by recording the contact force of each time step on the rigid structure, the force-time curve of the structure can be obtained. Based on the curve, key indicators such as impact peak force and total impact momentum (impulse) can be extracted for structure impact resistance design. In addition, when the driftwood and block stone are blocked in front of the structure, the method can calculate the blocking rate, discharge capacity and other parameters representing the interception efficiency of the structure to evaluate the functional performance of the structure. By changing the input working condition parameters, the inventors can also use the method to design multi-factor simulation tests: that is, a series of simulations are performed for different combinations of driftwood geometric size, driftwood quantity, block stone content, mud viscosity and structure spacing, and the influence law of each factor on the impact response and blocking behavior of the structure is compared and analyzed. The present application can thus provide a systematic parameter sensitivity analysis method for engineering. For example, it can reveal the influence trend of the ratio of driftwood length to structure spacing on the blocking probability, the influence degree of driftwood quantity and debris flow concentration on the maximum impact force, and the like, thereby providing a quantitative basis for the optimal design of protective structures.
[0158] In summary, the present application realizes precise simulation of the interaction process between debris flow containing block stones and driftwood and structures through innovative MPM-DEM semi-analytical coupling modeling and efficient calculation. The method can not only reproduce the formation mechanism of driftwood-block stone blocking structures, but also obtain the force response data of protective structures under different working conditions, filling the gap in the prior art and having important engineering application value.
Claims
1. A semi-analytical MPM-DEM coupled numerical simulation method for interaction between debris flow and structure, characterized in that, The process comprises the following steps: The fine-grained slurry is discretized as a continuum medium by using the material point method (MPM); meanwhile, the solid phase is modeled by using the discrete element method (DEM) with spherical particles and clustered particles; and the structure and the boundary are embedded in the calculation domain in the form of rigid walls; The fluid-solid multiphase coupling is synchronously solved through the bidirectional transmission of grid-particle information exchange and particle-grid coupling force; the fluid-solid contact force is calculated in a semi-analytical manner; and the DEM particles and the DEM particle-wall contact are precisely calculated based on the Hertz-Mindlin model; The semi-analytical coupling algorithm is used between the MPM and the DEM; Damping coefficient The Abraham piecewise empirical formula was used: , First segment is stokes regime, viscous dominated, particle Reynolds number Less than or equal to 1; middle segment is Abraham empirical, smooth transition, particle Reynolds number Between 1 and 1000; third segment Reynolds number greater than or equal to 1000, Tends to constant 0.44; Secondary damping law, drag force: , wherein, is the fluid density, is the particle radius, is the porosity, is the relative velocity, is the porosity-dependent exponent: 。 2. The semi-analytical MPM-DEM coupled numerical simulation method for debris flow and structure interaction according to claim 1, characterized in that, The process comprises the following steps: Step 1: The model size and the calculation parameters are determined through field investigation, measurement and evaluation, including the water tank inclination, the fluid volume, the stone content and the driftwood size; Step 2: The model is established; MPM is used to build the debris flow; SPHERE particles in the DEM are used to build the stones; CLUMP particles in the DEM are used to build the driftwood; The wall unit in the DEM is used to build the structure; Step 3: The fine-grained slurry and the water phase in the debris flow are regarded as a continuum medium, which is discretely represented by using the MPM; The continuum phase is controlled by the mass conservation equation and the momentum conservation equation: The mass conservation equation: , The momentum conservation equation: , wherein is the material density; is the material velocity; is the stress tensor; is the coupling force between DEM-MPM; stress tensor is determined by the constitutive relation of the fluid; is the time; is the gravitational acceleration; Step 4: The large-grained stones and the driftwood in the debris flow are regarded as discrete solids, which are discretely represented by using the DEM, and the DEM particle movement satisfies the Newton's second law: The DEM particle translational equation: , The DEM particle rotational equation: , where, is the DEM particle mass; is the DEM particle acceleration; is the total contact force on the particle (particle + wall); is the total contact force on the particle (particle + wall); is the DEM particle contact force; is the coupling force between DEM-MPM; is the particle moment of inertia; is the particle angular velocity; is the particle contact torque; is the coupling torque; Step 5: The semi-analytical coupling algorithm is used between the MPM and the DEM; Step 6: The position information and the velocity information of the solid phase and the liquid phase particles are updated; Step 7: The current simulation effect is visualized, the physical information is output, and the next calculation or the current calculation is ended.
Citation Information
Patent Citations
Sandy soil seepage failure simulation method and device and storage medium
CN114564899A
Method for simulating water and soil surface loss / underground leakage process in karst area
CN115544911A