Numerical simulation method for soil erosion caused by rainfall

Abstract

The invention discloses a numerical simulation method for soil erosion caused by rainfall. The method comprises the following steps: establishing a numerical model and boundary conditions of a simulation area; modeling a soil particle system by adopting a discrete element method; modeling a wind and fluid-solid interaction system by adopting computational fluid dynamics; modeling rainfall by adopting a discrete element method; performing CFD-DEM coupling solution and erosion process simulation, calculating movement of all soil particles and raindrops under related acting force in each time step, and updating positions and speeds of the movement; according to the updated particle positions, the porosity in each CFD grid unit is calculated; based on the relative velocity of the particles and the fluid; fluid acting force borne by each particle is calculated and aggregated into a source item in a CFD grid unit; substituting the calculated porosity and the source item into the model, and solving the updated fluid velocity field and pressure field; key interaction judgment and implementation are carried out; iteratively circulating until the set total simulation time is reached; and outputting related data of the soil particles in the simulation process.

Description

Technical Field

[0001] This invention relates to the field of soil erosion simulation technology, and in particular to a numerical simulation method for soil erosion caused by rainfall based on CFD-DEM coupling. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] In fields such as ecological restoration of mining areas, soil and water conservation, and prevention of geological disasters, accurate prediction and assessment of soil erosion processes caused by rainfall are crucial. High-precision numerical simulation is an effective tool for gaining a deeper understanding of erosion mechanisms, quantifying erosion effects, and guiding engineering practices.

[0004] Traditional numerical simulation methods for simulating rainfall-induced soil erosion typically employ a single theoretical framework and simplified models, which have significant limitations and distortions. 1. Model Simplification and Separation from Physical Processes: Existing methods often treat soil as a continuous medium, using methods such as the finite element method (FEM) to simulate slope stress, deformation, and stability; or they treat soil as a discrete collection of particles, using the discrete element method (DEM) to simulate particle transport; or they only use computational fluid dynamics (CFD) to simulate surface runoff. These methods fail to organically couple the "solid discrete phase (soil particles)" with the "fluid continuous phase (rainwater)," and cannot systematically reproduce the entire chain of dynamic erosion processes from raindrop splashing, particle initiation, runoff erosion to soil collapse.

[0005] 2. Lack of key interaction mechanisms: The discreteness of raindrop dynamics is ignored: Most models simplify rainfall to uniformly applied boundary conditions or continuous source terms, ignoring the discreteness, kinetic energy, and splashing and stripping effects of real raindrops on the soil surface, and failing to simulate the key mechanism by which raindrops collide with soil particles and cause their initial separation.

[0006] Simplified fluid-structure interaction: Traditional methods are difficult to accurately characterize the fluid-structure interaction forces of wind on discrete soil particles, resulting in a lack of physical basis for erosion dynamics simulation.

[0007] Changes in soil microstructure are not reflected: Existing continuous or simple discrete models are unable to describe the weakening and destruction of soil aggregates and bonding properties under raindrop impact at the particle scale.

[0008] 3. Insufficient simulation of multi-factor synergistic effects: For soil erosion under the coupled effects of complex terrain such as mining subsidence and slope changes and multiple forces such as wind and rainfall, traditional single methods or simple superposition models are difficult to reveal the erosion mechanism, which limits the application of the model in complex real-world scenarios.

[0009] Therefore, developing a numerical simulation method that can start from the micro-particle scale, couple discrete solid phases and continuous fluid phases, and completely simulate the entire process of rainfall erosion is of great scientific value and engineering significance for deepening the understanding of soil erosion mechanisms and improving prediction capabilities. Summary of the Invention

[0010] To address the aforementioned problems, this invention proposes a numerical simulation method for soil erosion caused by rainfall. This method couples computational fluid dynamics (CFD) with the discrete element method (DEM), and introduces key physical models such as soil particle bonding models, raindrop discretization models, and liquid bridge forces to construct a systematic numerical simulation framework, realizing the simulation of the entire process from raindrop splashing, particle separation, runoff erosion to soil collapse.

[0011] To achieve the above objectives, the present invention adopts the following technical solution: This invention proposes a numerical simulation method for soil erosion caused by rainfall, comprising the following steps: Step 1: Establish the numerical model and boundary conditions of the simulation region, and initialize the parameters; Step 2: Model the soil particle system using the discrete element method; Step 3: Model the wind and fluid-structure interaction system using computational fluid dynamics; Step 4: Model rainfall using the discrete element method; Step 5: Perform CFD-DEM coupled solution and erosion process simulation. Within each time step, execute the following sub-steps in sequence: Step 5.1 Based on Step 2, calculate the motion of all soil particles and raindrops under relevant forces, and update their positions and velocities; Step 5.2 Calculate the porosity within each CFD grid cell based on the updated particle positions; calculate the fluid force acting on each particle based on the relative velocity between the particles and the fluid, and aggregate them into source terms within the CFD grid cell; Step 5.3 Substitute the calculated porosity and source term into the model of Step 3 to solve for the updated fluid velocity field and pressure field; Step 5.4: Determine and implement key interactions; Step 5.5 Iterate through steps 5.1 to 5.4 until the set total simulation time is reached; Step 6: Output relevant data on soil particles during the simulation process to analyze the dynamic process, mechanism, and influencing factors of soil erosion.

[0012] As a further technical solution, in step 1, the physical properties, bonding parameters, and initial packing state of soil particles are defined; the properties and initial flow field of the fluid domain are defined; rainfall parameters are set, including raindrop size, velocity, intensity, and spatial distribution; and the simulation time step and total duration are set.

[0013] As a further technical solution, in step 2, the motion of each soil particle follows Newton's second law, where the particle translation equation is as follows:

[0014] The equation of motion for the particle is as follows:

[0015] in, m i , I i v i , ω i Particles i Mass, moment of inertia, translational velocity, and rotational velocity; f c , ij For interparticle contact force, f d,ij For damping force, f b,ij For parallel bond adhesion force, f cap,ij For capillary bridging force, f p f,i For particle-fluid interaction forces, m i g is gravity; rolling torque T ij Bonding torsional moment T b,ij Frictional torque M ij and bond moment M b,ij .

[0016] As a further technical solution, the aforementioned Based on the approximate solution of the Laplace-Young equation, it is used to calculate the capillary attraction between particles in unsaturated soil due to the presence of water. Its magnitude is related to the particle radius, surface tension, contact angle, and liquid bridge volume.

[0017] As a further technical solution, in step 2, a linear parallel bonding model is used to characterize the bonding strength between soil particles. This bonding strength can be gradually eroded and broken when raindrop impact or water flow shear force exceeds its shear threshold, thereby simulating the destruction of soil structure.

[0018] As a further technical solution, the bond fracture criterion in the linear parallel bonding model is based on the maximum tensile stress and shear stress. When the stress caused by raindrop impact or fluid shear exceeds the bond strength, the bond fractures.

[0019] As a further technical solution, in step 3, the motion of the fluid is described by the locally averaged Navier-Stokes equations, as follows: Continuity equation:

[0020] Momentum equation:

[0021] Among them, u, ρ f , p , ε These are fluid velocity, density, pressure, and porosity, respectively. τ F is the fluid viscous stress tensor; fp The opposing forces acting on a unit volume of fluid from all particles are key to the coupling of CFD and DEM, achieving bidirectional fluid-structure interaction.

[0022] As a further technical solution, in step 4, raindrops are also modeled using the discrete element method, and are regarded as spherical particles with mass, velocity, and diameter. Raindrop particles are generated from the top of the simulation domain according to the set intensity, velocity, and particle size distribution; the motion of raindrop particles also follows Newton's second law, and the forces acting on them mainly include gravity and collision forces with soil particles or other raindrops.

[0023] As a further technical solution, the specific process of step 5.4 is as follows: The DEM collision model is used to determine whether raindrop particles collide with soil particles. During the collision, the interaction force is calculated. This force may directly cause soil particles to splash or act on the parallel bonds between soil particles. When the impact stress exceeds the bond strength, it causes the bonding bonds to break, simulating the disintegration of aggregates. Through coupling force f p f , i The updated fluid velocity field generates forces that are applied to soil particles. When the fluid shear force exceeds the critical starting shear force of the particles or the residual cohesion between particles, the particles are activated and move with the fluid, simulating the scouring process. The migration of rainwater in unsaturated soil pores was simulated using a diffusion model to update the water content of soil particles; changes in the water content between particles then affect the liquid bridging force f. cap , ijThe size of the particles can be adjusted to change the effective stress and shear strength between them.

[0024] As a further technical solution, the raindrops are generated according to a given rainfall intensity, raindrop size distribution and terminal velocity, and the influence of wind is taken into account to correct their initial velocity and trajectory.

[0025] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Complete physical process and clear mechanism: This invention is the first to systematically couple the DEM raindrop model, the DEM soil particle model (including PBM bonding and liquid bridge force) with the CFD fluid model, and completely reproduce the whole chain of erosion physical process from "raindrop splashing and stripping soil particles" to "runoff scouring and transporting particles", breaking through the limitations of traditional single methods or simplified models in fragmented simulation of erosion process.

[0026] 2. Key interaction mechanisms are complete: The impact of discrete raindrops on the soil surface was accurately simulated using DEM-DEM collision, demonstrating the crucial role of raindrop kinetic energy in the initial stage of erosion.

[0027] The complex fluid-solid interaction forces between fluid (wind) and discrete soil particles were accurately characterized by bidirectional CFD-DEM coupling, providing a physical basis for particle initiation and transport.

[0028] By introducing the parallel bond model (PBM) and the liquid bridge force model, the changes in soil structural strength (dry and wet states) and their impact on erosion resistance were simulated at the particle scale, which is more consistent with the mechanical behavior of actual soil.

[0029] 3. High simulation accuracy and wide applicability: This method can take into account the synergistic effects of multiple factors such as subsidence, slope, wind force, and rainfall intensity. By adjusting the corresponding model parameters and boundary conditions, it can be flexibly applied to soil erosion simulation under different terrain and climate conditions, providing a powerful numerical tool for erosion prediction and assessment in complex scenarios.

[0030] 4. Providing new methods for mechanism research and engineering applications: This invention not only outputs macroscopic erosion data but also provides detailed information on particle-scale motion, greatly deepening our understanding of the microscopic mechanisms of erosion. Simultaneously, its simulation results can provide important data support and theoretical guidance for the design of soil and water conservation measures and the optimization of ecological restoration schemes. Attached Figure Description

[0031] Figure 1 This is an overall flowchart of the method proposed in this invention; Detailed Implementation The present invention will be further described below with reference to the embodiments.

[0032] It should be noted that the following detailed description is exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0033] It should be noted that the terminology used herein is for the purpose of describing particular implementations only and is not intended to limit the exemplary implementations of the present invention.

[0034] This embodiment discloses a numerical simulation method for soil erosion caused by rainfall, such as Figure 1 As shown, it includes the following steps: Step 1: Establish the numerical model and initialize the parameters.

[0035] Define the simulation region and boundary conditions, and define the physical properties of soil particles (density, particle size distribution, friction coefficient, etc.), bonding parameters (parallel bond strength, stiffness, etc.), and initial packing state. Define the properties of the fluid domain and the initial flow field. Set rainfall parameters, including raindrop size, velocity, intensity, and spatial distribution. Set calculation parameters such as simulation time step and total duration.

[0036] Step 2: Soil discrete phase modeling and mechanical behavior calculation.

[0037] The discrete element method (DEM) is used to model the soil particle system. The motion of each soil particle follows Newton's second law, and its translational and rotational equations are as follows: Particle translation equation:

[0038] Particle rotation equation:

[0039] in, mi , Ii v i , ωi Particles i The mass, moment of inertia, translational velocity, and rotational velocity of the particles. The forces involved include: interparticle contact force f. c , ij Damping force f d , ij Parallel bond adhesion force f b , ij (Used to simulate soil aggregates or cementation), capillary bridging force f cap , ij (Used to simulate the effect of moisture in unsaturated soil), particle-fluid interaction force f p f , i and gravitymi g. Torques include contact torque, bonding torque, etc. Specifically, a linear parallel bond model (PBM) is used to characterize the bond strength between soil particles. This bond can gradually erode and fracture when the shear force from raindrop impact or water flow exceeds its shear threshold, thus simulating soil structure failure. Furthermore, the linear parallel bond model (PBM) is used to simulate the bond between soil particles, and its bond fracture criterion is based on maximum tensile stress and shear stress. When the stress caused by raindrop impact or fluid shear exceeds the bond strength, the bond fractures. The aforementioned Based on the approximate solution of the Laplace-Young equation, it is used to calculate the interparticle capillary attraction caused by the presence of water in unsaturated soil. Its magnitude is related to the particle radius, surface tension, contact angle, and liquid bridge volume.

[0040] in, It is the harmonic mean radius of particles i and j, given by Give, It is the surface tension of the liquid. It is a dimensionless half-separation distance. It is the vector pointing from the center of particle i to j. , , and These parameters depend on the contact angle and the volume of the liquid bridge.

[0041] Step 3: Fluid continuous phase modeling and flow field calculation.

[0042] Computational fluid dynamics (CFD) is used to model the wind and fluid-structure interaction system. The motion of the fluid is described by the locally averaged Navier-Stokes equations: Continuity equation:

[0043] Momentum equation:

[0044] Among them, u, ρf , p , ε These are fluid velocity, density, pressure, and porosity, respectively. τ F is the fluid viscous stress tensor. fp The opposing forces acting on a unit volume of fluid from all particles are key to the coupling of CFD and DEM, achieving bidirectional fluid-structure interaction.

[0045] Step 4: Discretized modeling of rainfall; Raindrops are also modeled using the Discrete Element Method (DEM), treated as spherical particles with mass, velocity, and diameter. Raindrop particles are generated from the top of the simulation domain according to a set intensity, velocity, and particle size distribution. The motion of the raindrop particles also follows Newton's second law, and the forces acting on them primarily include gravity and collision forces with soil particles or other raindrops.

[0046] Step 5: CFD-DEM coupled solution and erosion process simulation. Within each time step, the following sub-steps are executed sequentially: 5.1: DEM Calculation: Based on the equations in step 2, calculate the motion of all soil particles and raindrop particles under the action of contact force, cohesion force, liquid bridge force, gravity, etc., and update their positions and velocities.

[0047] 5.2: Porosity and Coupling Force Calculation: Based on the updated particle locations, calculate the particle volume fraction (porosity) within each CFD mesh cell. ε Based on the relative velocity between the particles and the fluid, the fluid force f acting on each particle is calculated. p f , i (Mainly drag force), and aggregated into source term F within the CFD mesh cell. fp .

[0048] 5.3: CFD Calculation: The calculated porosity ε and coupling force source term F fp Substitute the fluid control equations from step 3 into the updated fluid velocity field u and pressure field. p .

[0049] 5.4: Key Interactions Identification and Realization: Raindrop-soil splashing: Using a DEM collision model, we determine whether raindrop particles collide with soil particles. During the collision, we calculate the interaction forces, which may directly cause the splashing of soil particles (especially loose surface particles) or act on the parallel bonds between soil particles. When the impact stress exceeds the bond strength, it causes the bonding bonds to break, simulating the disintegration of aggregates.

[0050] Fluid-soil erosion: via coupling force f p f , i The updated fluid velocity field generates a force (dragging force) that is applied to soil particles. When the fluid shear force exceeds the critical initiation shear force of the particles or the residual cohesion between particles, the particles are initiated and migrate with the fluid, simulating the scouring process.

[0051] Water infiltration and capillary action: A diffusion model was used to simulate the migration of rainwater in the pores of unsaturated soil, updating the water content of soil particles. Changes in the water content between particles subsequently affect the liquid bridging force f. cap , ij The size of the particles can be adjusted to change the effective stress and shear strength between them.

[0052] 5.5: Iterate through steps 5.1 to 5.4 until the set total simulation time is reached.

[0053] Step 6: Results Output and Analysis. Output data on soil particle moisture content, displacement, velocity, bond breakage, pressure, erosion rate (by statistically analyzing the residual strength of soil particles in the simulation area), and erosion morphology evolution during the simulation process. This data is used to analyze the dynamic process, mechanism, and influencing factors of soil erosion.

[0054] The liquid bridge force model is based on the approximate solution of the Laplace-Young equation and is used to calculate the capillary attraction between particles in unsaturated soil due to the presence of water. Its magnitude is related to the particle radius, surface tension, contact angle and liquid bridge volume.

[0055] The raindrops are generated according to a given rainfall intensity, raindrop size distribution, and terminal velocity, and the influence of wind is taken into account to correct their initial velocity and trajectory.

[0056] Specifically, this embodiment, based on the above method, provides a method for simulating the erosion process of a slope in a coal mining subsidence area under rainfall. The specific steps are as follows: Step 1: Model Construction and Initialization A two-dimensional computational domain is established to simulate a slope area that is 10m long and 3m high, with an initial slope of 15°. Fixed boundaries are set at the bottom and left and right ends, while the top is an open boundary.

[0057] Soil particles with a specified size distribution are generated using a DEM and then accumulated within the computational domain to form an initial slope. Linear parallel bonds are set between the soil particles, and their normal and tangential strengths and stiffnesses are assigned to simulate the aggregated structure of undisturbed soil.

[0058] Set up a CFD computational grid, with the fluids being water and air, and the initial flow field being static.

[0059] Set rainfall parameters: Set rainfall intensity and raindrop size. The final velocity of the raindrops is calculated based on the droplet size. Randomly generate raindrop DEM particles from the corresponding position at the top of the computation domain.

[0060] Step 2: Simulated execution (one time step example): DEM step: Calculates the forces acting on all soil particles and raindrop particles. The forces include: contact forces with other soil particles, cohesion forces, gravity, and coupling forces with the fluid (applied after the CFD step update).

[0061] Porosity and coupling force calculation: Based on the location of all particles, the solid volume fraction within each CFD mesh cell is calculated, and the porosity is then calculated. ε The drag force on each particle is calculated based on the difference between the current velocity of the particle and the velocity of the fluid at the previous time step.

[0062] CFD step: Apply the porosity calculated in the previous step... ε Substituting the reaction force of the particles on the fluid (coupled force source term) into the Navier-Stokes equations, we can calculate the new fluid velocity field u and pressure field. p .

[0063] Coupling force application and erosion assessment: The new flow field information calculated in the CFD step is interpolated to each particle location, updating the coupling force (drag force based on relative velocity) acting on the particle. For soil particles, it is determined whether the fluid shear force (from the coupling force) acting on them exceeds their erosion threshold. If it does, the particle's bonds will be eroded in the DEM step. The stress state of all bonds is checked; if the normal stress or shear stress of a bond exceeds its set strength, the bond breaks, indicating that the soil aggregate structure has been damaged at that location.

[0064] Moisture diffusion and capillary force update: Based on rainfall intensity and surface soil moisture status, the soil particle moisture content is updated using a diffusion model. Based on the updated moisture content, the interparticle liquid bridging forces are recalculated.

[0065] Step 3: Iteration and Output Repeat step 2 to gradually advance the simulation time.

[0066] Output soil particle location maps, velocity fields, bond fracture distribution, and landform evolution maps at different times, and analyze dynamic processes such as erosion gully development and particle transport paths.

[0067] Using the above method, this embodiment successfully simulated the complete process of soil erosion from raindrop splash and sheet erosion to gully erosion under different settlement depths (achieved by adjusting the initial slope morphology) and rainfall intensities. The numerical simulation results are consistent with the corresponding physical experimental observation trends, verifying the effectiveness of this method.

[0068] Furthermore, the method proposed in this invention can also be used in the field of aerospace propulsion and combustion chamber engineering. In this field, DEM is used to simulate high-temperature aluminum droplets or fuel droplets (as discrete liquid particles) and solid particles on the wall generated by combustion; CFD-DEM is used to simulate the fluid-structure interaction between these discrete droplets / particles and high-temperature, high-velocity combustion gas (continuous fluid phase), including the movement, collision, heat transfer and phase change of particles in the combustion gas, as well as the impact, erosion and deposition processes with the combustion chamber wall (solid structure).

[0069] The method proposed in this invention can also be used in the field of industrial spraying and surface coating. In this field, DEM is used to simulate molten or liquefied coating material droplets (discrete liquid particles) and the micro-roughness or contaminant particles (discrete solid particles) of the substrate surface; CFD-DEM is used to simulate the fluid-structure interaction between droplets and carrier gas (continuous fluid phase) to predict the acceleration, trajectory and breakup of droplets, and to calculate the spreading, splashing, solidification and their modification of surface morphology when droplets collide with the workpiece substrate (solid).

[0070] The method proposed in this invention can also be used in the field of defense and military industry (shaped charge armor penetration). In this field, DEM is used to simulate the high-speed metal jet formed by the explosion (considered as a fracture-stretchable quasi-liquid discrete body) and the armor target material (which may be discretized into solid particles to simulate damage); CFD-DEM is used to simulate the interaction between the jet and the ambient air (continuous fluid), as well as the fluid-structure interaction impact between the jet and the armor target (solid) under extreme conditions, and to reproduce the perforation process caused by fluid dynamic erosion.

[0071] The method proposed in this invention can also be used in the field of civilian industry (water jet cutting). In this field, DEM is used to simulate high-speed water jets or discrete water droplets (as liquid particles) and the processed materials (such as metals, stone, etc., which can be discrete into solid particles); CFD-DEM is used to simulate the interaction between the water jet and the surrounding air (continuous fluid), as well as the fluid-structure interaction erosion between the water jet and the material surface (solid), simulating material peeling and cutting achieved through dynamic pressure and cavitation effects.

[0072] The method proposed in this invention can also be used in the field of energy engineering (oil shaped charge perforation). In this field, DEM is used to simulate the high-speed metal jet (quasi-liquid discrete body) generated by downhole shaped charge perforation, as well as the casing, cement sheath and formation rock (which can be discretized into solid particles); CFD-DEM is used to simulate the movement of the jet in the wellbore fluid environment (continuous phase), as well as the fluid-structure interaction penetration between the jet and the multi-layer solid target, and to predict the formation and optimization of the perforation channel.

[0073] The method proposed in this invention can also be used in the field of rainfall-induced landslides. In this field, DEM is used to simulate slope soil particles (discrete solid aggregates) and raindrops (discrete liquid particles); CFD-DEM is used to simulate the fluid-structure interaction between infiltrated rainwater (as a continuous or two-phase fluid in the pores) and the soil particle skeleton, thereby reproducing the entire process of soil strength weakening, deformation and even instability and landslide.

[0074] The method proposed in this invention can also be used in the field of dam failure. In this field, DEM is used to simulate dam material particles (such as discrete solid particles such as rockfill and soil) and suspended sediment (discrete solid particles) carried by possible water flow; CFD-DEM is used to simulate the fluid-structure interaction between overflow or seepage water flow (continuous fluid) and dam particles, calculate the initiation and transport of particles by fluid drag force, and the progressive damage and failure process of dam structure caused by particle loss.

[0075] The method described in this invention is not limited to the above embodiments. By adjusting the model parameters, boundary conditions, and coupling details, it can be widely applied to various scientific research and engineering assessment scenarios related to rainfall erosion.

[0076] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A numerical simulation method for soil erosion caused by rainfall, characterized in that, Includes the following steps: Step 1: Establish the numerical model and boundary conditions of the simulation region, and initialize the parameters; Step 2: Model the soil particle system using the discrete element method; Step 3: Model the wind and fluid-structure interaction system using computational fluid dynamics; Step 4: Model rainfall using the discrete element method; Step 5: Perform CFD-DEM coupled solution and erosion process simulation. Within each time step, execute the following sub-steps in sequence: Step 5.1 Based on Step 2, calculate the motion of all soil particles and raindrops under relevant forces, and update the position and velocity of soil particles and raindrops; Step 5.2 Calculate the porosity within each CFD grid cell based on the updated particle locations; Based on the relative velocity between particles and fluid, the fluid force acting on each particle is calculated and aggregated into source terms within the CFD mesh cell. Step 5.3 Substitute the calculated porosity and source term into the model of Step 3 to solve for the updated fluid velocity field and pressure field; Step 5.4: Determine and implement key interactions; Step 5.5 Iterate through steps 5.1 to 5.4 until the set total simulation time is reached; Step 6: Output relevant data on soil particles during the simulation process to analyze the dynamic process, mechanism, and influencing factors of soil erosion.

2. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, In step 1, the physical properties, bonding parameters, and initial packing state of soil particles are defined; the properties and initial flow field of the fluid domain are defined; rainfall parameters, including raindrop size, velocity, intensity, and spatial distribution, are set; and the simulation time step and total duration are set.

3. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, In step 2, the motion of each soil particle follows Newton's second law, where the particle translation equation is as follows: The equation of motion for the particle is as follows: in, m i , I i v i , ω i Particles i Mass, moment of inertia, translational velocity, and rotational velocity; f c , ij For interparticle contact force, f d,ij For damping force, f b,ij For parallel bond adhesion force, f cap,ij For capillary bridging force, f p f,i For particle-fluid interaction forces, m i g is gravity; rolling torque T ij Bonding torsional moment T b,ij Frictional torque M ij and bond moment M b,ij .

4. The numerical simulation method for soil erosion caused by rainfall as described in claim 3, characterized in that, The aforementioned Based on the approximate solution of the Laplace-Young equation, it is used to calculate the capillary attraction between particles in unsaturated soil due to the presence of water. Its magnitude is related to the particle radius, surface tension, contact angle, and liquid bridge volume.

5. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, In step 2, a linear parallel bonding model is used to characterize the bonding strength between soil particles. This bonding strength can be gradually eroded and fractured when raindrop impact or water flow shear force exceeds its shear threshold, thereby simulating the destruction of soil structure.

6. The numerical simulation method for soil erosion caused by rainfall as described in claim 5, characterized in that, In the linear parallel bonding model, the bond fracture criterion is based on the maximum tensile stress and shear stress. When the stress caused by raindrop impact or fluid shear exceeds the bond strength, the bond fractures.

7. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, In step 3, the motion of the fluid is described by the locally averaged Navier-Stokes equations, as follows: Continuity equation: Momentum equation: Among them, u, ρ f , p , ε These are fluid velocity, density, pressure, and porosity, respectively. τ F is the fluid viscous stress tensor; fp The particle-fluid interaction force per unit volume is the key to the coupling of CFD and DEM, realizing bidirectional fluid-structure interaction.

8. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, In step 4, raindrops are also modeled using the discrete element method, and are considered as spherical particles with mass, velocity, and diameter. Raindrop particles are generated from the top of the simulation domain according to the set intensity, velocity, and particle size distribution; the motion of raindrop particles also follows Newton's second law, and the forces acting on them mainly include gravity and collision forces with soil particles or other raindrops.

9. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, The specific process of step 5.4 is as follows: The DEM collision model is used to determine whether raindrop particles collide with soil particles. During the collision, the interaction force is calculated. This force may directly cause soil particles to splash or act on the parallel bonds between soil particles. When the impact stress exceeds the bond strength, it causes the bonding bonds to break, simulating the disintegration of aggregates. Through coupling force f p f , i The updated fluid velocity field generates forces that are applied to soil particles. When the fluid shear force exceeds the critical starting shear force of the particles or the residual cohesion between particles, the particles are activated and move with the fluid, simulating the scouring process. The migration of rainwater in unsaturated soil pores was simulated using a diffusion model to update the water content of soil particles. Changes in the water content between particles thus affect the liquid bridge force f. cap , ij The size of the particles can be adjusted to change the effective stress and shear strength between them.

10. The numerical simulation method for soil erosion caused by rainfall as described in claim 1, characterized in that, The raindrops are generated according to a given rainfall intensity, raindrop size distribution, and terminal velocity, and the influence of wind is taken into account to correct their initial velocity and trajectory.

Images