A numerical method for simulation of latent corrosion based on semi-implicit two-phase material point method

Abstract

The application discloses a kind of simulation numerical method of latent erosion based on semi-implicit two-phase material point method, and relates to geotechnical engineering numerical simulation technical field, which comprises: based on mixed physics, non-saturated easy latent erosion accumulation soil is generalized as two-phase multi-component medium;Establish solid-liquid two-phase double-layer grid system;Establish the core control equation of internal seepage latent erosion of soil;The material point information of solid phase and liquid phase is transmitted to the corresponding background grid node in double-layer grid system;Boundary conditions are applied on the background grid of double-layer grid system, and the core control equation is semi-implicitly solved;The node velocity obtained is mapped to material point;Judge whether the relative velocity of solid phase and liquid phase exceeds seepage velocity threshold value;Update two-phase material point carrying state variable information;Correct permeability and saturation;Judge whether the calculation reaches specified time or whether soil reaches yield limit.The present application has important theoretical and engineering practical value for solving the simulation problem of accumulation soil latent erosion disaster.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology in geotechnical engineering, specifically to a numerical method for simulating erosion based on the semi-implicit two-phase material point method. Background Technology

[0002] In the field of geotechnical engineering, the phenomenon of erosion in embankment soil is one of the important causes of disasters such as slope instability, roadbed settlement, and foundation pit heave. The erosion process is essentially a complex multi-physics coupling process in which fine particles that can be eroded in soil are stripped and migrated by pore water under the action of pore water seepage, resulting in increased porosity of the soil skeleton, deterioration of mechanical properties, and ultimately soil deformation or even failure. It involves the interaction of seepage field, stress-strain field and fine particle migration field.

[0003] Traditional methods for simulating erosion have several limitations. Firstly, when using mesh-based methods such as the finite element method (FEM) and finite difference method (FD), mesh distortion easily occurs during simulations of abrupt changes in soil structure caused by large deformations and fine particle migration, leading to decreased computational accuracy or even computational interruption. Secondly, most numerical methods treat soil as a single-phase or simple two-phase medium, neglecting the distinction between erosive fine and coarse particles, and between pore water and suspended fine particles, making it difficult to accurately depict the entire process of fine particle stripping, migration, and soil degradation. Furthermore, existing methods often employ explicit time integration schemes when dealing with seepage-erosion-deformation coupling problems. To ensure computational stability, a small time step is required, resulting in low computational efficiency and making it difficult to meet the needs of large-scale, long-term erosion simulations in engineering projects.

[0004] The Material Point Method (MPM), a numerical method combining the advantages of Lagrangian and Eulerian descriptions, effectively overcomes the shortcomings of mesh-based methods in large deformation simulations by using material points to carry material properties and background meshes to calculate mechanical responses, providing an ideal framework for simulating erosion processes. However, existing studies on erosion simulation based on the MPM still have shortcomings: a dedicated mesh system for solid-liquid two-phase systems has not been constructed, resulting in low accuracy in solid-liquid information transfer; the dynamic correction of soil constitutive parameters by fine particle mass loss has not been considered, making it difficult to accurately reflect the deterioration law of soil mechanical properties; and the use of fully explicit or fully implicit time integration schemes cannot balance computational stability and computational efficiency.

[0005] Chinese invention patent CN110008599A provides a simulation method for water-soil coupled landslides based on a high-order dual-set two-phase material point method. Its solid-liquid co-location grid has overlapping velocity / pressure fields, which easily leads to "pressure oscillations" or "non-physical velocity abrupt changes" in the latent erosion region where "the relative velocity between solid and liquid is relatively large." This causes the calculation to diverge directly at large time steps (Δt>1.e-5), requiring only extremely small time steps and resulting in extremely low computational efficiency. Secondly, it only uses porosity to correlate strength parameters, ignoring the "strength-enhancing effect of changes in the particle size distribution of the solid skeleton (increased proportion of coarse particles) after the loss of fine particles." This results in a higher degree of soil strength weakening in the simulation than in actual experimental results, leading to a dangerously high safety factor calculation. Furthermore, it cannot simulate common engineering objects such as layered soils and heterogeneous foundations, and is only suitable for laboratory-level simulations of ideal homogeneous soils, making it difficult to directly apply to the latent erosion risk assessment of actual slopes and foundation pits. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides a numerical method for simulating latent erosion based on a semi-implicit two-phase material point method. This method solves the problems of low computational efficiency and accuracy, difficulty in accurately depicting the entire process of fine particle stripping-migration-soil deterioration, difficulty in meeting the needs of large-scale, long-term latent erosion simulation in engineering, and difficulty in dynamically correcting soil mechanical parameters.

[0007] To achieve the aforementioned objectives, the technical solution adopted by this invention is: a numerical method for simulating erosion based on a semi-implicit two-phase material point method, comprising:

[0008] Based on the mixture theory, unsaturated easily eroded embankments are generalized into two-phase, two-point, multi-component media, with the two phases being a solid phase and a liquid phase.

[0009] A dedicated two-layer mesh system for solid-liquid two-phase systems was established, including a dedicated co-location mesh for the solid phase and a dedicated marker-cell staggered mesh for the liquid phase.

[0010] Establish the core control equations for the seepage-undercutting-deformation coupling, including the solid phase control equation, the incompressible Navier-Stokes equation for the liquid phase, and the mass conservation equation for fine particles in the liquid phase.

[0011] The material point information of the solid and liquid phases is transferred to the corresponding grid nodes or cell surfaces in the two-layer grid system by using second-order B-spline basis functions and affine particle mapping within the grid.

[0012] Boundary conditions are applied to the background mesh of the two-layer mesh system, and the core control equations are solved semi-implicitly based on the boundary conditions to obtain the nodal velocities;

[0013] Determine whether the relative velocity between the solid and liquid phases exceeds the seepage velocity threshold. If it does, initiate the undercut calculation; otherwise, proceed directly to the next step.

[0014] The obtained nodal velocities are mapped to material points, and the velocity and position information of the material points are updated. Based on the soil constitutive model, the stress, strain, density and porosity of the solid phase material points are updated, and the porosity, volume and pore water pressure of the liquid phase material points are updated.

[0015] Determine whether the calculation has reached the specified time or whether the soil has reached the yield limit. If so, terminate the calculation and output the erosion evolution process, soil deformation distribution, and fine particle migration trajectory. Otherwise, proceed to the next time step for calculation.

[0016] The beneficial effects of this invention are as follows:

[0017] By constructing a dedicated mesh system specifically for solid-liquid two-phase systems, the lossless transfer of physical information between the solid skeleton and the pore fluid in space is realized, effectively overcoming the data distortion problem caused by scale differences in traditional single-medium models, and significantly improving the simulation accuracy of stress, strain and seepage parameters at the interface.

[0018] By coupling fine-particle loss mass parameters that vary with the erosion process, real-time dynamic correction of key mechanical properties such as soil stiffness, strength, and permeability can be achieved. This allows for accurate simulation of the progressive weakening process of soil under erosion, revealing the intrinsic mechanism of microparticle migration leading to macroscopic mechanical property degradation.

[0019] A semi-implicit method is employed to solve the governing equations, which effectively reduces the high computational cost of fully implicit methods while ensuring the stability of solving complex nonlinear problems. This method is suitable for large-scale, long-term simulations of erosion evolution, significantly improving the engineering applicability of numerical simulations. This numerical method for erosion simulation, capable of accurately characterizing the properties of multi-component media, efficiently handling multi-physics coupling, and dynamically correcting soil mechanical parameters, has significant theoretical and practical engineering value for solving the challenge of simulating erosion hazards in embankments. Attached Figure Description

[0020] Figure 1 The flowchart of a numerical method for simulating erosion based on a semi-implicit two-phase material point method is provided for an embodiment. Detailed Implementation

[0021] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0022] like Figure 1As shown, in one embodiment of the present invention, a numerical method for simulating erosion based on the semi-implicit two-phase material point method includes the following steps:

[0023] S1. Based on mixture theory, unsaturated, easily eroded embankments are generalized as two-phase, two-point, multi-component media, with the two phases being a solid phase and a liquid phase. Coarse-grained components are distinguished within the solid phase framework. and corrosive fine particulate components Distinguishing water components in the liquid phase and liquid fine particulate components

[0024] S2. Establish a dedicated two-layer mesh system for the solid-liquid two-phase system, including a dedicated co-location mesh for the solid phase and a dedicated marker-interlaced (MAC) mesh for the liquid phase. The dedicated co-location mesh for the solid phase is used to discretize the solid phase governing equations and store solid phase information. The dedicated MAC mesh for the liquid phase is used to discretize velocity at the element faces and pressure at the element centers, solving the incompressible Navier-Stokes equations and storing fluid information. Solid phase particles carry the mass, porosity, and constitutive parameters of the solid fine particles, while fluid particles carry the velocity, pressure, and concentration of suspended fine particles. Coupling occurs only through the mass exchange of fine particles; liquid phase particles can flow into and out of solid phase particles. The flow triggering condition is the relative velocity between the solid and liquid phases. ( (This is the seepage velocity threshold).

[0025] This invention adopts a decoupled architecture of "solid phase co-location grid + liquid phase marker-interlaced grid (MAC grid)". The MAC grid is naturally adapted to the Navier-Stokes equations for incompressible fluids, which reduces the pressure-velocity coupling error by more than 60% (significantly improving numerical stability).

[0026] S3. Establish the core control equations for the seepage-undercutting-deformation coupling, including the discrete solid phase control equation, the incompressible liquid phase Navier-Stokes equation, and the liquid phase fine particle mass conservation equation.

[0027] The discrete solid-phase control equations include:

[0028] The overall mass of the solid framework is conserved: Assuming that the solid particles (coarse particles + permeable fine particles) are incompressible, the porosity evolution is driven by both framework deformation and fine particle permeation, and its discrete form is as follows:

[0029]

[0030] In the formula, For matter points porosity, For Hamiltonian operators, For matter points The quality of solid-phase fine particles, Representing a point of matter The quality of solid-phase fine particles, This represents the inherent density of fine solid particles. Let be the Jacobian determinant of the solid deformation gradient. Let be the deformation gradient tensor of the solid phase. To follow the material derivative of the solid framework.

[0031] The governing equation for the final mass fraction of fine solid particles under long-term seepage is as follows:

[0032]

[0033] in, For reference seepage velocity, This represents the initial mass fraction of the permeable fine particles. To provide a reference for the final mass fraction of permeable fine particles at the reference seepage velocity, This represents the infiltration effect coefficient.

[0034] The momentum conservation equation for a solid framework, in discretized form, is as follows:

[0035]

[0036] In the formula, For matter points The apparent density of the solid phase, For solid density, For matter points solid-state acceleration, For effective stress, For matter points The solid volume fraction, For matter points Liquid pressure, For matter points Gravitational volume force, For matter points The solid-liquid drag force.

[0037] For solid-liquid drag force, the Ergun model is used to represent it:

[0038]

[0039] , , For fluid viscosity, Characteristic particle size of solid particles, This represents the fluid volume fraction.

[0040] The expression for the mass conservation equation for fine particles in the liquid phase is:

[0041]

[0042] In the formula, For matter points The concentration of liquid fine particles, For matter points The fluid volume fraction.

[0043] S4. Second-order B-spline basis functions and affine particle mapping within the mesh are used to transfer the material point information of the solid and liquid phases to the corresponding mesh nodes or unit surfaces in the two-layer mesh system, thereby realizing particle-mesh information transfer.

[0044] The mapping process uses the weights of the second-order B-spline basis functions. (Matter Point) For mesh nodes / cell surfaces The contribution weight is the core, divided into two directions: "matter point → mesh" and "mesh → matter point", and the solid and liquid phases follow a unified logic:

[0045] Material point → Mesh: The physical quantity of a mesh node / cell surface is equal to the physical quantities of all material points within its influence domain and their corresponding values. The weighted sum. For example, the momentum of a solid co-located grid node is calculated by multiplying the momentum of solid matter points within the influence domain by... The results are obtained by summing them up.

[0046] Mesh → Material Point: The physical quantity of a material point is equal to the physical quantities of all mesh nodes / cell surfaces within its influence domain and their corresponding physical quantities. The weighted sum. For example, the acceleration of a solid point is calculated by multiplying the mechanical response of the co-located grid nodes within the influence domain by... The results are obtained by summing them up.

[0047] The specific method is as follows:

[0048] The instability of MPM particles across units is eliminated using second-order B-spline basis functions, the expression of which is:

[0049]

[0050] In the formula, Represents grid nodes The corresponding second-order B-spline basis functions, For normalized distance, Let these be the coordinates of the material point. These are the coordinates of the grid nodes. This refers to the grid size;

[0051] Intra-mesh affine particle mapping (APIC) is used to map the momentum of each particle to nodes of the background mesh. The mapping process is based on the affine matrix of the particle, such that the momentum mapped to the mesh node includes not only the translational velocity of the particle, but also the local velocity gradient information described by the affine matrix.

[0052] Affine matrices include solid-phase affine matrices and liquid-phase affine matrices, where the expression for the solid-phase affine matrix is:

[0053]

[0054] in,

[0055] In the formula, This represents the velocity gradient in the region where the solid phase material point is located. solid phase material points The corresponding momentum correlation matrix; The inertia matrix; Represents solid matter points Associated grid node index; For the first At time step, grid nodes Corresponding solid phase material points The shape function values; For grid nodes Spatial coordinates; For the first With each time step, solid matter points Spatial coordinates; the superscript T indicates the transpose of the matrix;

[0056] Matter Point The solid affine matrix The expression is solved as follows:

[0057]

[0058] in solid phase material points quality The weights are the second-order B-spline basis functions. For grid nodes With matter point The coordinate difference For tensor product, For co-location grid nodes quality For co-location grid nodes The speed.

[0059] The liquid phase affine matrix is ​​decomposed into horizontal and vertical gradient components, and its expression is:

[0060]

[0061]

[0062] In the formula, Indicates the first At each time step, liquid phase material points exist The velocity gradient component in the direction; Indicates the first At each time step, liquid phase material points exist Momentum correlation matrix in direction; Indicates the first At each time step, liquid phase material points exist The inertia matrix of the direction; Indicates the first At each time step, liquid phase material points exist The velocity gradient component in the direction; Indicates the first At each time step, liquid phase material points exist Momentum correlation matrix in direction; Indicates the first At each time step, liquid phase material points exist The inertia matrix of the direction.

[0063] Unlike solid-phase affine matrices, fluid affine matrices need to be adjusted for the "velocity-position separation" characteristics of the MAC mesh: horizontal velocity Stored in Directional unit plane (perpendicular to) (plane of axis), vertical velocity Stored in Directional unit plane (perpendicular to) (axial surface), pressure Stored at the center of the unit.

[0064] Directional unit surface material point solid affine matrix The solution is:

[0065] ;

[0066] Directional unit surface material point solid affine matrix The solution is:

[0067] .

[0068] S5. Apply boundary conditions to the background grid of the double-layer grid system, including displacement boundary, pressure boundary, and seepage boundary (such as free boundary at the top of the slope, fixed boundary at the bottom, and seepage boundary on the side). Based on the boundary conditions, perform a semi-implicit solution to the core control equations to obtain the nodal velocities.

[0069] The semi-implicit stepwise method is used to solve the Navier-Stokes equations for incompressible liquid phases. The solution process includes:

[0070] The expression for calculating the intermediate velocity is:

[0071]

[0072] In the formula, For unit surface The apparent density of the liquid phase, , For unit surface saturation For unit surface The corresponding real-time porosity, The density of the liquid phase; For unit surface The velocity in the middle of the liquid phase, For unit surface The velocity in the middle of the liquid phase, For time intervals, For unit surface Liquid phase gravitational volume force, It is a viscous force. For fluid viscosity, For unit surface The liquid phase at the first The speed of the time step For unit surface The liquid phase in the first The speed of the time step For liquid phase specific marking - cell size of interleaved grid, For solid-liquid drag force;

[0073] The pressure Poisson equation is constructed and solved based on the intermediate velocity to obtain the pressure field at the next time step. Discretized at the center of the MAC grid cells, the solution is obtained using the MGPCG method, avoiding coefficient matrix assembly and reducing computational cost. The pressure Poisson equation is:

[0074]

[0075] In the formula, , They are unit surfaces Solid volume fraction, unit surface The liquid phase volume fraction; Representing a unit surface Corresponding node The volume fraction of the liquid phase at that point; Indicates the first At time step, unit surface corresponding nodes Liquid phase pressure at the location; Representing a unit surface The velocity in the solid phase intermediate;

[0076] The intermediate velocity is corrected using the pressure field of the next time step, and its expression is as follows:

[0077]

[0078] In the formula, Represents the corrected unit surface The next time step liquid phase intermediate velocity.

[0079] S6. Determine whether the relative velocity between the solid and liquid phases exceeds the seepage velocity threshold. If it does, the water flow velocity is sufficient to overcome the adhesion force and gravity between the fine particles and the solid phase skeleton, satisfying the kinetic conditions for the occurrence of undercutting. Enter the undercutting process and start the undercutting calculation. Otherwise, the water flow only seeps without fine particle peeling. Skip the undercutting calculation and enter the "conventional material point update process" to directly update the velocity and position information of the material points.

[0080] The specific steps for calculating erosion are as follows:

[0081] Solving the governing equation for the penetrating corrosion of solid fine particles yields the penetration rate, which is expressed as follows:

[0082]

[0083] In the formula, For the rate of erosion, For erosion sensitivity parameters, For matter points The mass fraction of solid fine particles, For long-term seepage material points The final mass fraction of solid-phase fine particles, For matter points The relative velocity modulus of solid and liquid, For matter points liquid phase velocity, For matter points solid phase velocity, For matter points Volume;

[0084] The initial porosity and real-time porosity are obtained by solving the overall mass conservation equation of the solid skeleton based on the erosion rate.

[0085] Permeability and saturation are corrected using a porosity-permeability correlation model based on initial porosity and real-time porosity.

[0086] The specific method for correcting permeability and saturation using the porosity-permeability correlation model is as follows:

[0087] A permeability correction was performed using a permeability model improved based on the Kozeny-Carman theory. The expression for this model is as follows:

[0088]

[0089] In the formula, The current permeability coefficient, The initial permeability coefficient, For real-time porosity, , The increase in porosity is caused by the erosion of fine particles. The initial porosity, For pore tortuosity, The initial tortuosity, The tortuosity-porosity correlation coefficient;

[0090] Saturation correction is performed based on the relationship between saturation and time. The expression for this relationship is:

[0091]

[0092] In the formula, For saturation, This represents the volume fraction of the liquid phase.

[0093] This invention significantly reduces the fitting error with experimental data by using multi-parameter coupling correction of "fine particle mass fraction + porosity + saturation" (e.g., the internal friction angle is simultaneously affected by fine particle loss and porosity).

[0094] S7. Map the obtained nodal velocities to material points and update the velocity and position information of the material points; update the stress, strain, density, and porosity of the solid phase material points based on the soil constitutive model; update the porosity, volume, and pore water pressure of the liquid phase material points.

[0095] The soil constitutive model is the Mohr-Coulomb elastoplastic model, with its yield function and plastic potential function as follows:

[0096]

[0097]

[0098] in, This represents the yield function value. , For effective principal stress, To correct the internal friction angle, , The initial internal friction angle, The coefficient representing the impact of erosion is denoted as . This represents the mass of the fine solid particles in the initial state. This represents the mass of the fine solid particles in their current state. To correct cohesion, , The porosity influence coefficient is... , These are the initial porosity and the real-time porosity, respectively. This is the plastic potential value. To cut the expansion angle.

[0099] S8. Determine whether the calculation has reached the specified time or whether the soil has reached the yield limit (e.g., the solid stress reaches the Mohr-Coulomb yield function value). If yes, terminate the calculation and output the erosion evolution process, soil deformation distribution, and fine particle migration trajectory. Otherwise, return to step S4 to enter the calculation of the next time step.

[0100] This invention achieves accurate and efficient simulation of the seepage-impact erosion-deformation coupling process in embankment soil by designing a numerical method for simulating imperforate erosion that can accurately characterize the properties of multi-component media, efficiently handle multi-physics field coupling, and dynamically correct soil mechanical parameters. This provides a reliable numerical tool for the prediction and prevention of imperforate erosion disasters in geotechnical engineering.

Claims

1. A numerical method for simulating erosion based on the semi-implicit two-phase material point method, characterized in that, include: Based on the mixture theory, unsaturated easily eroded embankments are generalized into two-phase, two-point, multi-component media, with the two phases being a solid phase and a liquid phase. A solid-liquid two-phase, two-layer mesh system was established, including a solid phase co-location mesh and a liquid phase marker-cell staggered mesh; Establish the core control equations for the seepage-undercutting-deformation coupling, including the solid phase control equation, the incompressible Navier-Stokes equation for the liquid phase, and the mass conservation equation for fine particles in the liquid phase. The material point information of the solid and liquid phases is transferred to the corresponding grid nodes or cell surfaces in the two-layer grid system by using second-order B-spline basis functions and affine particle mapping within the grid. Boundary conditions are applied to the background mesh of the two-layer mesh system, and the core control equations are solved semi-implicitly based on the boundary conditions to obtain the nodal velocities; Determine whether the relative velocity between the solid and liquid phases exceeds the seepage velocity threshold. If it does, initiate the undercut calculation; otherwise, proceed directly to the next step. The obtained nodal velocities are mapped to material points, and the velocity and position information of the material points are updated. Based on the soil constitutive model, the stress, strain, density and porosity of the solid phase material points are updated, and the porosity, volume and pore water pressure of the liquid phase material points are updated. Determine whether the calculation has reached the specified time or whether the soil has reached the yield limit. If yes, terminate the calculation and output the erosion evolution process, soil deformation distribution, and fine particle migration trajectory. Otherwise, proceed to the next time step for calculation. The semi-implicit stepwise method is used to solve the incompressible Navier-Stokes equations for the liquid phase. The solution process includes: The expression for calculating the intermediate velocity is: In the formula, For unit surface The apparent density of the liquid phase, , For unit surface saturation For unit surface The corresponding real-time porosity, The density of the liquid phase; For unit surface The velocity in the middle of the liquid phase, For unit surface The liquid phase at the first The speed of the time step The time interval between adjacent time steps. For unit surface Liquid phase gravitational volume force, It is a viscous force. The dynamic viscosity of the liquid phase. For unit surface The liquid phase at the first The speed of the time step For unit surface The liquid phase in the first The speed of the time step For liquid phase labeling - cell size of interleaved mesh, For unit surface Solid-liquid drag force; Based on the intermediate velocity, construct and solve the pressure Poisson equation to obtain the pressure field at the next time step. The pressure Poisson equation is as follows: In the formula, , They are unit surfaces Solid volume fraction, unit surface The liquid volume fraction, Representing a unit surface Corresponding node The liquid volume fraction at that point Indicates the first At time step, unit surface corresponding nodes The liquid phase pressure at that location, Representing a unit surface The velocity in the solid phase intermediate; The intermediate velocity is corrected using the pressure field of the next time step, and its expression is as follows: In the formula, Represents the corrected unit surface The velocity in the middle of the liquid phase.

2. The method according to claim 1, characterized in that, The solid phase includes coarse-particle components. and corrosive fine particulate components The liquid phase includes an aqueous component. and liquid fine particulate components .

3. The method according to claim 2, characterized in that, The solid-phase governing equations include: The overall mass conservation equation for the solid framework is as follows: In the formula, For matter points porosity, For Hamiltonian operators, For matter points The quality of solid-phase fine particles, Representing a point of matter The quality of solid-phase fine particles, This represents the inherent density of fine solid particles. Let be the Jacobian determinant of the solid deformation gradient. Let be the deformation gradient tensor of the solid phase. To follow the material derivative of the solid framework; Momentum conservation equation for solid framework: In the formula, For matter points The apparent density of the solid phase, For solid density, For matter points solid-state acceleration, For effective stress, For matter points The solid volume fraction, For matter points Liquid pressure, For matter points Gravitational volume force, For matter points The solid-liquid drag force.

4. The method according to claim 3, characterized in that, The expression for the mass conservation equation for fine particles in the liquid phase is: In the formula, For matter points The concentration of liquid fine particles, For matter points The fluid volume fraction.

5. The method according to claim 4, characterized in that, The specific method for transferring the material point information of the solid and liquid phases to the corresponding grid nodes or element surfaces in the two-layer grid system using second-order B-spline basis functions and in-mesh affine particle mapping is as follows: The instability of MPM particles across units is eliminated using second-order B-spline basis functions, the expression of which is: In the formula, Represents grid nodes The corresponding second-order B-spline basis functions, For normalized distance, Let these be the coordinates of the material point. These are the coordinates of the grid nodes. This refers to the grid size; A affine particle mapping within the mesh is used to map the momentum of each particle to a node of the background mesh. The mapping process is based on the affine matrix of the particle, so that the momentum mapped to the mesh node includes not only the translational velocity of the particle, but also the local velocity gradient information described by the affine matrix.

6. The method according to claim 5, characterized in that, Affine matrices include solid-phase affine matrices and liquid-phase affine matrices, where the expression for the solid-phase affine matrix is: in, In the formula, This represents the velocity gradient in the region where the solid phase material point is located. solid phase material points The corresponding momentum correlation matrix; The inertia matrix; Represents solid matter points Associated grid node index; For the first At time step, grid nodes Corresponding solid phase material points The shape function values; For grid nodes Spatial coordinates; For the first With each time step, solid matter points Spatial coordinates; the superscript T indicates the transpose of the matrix; The liquid phase affine matrix is ​​decomposed into horizontal and vertical gradient components, and its expression is: In the formula, Indicates the first At time step, liquid phase substance points exist The velocity gradient component in the direction; Indicates the first At time step, liquid phase substance points exist Momentum correlation matrix in direction; Indicates the first At time step, liquid phase substance points exist The inertia matrix of the direction; Indicates the first At time step, liquid phase substance points exist The velocity gradient component in the direction; Indicates the first At time step, liquid phase substance points exist Momentum correlation matrix in direction; Indicates the first At time step, liquid phase substance points exist The inertia matrix of the direction.

7. The method according to claim 6, characterized in that, The soil constitutive model is a Mohr-Coulomb elastoplastic model, with its yield function and plastic potential function as follows: in, This represents the yield function value. , For effective principal stress, To correct the internal friction angle, , The initial internal friction angle, The coefficient representing the impact of erosion is denoted as . This represents the mass of the fine solid particles in the initial state. This represents the mass of the fine solid particles in their current state. To correct cohesion, , The porosity influence coefficient is... , These are the initial porosity and the real-time porosity, respectively. The value of the plastic potential function. To cut the expansion angle.

8. The method according to claim 7, characterized in that, The specific steps for calculating erosion are as follows: Solving the governing equation for the penetrating corrosion of solid fine particles yields the penetration rate, which is expressed as follows: In the formula, For the rate of erosion, For erosion sensitivity parameters, For matter points The mass fraction of solid fine particles, For long-term seepage material points The final mass fraction of solid-phase fine particles, For matter points The relative velocity modulus of solid and liquid, For matter points liquid phase velocity, For matter points solid phase velocity, For matter points Volume; The initial porosity and real-time porosity are obtained by solving the overall mass conservation equation of the solid skeleton based on the erosion rate. Permeability and saturation are corrected using a porosity-permeability correlation model based on initial porosity and real-time porosity.

9. The method according to claim 8, characterized in that, The specific method for correcting permeability and saturation using the porosity-permeability correlation model is as follows: A permeability correction was performed using a permeability model improved based on the Kozeny-Carman theory. The expression for this model is as follows: In the formula, The current permeability coefficient, The initial permeability coefficient, For real-time porosity, , The increase in porosity is caused by the erosion of fine particles. The initial porosity, For pore tortuosity, The initial tortuosity, The tortuosity-porosity correlation coefficient; Saturation correction is performed based on the relationship between saturation and time. The expression for this relationship is: In the formula, For saturation, This represents the volume fraction of the liquid phase.

Images