Subsurface erosion simulation numerical method based on semi-implicit two-phase substance point method

Abstract

The invention discloses a subsurface erosion simulation numerical method based on a semi-implicit two-phase substance point method, and relates to the technical field of geotechnical engineering numerical simulation, and the method comprises the following steps: generalizing unsaturated easily subsurface-eroded accumulated soil into a two-phase multi-component medium based on a mixture theory; establishing a solid-liquid double-phase double-layer grid system; establishing a core control equation of seepage subsurface erosion in the soil body; transferring the material point information of the solid phase and the liquid phase to corresponding background grid nodes in the double-layer grid system; applying a boundary condition on a background grid of the double-layer grid system, and carrying out semi-implicit solution on the core control equation; mapping the obtained node speed to a substance point; judging whether the relative velocity of the solid phase and the liquid phase exceeds a seepage velocity threshold value or not; state variable information carried by the two-phase substance points is updated; correcting permeability and saturation; and judging and calculating whether the specified time is up or whether the soil body reaches the yield limit. The method provided by the invention has important theoretical and engineering practical values for solving the difficulty in simulating the subsurface erosion disaster of the accumulated soil.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geotechnical engineering numerical simulation, and particularly relates to a simulation numerical method for latent erosion based on a semi-implicit two-phase material point method. BACKGROUND

[0002] In the field of geotechnical engineering, the latent erosion phenomenon in accumulated soil is one of the important reasons for triggering slope instability, roadbed settlement, foundation pit gushing and other disasters. The latent erosion process is essentially a complex multi-physical field coupling process in which the erodible fine particles in the soil are stripped and migrated by the seepage of pore water, resulting in an increase in the porosity of the soil skeleton and a deterioration of the mechanical properties of the soil, and further leading to deformation and even failure of the soil. The process involves the interaction of seepage field, stress-strain field and fine particle migration field.

[0003] Traditional latent erosion simulation methods have many limitations: on the one hand, when using grid-based methods such as finite element method and finite difference method, grid distortion easily occurs in the simulation of large deformation and soil structure mutation caused by fine particle migration, resulting in a decrease in calculation accuracy or even calculation interruption; on the other hand, most numerical methods treat soil as a single-phase medium or a simple two-phase medium, ignoring the distinction between erodible fine particles and coarse particles, and the distinction between pore water and suspended fine particles, making it difficult to accurately depict the whole process of fine particle stripping-migration-soil deterioration. In addition, existing methods mostly use explicit time integration format when dealing with seepage-erosion-deformation coupling problems, and a small time step needs to be selected to ensure calculation stability, resulting in low calculation efficiency and difficulty in meeting the needs of large-scale and long-time latent erosion simulation in engineering.

[0004] Material point method (MPM) is a numerical method that combines the advantages of Lagrangian description and Eulerian description, and effectively overcomes the defects of grid-based methods in large deformation simulation by using material points to carry material properties and background grids to calculate mechanical responses, providing an ideal framework for latent erosion process simulation. However, existing latent erosion simulation research based on material point method still has some deficiencies: a special grid system for solid-liquid two-phase has not been constructed, the accuracy of solid-liquid information transmission is low; the dynamic correction of soil constitutive parameters due to fine particle mass loss is not considered, making it difficult to accurately reflect the deterioration law of soil mechanical properties; and the fully explicit or fully implicit time integration format is used, which cannot balance calculation stability and calculation 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: 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 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. 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 node velocity obtained is mapped to the material points, the material point velocity and position information are updated, the stress, strain, density and porosity of the solid phase material points are updated based on the soil constitutive model, and the porosity, volume and pore water pressure of the liquid phase material points are updated; If it is judged that the calculation reaches a specified time or the soil reaches a yield limit, the calculation is terminated, and the latent erosion evolution process, the soil deformation distribution and the fine particle migration trajectory are output, otherwise, the calculation of the next time step is entered.

[0008] The present application has the following advantages: By constructing a special grid system specially for the solid-liquid two-phase, the physical information between the solid skeleton and the pore fluid is realized in space without loss, the data distortion problem caused by the scale difference in the traditional single medium model is effectively overcome, and the simulation accuracy of the stress, strain and seepage parameters at the interface is significantly improved.

[0009] By coupling the fine particle loss mass parameter varying with the latent erosion process, the key mechanical indexes such as the stiffness, strength and permeability of the soil are realized in real-time dynamic correction. The progressive weakening process of the soil under the action of latent erosion can be accurately simulated, and the internal mechanism from the micro-particle migration to the macro-mechanical property degradation is revealed.

[0010] The semi-implicit method is used to solve the control equation, which not only ensures the stability of the solution of the complex nonlinear problem, but also effectively reduces the high calculation cost of the fully implicit method, is suitable for large-scale and long-time latent erosion evolution simulation, and significantly improves the engineering applicability of the numerical simulation. The latent erosion simulation numerical method which can accurately depict the characteristics of multi-component medium, efficiently process multi-physical field coupling and dynamically correct the mechanical parameters of the soil has important theoretical significance and engineering practical value for solving the simulation problem of the latent erosion disaster of the accumulated soil. BRIEF DESCRIPTION OF DRAWINGS

[0011] Figure 1 A latent erosion simulation numerical method based on a semi-implicit two-phase material point method is provided for the embodiment. DETAILED DESCRIPTION

[0012] The specific embodiments of the present application are described below to facilitate the understanding of the present application by those skilled in the art, but it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, it is obvious that various changes are within the spirit and scope of the present application defined and limited by the appended claims, and all the inventions utilizing the concept of the present application are within the scope of protection.

[0013] As shown in the drawings, Figure 1 In one embodiment of the present application, a latent erosion simulation numerical method based on a semi-implicit two-phase material point method includes the following steps: 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

[0014] 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).

[0015] 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).

[0016] 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.

[0017] The discrete solid-phase control equations include: 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:

[0018] 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.

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

[0020] 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 required seepage velocity, This represents the infiltration effect coefficient.

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

[0022] 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.

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

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

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

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

[0027] 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.

[0028] 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: 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.

[0029] 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.

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

[0031] 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; 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.

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

[0033] wherein,

[0034] wherein, is the velocity gradient of the region where the solid phase material point is located; is the solid phase material point corresponding momentum-related matrix; is the inertia matrix; represents the grid node index associated with the solid phase material point ; is the grid node at the first time step; corresponding solid phase material point ; is the spatial coordinate of the grid node ; is the spatial coordinate of the solid phase material point at the first time step; the superscript T represents the transpose of the matrix; is the solid phase affine matrix (affine matrix) of the material point ; the expression of the solution is:

[0035] wherein, is the mass of the solid phase material point ; is the weight of the second-order B-spline basis function, is the coordinate difference between the grid node and the material point ; is the tensor product, is the mass of the co-located grid node ; is the velocity of the co-located grid node .

[0036] The liquid phase affine matrix is decomposed into a horizontal gradient component and a vertical gradient component, and the expression is:

[0037]

[0038] wherein, represents the velocity gradient component of the liquid phase material point in the direction at the first time step; represents the momentum-related matrix of the liquid phase material point in the direction at the first time step; represents the first liquid phase material point at the time step direction inertia matrix; represents the first liquid phase material point at the time step direction velocity gradient component; represents the first liquid phase material point at the time step direction momentum related matrix; represents the first liquid phase material point at the time step direction inertia matrix.

[0039] Unlike the solid affine matrix, the fluid affine matrix needs to be adjusted for the "velocity-position separation" characteristics of the MAC grid: the horizontal velocity is stored in the direction cell face (the face perpendicular to the axis), the vertical velocity is stored in the direction cell face (the face perpendicular to the axis), and the pressure is stored in the cell center.

[0040] The solid affine matrix of the direction cell face material point is solved as: ; The solid affine matrix of the direction cell face material point is solved as: .

[0041] S5, boundary conditions are applied on 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), and the core control equation is solved semi-implicitly based on the boundary conditions to obtain the node velocity.

[0042] The incompressible Navier-Stokes equation of the liquid phase is solved by using a semi-implicit step-by-step method, and the solving process includes: The intermediate velocity is calculated, and its expression is:

[0043] In the formula, the liquid phase apparent density of the unit face , , the saturation of the unit face , the corresponding real-time porosity of the unit face , the liquid phase density; the liquid phase intermediate velocity of the unit face , the liquid phase intermediate velocity of the unit face , the time interval, the liquid phase gravity volume force of the unit face , the viscous force, the fluid viscosity, the velocity of the liquid phase at the unit face at the first time step, the velocity of the liquid phase at the unit face at the first time step, the velocity of the liquid phase at the unit face at the first time step, the unit size of the unit staggered grid specially for the liquid phase, the solid-liquid drag force; The pressure Poisson equation is constructed and solved according to the intermediate velocity to obtain the pressure field of the next time step, which is discretized at the center of the MAC grid unit and solved by the MGPCG method to avoid the assembly of the coefficient matrix and reduce the calculation cost. The pressure Poisson equation is:

[0044] In the formula, , respectively the solid phase volume fraction of the unit face , the liquid phase volume fraction of the unit face ; denote the liquid phase volume fraction at the corresponding node of the unit face ; denote the liquid phase pressure at the corresponding node of the unit face at the first time step; denote the solid phase intermediate velocity of the unit face ; The intermediate velocity is corrected by using the pressure field of the next time step, and the expression is:

[0045] In the formula, denote the corrected unit face ​the next time step liquid phase intermediate velocity.

[0046] S6, judging whether the relative velocity of the solid phase and the liquid phase exceeds a percolation velocity threshold, if yes, the water flow velocity is sufficient to overcome the adhesion between the fine particles and the solid phase skeleton and the gravity, and meets the dynamic conditions for the occurrence of latent erosion, and the latent erosion calculation is started, otherwise, the water flow only occurs percolation, no fine particle peeling, and the velocity and position information of the material points are directly updated by skipping the latent erosion calculation and entering the "regular material point updating process".

[0047] The specific steps of the latent erosion calculation are as follows: S7, solving the latent erosion control equation of the solid phase fine particles to obtain a latent erosion rate, and the expression is as follows:

[0048] In the formula, V is the latent erosion rate, V is a latent erosion sensitivity parameter, V is the solid-liquid relative velocity module of the material point, V is the liquid phase velocity of the material point, V is the solid phase velocity of the material point, and V is the volume of the material point. In the formula, V is the latent erosion rate, V is a latent erosion sensitivity parameter, V is the solid-liquid relative velocity module of the material point, V is the liquid phase velocity of the material point, V is the solid phase velocity of the material point, and V is the volume of the material point. In the formula, V is the latent erosion rate, V is a latent erosion sensitivity parameter, V is the solid-liquid relative velocity module of the material point, V is the liquid phase velocity of the material point, V is the solid phase velocity of the material point, and V is the volume of the material point. In the formula, V is the latent erosion rate, V is a latent erosion sensitivity parameter, V is the solid-liquid relative velocity module of the material point, V is the liquid phase velocity of the material point, V is the solid phase velocity of the material point, and V is the volume of the material point. According to the latent erosion rate, the overall mass conservation equation of the solid phase skeleton is solved to obtain the initial porosity and the real-time porosity. Based on the initial porosity and the real-time porosity, the permeability and the saturation are corrected by using a porosity-permeability correlation model.

[0049] In the formula, the specific method for correcting the permeability and the saturation by using the porosity-permeability correlation model is as follows: A permeability model based on the improved Kozeny-Carman theory is used for permeability correction, and the expression of the model is as follows:

[0050] In the formula, K is the current permeability coefficient, K is the initial permeability coefficient, e is the real-time porosity, is the initial porosity, and is the pore tortuosity. In the formula, K is the current permeability coefficient, K is the initial permeability coefficient, e is the real-time porosity, is the initial porosity, and is the pore tortuosity. In the formula, K is the current permeability coefficient, K is the initial permeability coefficient, e is the real-time porosity, is the initial porosity, and is the pore tortuosity.​​​​​​​​​​​​​​​​ is the initial tortuosity, is the tortuosity-porosity correlation coefficient; The saturation is corrected based on the change relationship of the saturation over time, and the expression of the relationship is:

[0051] In the formula, is the saturation, is the volume fraction of the liquid phase.

[0052] The present application greatly reduces the fitting error with the test data by the multi-parameter coupling correction of "fine particle mass fraction + pore tortuosity + saturation" (such as the internal friction angle being simultaneously affected by the fine particle loss amount, the pore tortuosity).

[0053] S7, map the obtained node velocity to the material points, update the material point velocity and position information; 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.

[0054] The soil constitutive model is a Mohr-Coulomb elastic-plasticity model, and the yield function and plastic potential function thereof are respectively:

[0055]

[0056] wherein, represents the yield function value, , is the effective principal stress, is the corrected internal friction angle, , is the initial internal friction angle, is the latent erosion influence coefficient, is the mass of the solid phase fine particles in the initial state, is the mass of the solid phase fine particles in the current state, is the corrected cohesion, , is the pore influence coefficient, , are the initial porosity and real-time porosity respectively, is the plastic potential value, is the shear dilation angle.

[0057] S8, judging whether the calculation reaches a specified time or whether the soil body reaches a yield limit (such as a solid phase stress reaching a Mohr-Coulomb yield function value), if yes, terminating the calculation and outputting a latent erosion evolution process, a soil body deformation distribution, and a fine particle migration trajectory, otherwise, returning to step S4 to enter a calculation of a next time step.

[0058] The present application realizes accurate and efficient simulation of a seepage-latent erosion-deformation coupling process in accumulated soil by designing a latent erosion simulation numerical method capable of accurately depicting multi-component medium characteristics, efficiently processing multi-physics field coupling, and dynamically correcting soil mechanical parameters, thereby providing a reliable numerical tool for latent erosion disaster prediction and prevention and control in geotechnical engineering.

Claims

1. A numerical method for simulation of latent corrosion based on a semi-implicit two-phase material point method, characterized in that, The application relates to a method for simulating the evolution of non-saturated and easily-erodible accumulated soil. The non-saturated and easily-erodible accumulated soil is generalized as a two-phase and two-point multi-component medium based on a mixing theory, and the two phases are a solid phase and a liquid phase; A solid-liquid double-layer grid system is established, including a solid-phase co-location grid and a liquid-phase marker-cell staggered grid; Core control equations of seepage-erosion-deformation coupling are established, including a solid-phase control equation, a liquid-phase incompressible Navier-Stokes equation and a liquid-phase fine particle mass conservation equation; Second-order B-spline basis functions and in-grid affine particle mapping are adopted to transfer the information of solid-phase and liquid-phase material points to corresponding grid nodes or element surfaces in the double-layer grid system; Boundary conditions are applied to the background grid of the double-layer grid system, and the core control equations are solved semi-implicitly based on the boundary conditions to obtain node velocities; It is judged whether the relative velocities of the solid phase and the liquid phase exceed a seepage velocity threshold, if yes, the erosion calculation is started, otherwise, the next step is directly entered; The obtained node velocities are mapped to the material points, the velocity and position information of the material points are updated, the stress, strain, density and porosity of the solid-phase material points are updated based on a soil constitutive model, and the porosity, volume and pore water pressure of the liquid-phase material points are updated; It is judged whether the calculation reaches a specified time or the soil reaches a yield limit, if yes, the calculation is terminated, and the erosion evolution process, soil deformation distribution and fine particle migration trajectory are output, otherwise, the calculation of the next time step is entered.

2. The method of claim 1, wherein, The solid phase comprises a coarse particle component and a potentially erodible fine particle component The liquid phase comprises a water component and a liquid fine particle component .

3. The method of claim 2, wherein, The solid-phase control equation includes: A solid-phase skeleton overall mass conservation equation: wherein porosity of the material points porosity of the material points Hamiltonian operator solid phase fine particle mass of the material points solid phase fine particle mass of the material points solid phase fine particle mass of the material points solid phase fine particle mass of the material points intrinsic density of the solid phase fine particles Jacobian of the deformation gradient of the solid phase deformation gradient tensor of the solid phase material derivative following the solid phase skeleton A solid-phase skeleton momentum conservation equation: 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 of claim 3, wherein, The expression of the liquid-phase fine particle mass conservation equation is: wherein is the liquid fine particle concentration of the material points is the fluid volume fraction of the material points is the liquid fine particle concentration of the material points is the fluid volume fraction of the material points 5. The method of claim 4, wherein, The liquid-phase incompressible Navier-Stokes equation is solved by using a semi-implicit step-by-step method, and the solving process includes: An intermediate velocity is calculated, and the expression 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; A pressure Poisson equation is constructed and solved according to the intermediate velocity to obtain a pressure field of the next time step, and the pressure Poisson equation is: wherein, , are the solid volume fraction of the cell face , the liquid volume fraction of the cell face , denotes the liquid volume fraction of the cell face at the corresponding node , denotes the liquid pressure at the corresponding node of the cell face at the time step, denotes the solid intermediate velocity of the cell face ; The intermediate velocity is corrected by using the pressure field of the next time step, and the expression is: wherein represents the modified cell face of the liquid phase intermediate velocity.

6. The method of claim 5, wherein, The specific method for transferring the information of solid-phase and liquid-phase material points to corresponding grid nodes or element surfaces in the double-layer grid system by using second-order B-spline basis functions and in-grid affine particle mapping is: Second-order B-spline basis functions are adopted to eliminate the instability of MPM particles when crossing elements, and the expression is: wherein representing a mesh node corresponding second order B-spline basis function, is a normalized distance, is a material point coordinate, is a mesh node coordinate, is a mesh size; In-grid affine particle mapping is adopted to map the momentum of each particle to the nodes of the background grid, wherein the mapping process is based on the affine matrix of the particle, so that the momentum mapped to the grid nodes contains not only the translational velocity of the particle, but also the local velocity gradient information described by the affine matrix.

7. The method of claim 6, wherein, The affine matrix includes a solid-phase affine matrix and a liquid-phase affine matrix, and the expression of the solid-phase affine matrix is: wherein wherein is the velocity gradient of the region where the solid phase material point is located; is the solid phase material point corresponding momentum related matrix; is the inertia matrix; denotes the shape function value of the grid node associated with the solid phase material point is the grid node at the time step; superscript T denotes the transpose of a matrix; corresponding solid phase material point is the spatial coordinate of the grid node ; superscript T denotes the transpose of a matrix; is the spatial coordinate of the solid phase material point at the time step; superscript T denotes the transpose of a matrix; The liquid-phase affine matrix is decomposed into a horizontal gradient component and a vertical gradient component, and the expressions are: wherein represents the velocity gradient component in the direction at the time step ; represents the moment of inertia matrix in the direction at the time step ; represents the moment of inertia matrix in the direction at the time step ; represents the velocity gradient component in the direction at the time step ; represents the moment of inertia matrix in the direction at the time step ; represents the moment of inertia matrix in the direction at the time step .

8. The method of claim 7, wherein, The soil constitutive model is a Mohr-Coulomb elastoplasticity model, and the yield function and plastic potential function are respectively: wherein, represents a yield function value, , is an effective principal stress, is a corrected internal friction angle, , is an initial internal friction angle, is a latent erosion influence coefficient, is a mass of solid-phase fine particles in an initial state, is a mass of solid-phase fine particles in a current state, is a corrected cohesion, , is a pore influence coefficient, , are an initial porosity and a real-time porosity, respectively, is a plastic potential function value, is a shear dilation angle.

9. The method of claim 8, wherein, The specific steps of the erosion calculation are: The solid-phase fine particle erosion control equation is solved to obtain an erosion rate, and the expression is: wherein, is the potential corrosion rate, is the potential corrosion sensitivity parameter, is the solid phase fine particle mass fraction of the material point , is the final solid phase fine particle mass fraction of the material point under long-term seepage, is the solid-liquid phase relative velocity module of the material point , is the liquid phase velocity of the material point , is the solid phase velocity of the material point , is the volume of the material point . The solid-phase skeleton overall mass conservation equation is solved according to the erosion rate to obtain an initial porosity and a real-time porosity. Based on the initial porosity and real-time porosity, the porosity-permeability correlation model is used to correct the permeability and saturation.

10. The method of claim 9, wherein, The specific method for correcting the permeability and saturation by using the porosity-permeability correlation model is as follows: The permeability is corrected by using a permeability model improved based on the Kozeny-Carman theory, and the expression of the model is as follows: wherein Kp is the current permeability, Kpo is the initial permeability, φ is the real-time porosity, , Δφ is the porosity increment due to fine particle induced erosion, φo is the initial porosity, C is the pore tortuosity, C0 is the initial tortuosity, is the tortuosity-porosity correlation coefficient; The saturation is corrected based on the change relationship of the saturation with time, and the expression of the relationship is as follows: The saturation is corrected based on the change relationship of the saturation with time, and the expression of the relationship is as follows: wherein is the saturation, is the volume fraction of liquid phase.

Images