Solid-liquid two-phase flow numerical simulation method based on material point method

Abstract

The invention relates to the technical field of fluid mechanics and solid mechanics, and provides a solid-liquid two-phase flow numerical simulation method based on a material point method, which comprises the following steps of: 1, interpolating the mass and momentum of a mass point to a background grid node; 2, obtaining a background grid node liquid phase addition speed; 3, obtaining a background grid node solid phase addition speed; 4, updating the particle speed; 5, updating the node speed; 6, calculating the strain increment and the rotation increment of the solid phase of the mass point and the strain increment of the liquid phase by using the shape function and the derivative of the shape function; 7, updating the pore water pressure of the mass point; and 8, updating the position information of the mass point according to the new mass point solid phase velocity, and updating the volume, porosity and permeability coefficient of the mass point according to the solid phase strain increment of the mass point. According to the invention, solid-liquid two-phase flow numerical simulation can be well carried out.

Description

Technical Field

[0001] The present invention relates to the technical fields of fluid mechanics and solid mechanics, and in particular, to a numerical simulation method for solid-liquid two-phase flow based on the material point method. Background Art

[0002] Solid-liquid two-phase flow widely exists in the field of geological disasters, such as debris flow, rainfall-induced landslides, submarine landslides, dam seepage failures, etc. The complexity of such problems is mainly reflected in the combined action of physical phenomena such as the coupling effect between solid-phase particles and the liquid phase, and the large deformation process. These phenomena involve complex multi-scale and multi-physics field coupling problems, such as particle collisions and frictions, fluid shear effects, and dynamic interactions at the solid-liquid interface. Therefore, the numerical simulation method for solid-liquid two-phase flow has always been the focus and difficulty of current research. Traditional solid-liquid two-phase flow simulation methods mainly include the Euler method and the Lagrangian method. These two methods perform well in specific applications but also have their respective limitations:

[0003] (1) Euler Method

[0004] The Euler method uses a fixed grid to describe the motion of the liquid phase and the solid phase, and is suitable for dealing with large-scale flow problems and flow fields with strong integrity. Since the grid is fixed and does not move with the deformation of the material, this method has high efficiency in solving conservation equations and dealing with complex boundary conditions. However, when there are significant interface deformations in the flow (such as the solid-liquid mixing interface or the free surface), it is difficult for the Euler method to accurately capture the interface details, and usually an additional interface tracking algorithm needs to be introduced, which increases the complexity of the algorithm. In addition, the motion and distribution of solid-phase particles require additional model approximations, which makes the simulation results highly dependent on the empirical formula and parameter selection, and it is difficult to accurately reflect the actual physical process.

[0005] (2) Lagrangian Method

[0006] The Lagrangian method directly tracks the motion of the solid-liquid interface based on the body-fitted grid and can accurately describe the interface deformation and free surface behavior. However, when dealing with problems involving large deformations and large displacements, the body-fitted grid is prone to severe distortion, resulting in a decrease in grid quality, a significant reduction in numerical calculation efficiency, and even numerical instability. In addition, due to the grid distortion problem, the applicability of the Lagrangian method in long-term dynamic simulations is limited, especially in cases where accurate description of particle-particle interactions is required.

[0007] The Material Point Method (MPM) combines the advantages of the Eulerian method and the Lagrangian method while overcoming their limitations. It adopts a dual description of Lagrangian material points and Eulerian background grids, making it perform excellently in dealing with large deformation and multiphase flow problems. In MPM, the continuum is discretized into a set of material points. Each material point represents a material region and carries all the physical information of that region (such as mass, velocity, stress, and strain). The motion of the material points follows the Lagrangian method, which can naturally track interfaces and particle trajectories while avoiding the numerical errors of the convective terms in the Eulerian method. The background grid is used to solve the momentum conservation equation and calculate spatial derivatives but does not carry physical quantities. The characteristic of this method is that there is no relative motion between the material points and the background grid within each time step, thus avoiding the numerical difficulties of the nonlinear convective terms in the traditional Eulerian method. At the same time, after each time step, the background grid is reconstructed, thus overcoming the numerical instability problem caused by grid distortion in the Lagrangian method. The Material Point Method shows significant advantages in the field of geological disasters and has gradually been applied to the simulation calculation of geological disasters such as landslides in recent years.

[0008] With the continuous improvement of the Material Point Method algorithms (such as the Generalized Interpolation Material Point Method GIMP, the fluid phase coupling algorithm, etc.), its application scope has gradually expanded to more complex multi-physical field problems. In recent years, the improved methods based on the multi-phase Material Point Method have further enhanced its applicability in solid-liquid two-phase coupling problems, providing an efficient tool for describing complex processes such as particle flow, foundation liquefaction, and debris flow simulation. When constructing the control equations of the hydro-mechanical coupled two-phase Material Point Method, the following two formats can be adopted: the u-p format (solid phase displacement - liquid phase pressure), which ignores the acceleration of the liquid phase and was widely used in the early stage; the v-w format (solid phase velocity - liquid phase velocity), which includes all acceleration terms and can capture more complex dynamic responses and has been a research hotspot in recent years. In the multi-phase Material Point Method, the interaction between the solid phase and the liquid phase can be represented in two different ways: one is the format based on a set of material points (two-phase single-point format), that is, each material point carries all the information of the saturated porous medium, and the pore water pressure is considered an additional variable, and the motion of the material points is represented by the motion of the solid phase; the other is the format based on two sets of material points (two-phase double-point format), that is, the solid phase and the liquid phase are represented by two sets of material points respectively, which can be regarded as the superposition of two continuous media, and the material points representing the liquid phase can simulate pore water and free water. When the liquid phase flow velocity is small and the spatial variation of the solid phase volume fraction can be ignored, these two formats are equivalent and both are applicable to seepage problems. The two-phase double-point format can be extended to the case of non-laminar pore fluid, but the calculation efficiency is low, and the number of material points required in the calculation is twice that of the two-phase single-point format. The two-phase single-point format cannot simulate free water, but its expression form is simple and can be extended to the case of unsaturated pore media.

[0009] The main disadvantages of the prior art are as follows:

[0010] (1) The fixed grid cannot dynamically adapt to the movement of solid-phase particles, resulting in the interface simulation depending on additional models and empirical parameters, reducing the calculation progress of the Euler method and increasing the uncertainty of the model;

[0011] (2) In the case of large deformations, the convective grid in the Lagrangian method is prone to distortion, which will not only reduce the calculation efficiency but also may lead to the instability of the numerical solution. Summary of the Invention

[0012] The content of the present invention is to provide a numerical simulation method for solid-liquid two-phase flow based on the material point method, which preferably simulates the numerical value of solid-liquid two-phase flow.

[0013] A numerical simulation method for solid-liquid two-phase flow based on the material point method according to the present invention includes the following steps:

[0014] Step 1: Initialize the background grid and particle information. According to the position of the particles in the background grid, interpolate the mass and momentum of the particles to the background grid nodes through the shape function;

[0015] Step 2: Calculate the nodal force of the liquid phase on the background grid nodes, and combine with the boundary conditions to solve the liquid phase momentum equation, so as to obtain the liquid phase acceleration of the background grid nodes;

[0016] Step 3: Calculate the nodal force of the mixture on the background grid nodes, and combine with the boundary conditions and the nodal liquid phase acceleration to solve the mixture momentum equation to obtain the solid phase acceleration of the background grid nodes;

[0017] Step 4: Update the particle velocity according to the obtained background grid node acceleration;

[0018] Step 5: Interpolate the particle momentum back to the background grid nodes and update the nodal velocity;

[0019] Step 6: Use the shape function and its derivative to calculate the strain increment and vorticity increment of the solid phase of the particles, and the strain increment of the liquid phase;

[0020] Step 7: Update the internal force state of the solid phase of the particles according to the stress-strain relationship in the constitutive model with the strain increment, and update the pore water pressure of the particles;

[0021] Step 8: Update the position information of the particles according to the new solid phase velocity of the particles, and update its volume, porosity, and permeability coefficient from the solid phase strain increment of the particles;

[0022] Step 9: Discard the deformed background grid and initialize the calculation grid for the next calculation.

[0023] Preferably, in step 1, the formula for calculating the node mass is:

[0024]

[0025] where represents the liquid-phase node mass, represents the solid-phase node mass; the superscript k represents the calculation step, the quantity with subscript I represents the variable of the background grid node, the quantity with subscript p represents the variable carried by the particle, and N Ip is the shape function value of node I at particle p, and N p is the total number of particles; the mass of each phase of the particle is calculated according to m wp = n p ρ w V p , m sp = (1 - n p )ρ s V p to calculate, m wp is the liquid-phase particle mass, m sp is the solid-phase particle mass, V p is the particle volume, n p is the particle porosity, ρ s , ρ w are the solid-phase density and the liquid-phase density, respectively;

[0026] The formula for calculating the node momentum is:

[0027]

[0028] where represents the liquid-phase node momentum, represents the solid-phase node momentum, represents the liquid-phase velocity of the particle, represents the solid-phase velocity of the particle, and the subscript i represents the spatial coordinate component.

[0029] Preferably, in step 2, the formula for calculating the node force of the liquid phase on the background grid node is:

[0030]

[0031] where represents the external force of the node liquid phase, represents the internal force of the node liquid phase, represents the solid-liquid drag force of the node, b ip represents the gravity component of the particle, represents the pressure boundary, represents the boundary pressure, S represents the boundary region, the subscript j represents the spatial coordinate component, m represents the Kronecker function vector, p w represents the pore water pressure, Vwp represents the liquid-phase volume, k w is the permeability coefficient, γ w is the specific weight of the liquid phase;

[0032] The liquid-phase momentum equation is:

[0033]

[0034] where, represents the liquid-phase acceleration of the background grid node.

[0035] Preferably, in step 3, the calculation formula for the nodal force of the mixture on the background grid node is:

[0036]

[0037] where, represents the internal force of the nodal mixture, represents the external force of the nodal mixture, τ i represents the boundary surface force, represents the surface force boundary, m p represents the total mass of the particle, σ represents the total stress, V p represents the total volume of the particle;

[0038] The mixture momentum equation is:

[0039]

[0040] where, represents the solid-phase acceleration of the background grid node.

[0041] Preferably, in step 4, the particle velocity update formula is:

[0042]

[0043] where, Δt is the time step, N g is the total number of nodes.

[0044] Preferably, in step 5, the particle momentum is re-interpolated back to the background grid node through the following formula:

[0045]

[0046]

[0047] where, the superscript k + 1 represents the calculation step.

[0048] Preferably, in step 6, the strain increment Δε sijp and the vorticity increment Δω sijp , and the strain increment Δε of the liquid phase wijpThe calculation formula is as follows:

[0049]

[0050] Wherein, represents the solid-phase velocity of the node, represents the liquid-phase velocity of the node, and the subscripts i and j represent the spatial coordinate components.

[0051] Preferably, in step 7, the update formula for the pore water pressure p of the particle is:

[0052]

[0053] Wherein, n represents the porosity, t represents the time, v and w respectively represent the solid-phase velocity and the true liquid-phase velocity, represents the gradient operator.

[0054] Preferably, in step 8, the update formulas for the volume V p , porosity n p , permeability coefficient k w of the particle are:

[0055]

[0056] Wherein, x ip is the particle position coordinate, Δε iip is the particle volume strain increment, and C 1 is an empirical parameter.

[0057] The beneficial effects of the present invention are as follows:

[0058] (1) Improved accuracy

[0059] Traditional methods are difficult to effectively capture the large deformation behavior and solid-liquid interaction during the slope instability process caused by rainfall. The two-phase material point method proposed by the present invention can accurately describe the coupling effect between solid particles and liquid. Especially in the landslide initiation, acceleration and movement stages, the model can accurately depict the key processes such as the rapid accumulation and dissipation of pore pressure, and the formation and evolution of shear zones. The simulation results are more consistent with the actual dynamic performance of geological disasters.

[0060] (2) Enhanced applicability

[0061] The present invention does not rely on fixed grids, avoiding the problem of grid distortion in large deformation scenarios. At the same time, through an adaptable background grid update mechanism, it can flexibly handle complex terrain and boundary conditions, such as irregular slopes and non-uniform rainfall fields.

[0062] (3) Reproduction of the dynamic evolution process

[0063] The simulation results show that the method of the present invention can comprehensively reproduce the whole process of rainfall-induced landslides from initial infiltration, pore pressure increase, shear instability to landslide movement, and provide key data such as landslide body velocity and accumulation range, which can be directly used for in-depth study of the dynamic mechanism of rainfall-induced landslides. Through the sensitivity analysis of different rainfall intensities, durations and slope structure parameters, the main controlling factors of landslide occurrence can be revealed, providing theoretical support for geological disaster warning.

[0064] (4) Engineering application prospects

[0065] By simulating the movement process of rainfall-induced landslides, the present invention provides a reliable tool for disaster risk assessment and optimization of prevention and control measures. For example, in actual cases, the scale and influence range of potential landslides are predicted through numerical simulation, providing a scientific basis for the design of disaster prevention plans and significantly reducing the risk of potential economic losses and casualties.

[0066] In summary, the method of the present invention effectively makes up for the deficiencies of traditional solid-liquid two-phase flow simulation methods in the field of geological disasters, has higher calculation accuracy, wider application range and stronger practical application value, provides comprehensive and reliable technical support for geological disaster prevention and control, and has important theoretical and practical significance. Description of the drawings

[0067] Figure 1 It is a flow chart of a numerical simulation method for solid-liquid two-phase flow based on the material point method in an embodiment;

[0068] Figure 2 It is a schematic diagram of particle volume calculation in an embodiment. Detailed implementation manners

[0069] To further understand the content of the present invention, the present invention will be described in detail in combination with the drawings and embodiments. It should be understood that the embodiments are only for explaining the present invention and not for limiting it.

[0070] As Figure 1 shown, this embodiment provides a numerical simulation method for solid-liquid two-phase flow based on the material point method, which includes the following steps:

[0071] Step 1, initialize the background grid and particle information, and interpolate the mass and momentum of the particles to the background grid nodes through the shape function according to the positions of the particles in the background grid;

[0072] Step 2, calculate the nodal force of the liquid phase on the background grid nodes, and solve the liquid phase momentum equation in combination with the boundary conditions to obtain the liquid phase acceleration of the background grid nodes;

[0073] Step 3: Calculate the nodal forces of the mixture at the background grid nodes, and solve the momentum equation of the mixture by combining the boundary conditions and the nodal liquid-phase acceleration to obtain the solid-phase acceleration at the background grid nodes;

[0074] Step 4: Update the particle velocity according to the obtained acceleration at the background grid nodes;

[0075] Step 5: Interpolate the particle momentum back to the background grid nodes and update the nodal velocity;

[0076] Step 6: Calculate the strain increment and vorticity increment of the solid phase of the particle, and the strain increment of the liquid phase by using the shape function and its derivative;

[0077] Step 7: Update the internal force state of the solid phase of the particle and the pore water pressure of the particle according to the stress-strain relationship in the constitutive model based on the strain increment;

[0078] Step 8: Update the position information of the particle according to the new solid-phase velocity of the particle, and update its volume, porosity, and permeability coefficient according to the solid-phase strain increment of the particle;

[0079] Step 9: Discard the deformed background grid and initialize the computational grid for the next calculation.

[0080] In this embodiment, the control equations adopt the v - w format (solid-phase velocity - liquid-phase velocity). The calculation model is based on the two-phase single-point format. A set of particles is used to represent the saturated porous medium, and the motion of the pore medium is described by the motion of the solid phase. Each particle represents a certain volume of saturated pore medium and carries the basic physical quantities of the solid phase and the liquid phase, such as velocity, stress, strain, pore water pressure, pore water velocity, etc. The volume V represented by the particle p can be expressed as the sum of the solid-phase volume V sp and the liquid-phase volume V wp , that is, V p = V sp + V wp , as shown in Figure 2 . The motion of the particle is consistent with that of the solid phase and is described by the Lagrangian description. The motion of the liquid phase is described according to its relative motion with respect to the solid phase. Its control equations are as follows:

[0081] (A) Mass conservation equations of the solid phase and the liquid phase:

[0082]

[0083] (B) Momentum conservation equations of the mixture and the liquid phase:

[0084]

[0085] In the formula: ρ s , ρ wThey are the solid density and the liquid density respectively. v and w represent the solid velocity and the true liquid velocity respectively. n represents the porosity, and a s and a w represent the accelerations of the solid phase and the liquid phase respectively. b represents the body force, and k w is the permeability coefficient, γ w is the specific weight of the liquid phase, and t represents time. According to the principle of effective stress, σ = σ′ - p w m, where σ represents the total stress, σ′ represents the effective stress, p w represents the pore water pressure, and m is the Kronecker delta vector of the Kronecker function ([111000]T). It is stipulated that the stress tension is positive and the pore water pressure is positive.

[0086] Assume that the solid particles are incompressible, neglect the spatial variations of the density and the porosity, and according to the continuity equation of fluid mechanics, Eqs. (1) and (2) can be simplified to:

[0087]

[0088] Substitute Eq. (5) into Eq. (6) and eliminate the dn / dt term, and we can get:

[0089]

[0090] Assume that the liquid phase has the positive pressure characteristic (the liquid phase is compressible), that is:

[0091]

[0092] where K w is the bulk modulus of compressibility of the liquid phase, and finally we can get:

[0093]

[0094] For the solution of the control equations, the continuum is discretized into Lagrangian particles in the material point method, and the field variables of the particles can be interpolated from the values of their adjacent grid nodes. For example, the displacement and its derivative of particle p are expressed as:

[0095]

[0096] where, u iI is the displacement at grid node I; N Ip = N I(xp) is the shape function of particle p at grid node I; x p represents the coordinate of particle p; N Ip,j is the derivative of the shape function.

[0097] In step 1, the calculation formula for the node mass is:

[0098]

[0099] Among them, represents the mass of the liquid-phase node, represents the mass of the solid-phase node; the superscript k represents the calculation step, the quantity with subscript I represents the variable of the background grid node, and the quantity with subscript p represents the variable carried by the particle, N Ip is the shape function value of node I at particle p, N p is the total number of particles; the mass of each phase of the particle is calculated according to m wp = n p ρ w V p , m sp = (1 - n p )ρ s V p to calculate, m wp is the mass of the liquid-phase particle, m sp is the mass of the solid-phase particle, V p is the particle volume, n p is the particle porosity, ρ s , ρ w are the solid-phase density and the liquid-phase density respectively;

[0100] The calculation formula for the node momentum is:

[0101]

[0102] Among them, represents the liquid-phase node momentum, represents the solid-phase node momentum, represents the liquid-phase velocity of the particle, represents the solid-phase velocity of the particle, and the subscript i represents the spatial coordinate component.

[0103] In step 2, the calculation formula for the node force of the liquid phase on the background grid node is:

[0104]

[0105]

[0106] Among them, represents the external force of the liquid phase of the node, represents the internal force of the liquid phase of the node, represents the solid-liquid drag force of the node, represents the pressure boundary, represents the boundary pressure, S represents the boundary region, the subscript j represents the spatial coordinate component, m represents the Kronecker function vector, p w represents the pore water pressure, V wp represents the liquid-phase volume, k w is the permeability coefficient, γw is the liquid-phase specific weight, b ip represents the gravitational component of the particle.

[0107] The liquid-phase momentum equation is:

[0108]

[0109] where, represents the liquid-phase acceleration of the background grid node.

[0110] In step 3, the calculation formula for the nodal force of the mixture on the background grid node is:

[0111]

[0112] where, represents the internal force of the nodal mixture, represents the external force of the nodal mixture, τ i represents the boundary surface force, represents the surface force boundary, m p represents the total mass of the particle, σ represents the total stress, V p represents the total volume of the particle.

[0113] The mixture momentum equation is:

[0114]

[0115] where, represents the solid-phase acceleration of the background grid node.

[0116] In step 4, the particle velocity update formula is:

[0117]

[0118] where, Δt is the time step, N g is the total number of nodes.

[0119] In step 5, the particle momentum is re-interpolated back to the background grid node through the following formula:

[0120]

[0121] where, the superscript k + 1 represents the calculation step.

[0122] In step 6, the strain increment Δε sijp and the vorticity increment Δω sijp , and the strain increment Δε wijp of the liquid phase are calculated by the formula:

[0123]

[0124] where, represents the solid-phase velocity of the node, represents the liquid-phase velocity of the node, and the subscripts i and j represent the spatial coordinate components.

[0125] In step 7, the update formula for the pore water pressure p of the particle is:

[0126]

[0127] where n represents the porosity, t represents time, v and w represent the solid-phase velocity and the true liquid-phase velocity respectively, represents the gradient operator.

[0128] In step 8, the volume V of the particle p , porosity n p , permeability coefficient k w update formulas are:

[0129]

[0130]

[0131] where x ip is the particle position coordinate, Δε iip is the increment of the particle volumetric strain, and C 1 is an empirical parameter.

[0132] This embodiment adopts a modular design, and its core functions can be divided into the following several modules. Each module cooperates with each other to complete the entire simulation calculation process:

[0133] 1) Input module

[0134] This module is responsible for reading the information required for the simulation, including the number of particles, particle position information, mesh range and size, material parameters, calculation step size, boundary conditions, physical quantities to be output, selection of shape functions, etc.;

[0135] 2) Data module

[0136] This module stores the calculation data information, including the data information of various physical quantities on the particles, mesh information (mesh range and size, shape function calculation, mesh node information), and the basic physical quantity information of various material models involved in the calculation;

[0137] 3) Constitutive module

[0138] This module is the implementation of the constitutive model, used to calculate the stress-strain relationship of the particles and update the internal force state of the particles;

[0139] 4) Calculation module

[0140] This module is the core of the program and is responsible for the numerical solution based on the Material Point Method (MPM), including: ① Particle-to-Grid Mapping (P2G): Mapping the physical quantities of material points (such as mass, momentum, etc.) to the background grid for calculating the internal and external forces and interaction forces between the solid and liquid phases; ② Update of Grid Physical Information: Performing mechanical calculations on the background grid to update the momentum and velocity of the solid and liquid phases at the grid nodes; ③ Grid-to-Particle Mapping (G2P): Interpolating the updated grid information back to the particles to update the velocity and position of the material points and complete the motion simulation of the particles.

[0141] 5) Output Module

[0142] This module is responsible for outputting the simulation results. The program saves the physical states of the material points (such as position, velocity, stress, strain, etc.) at each time step, and outputs physical quantities such as the velocity, pressure, density, and temperature of the solid and liquid phases according to the requirements of the output file.

[0143] The numerical simulation method of solid-liquid two-phase flow based on the Material Point Method can reproduce the flow process of solid-liquid two-phase flow and has the following two functions:

[0144] (a) Mechanical simulation of liquid and solid phases: This method can accurately calculate the interaction forces between the solid and liquid phases, including drag force, buoyancy force, collision force, etc., considering the influence of the liquid phase on the solid phase and the reaction force of the solid phase on the liquid phase;

[0145] (b) Flow state tracking: It supports the tracking of the motion trajectories of fluids and solid particles, can output physical quantities such as the velocity, acceleration, pressure, and density of the material points during the motion process, and the output data file can display the flow states of the fluid and solid particles through visualization tools.

[0146] The method of this embodiment effectively makes up for the deficiencies of traditional solid-liquid two-phase flow simulation methods in the field of geological disasters, has higher calculation accuracy, a wider application range, and stronger practical application value, provides comprehensive and reliable technical support for geological disaster prevention and control, and has important theoretical and practical significance.

[0147] The above schematically describes the present invention and its implementation manners. This description is not restrictive, and what is shown in the drawings is only one of the implementation manners of the present invention. The actual structure is not limited thereto. Therefore, if those of ordinary skill in the art are inspired by it and design similar structural manners and embodiments without creative efforts without departing from the spirit of the present invention, they shall fall within the protection scope of the present invention.

Claims

1. A numerical simulation method for solid-liquid two-phase flow based on the material point method, characterized by: The following steps are involved: Step 1: Initialize the background grid and particle information, and interpolate the mass and momentum of the particle to the background grid node through the shape function according to the position of the particle in the background grid; Step 2: Calculate the nodal force of the liquid phase on the background grid node, solve the liquid phase momentum equation in combination with the boundary conditions, and thus obtain the liquid phase acceleration of the background grid node; Step 3, calculate the node force of the mixture on the background grid node, combine the boundary conditions and the node liquid phase acceleration, solve the momentum equation of the mixture, and obtain the solid phase acceleration of the background grid node; Step 4: Update the particle velocity according to the obtained background grid node acceleration; Step 5: interpolate the particle momentum back to the background grid node and update the node velocity; Step 6, using the shape function and its derivatives to calculate the strain increment and curl increment of the solid phase of the particle, as well as the strain increment of the liquid phase; Step 7: The strain increment updates the internal force state of the solid phase of the particle according to the stress-strain relationship in the constitutive model, and updates the pore water pressure of the particle; Step 8: Update the position information of the particle according to the new solid phase velocity of the particle, and update its volume, porosity and permeability coefficient according to the solid phase variation increment of the particle; Step 9: Discard the deformed background grid and initialize the computational grid for the next step of calculation.

2. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 1, characterized in that: In step 1, the node quality calculation formula is: in, represents the quality of the liquid node, represents the mass of the solid node; the superscript k represents the calculation step, the quantity with the subscript I represents the variable of the background grid node, and the quantity with the subscript p represents the variable carried by the particle, N Ip is the shape function value of node I at mass point p, N p is the total number of particles; the mass of each phase of particles is based on m wp =n p ρ w V p , m sp =(1-n p )ρ s V p To calculate, m wp is the mass of liquid particles, m sp is the mass of the solid phase particle, V p is the volume of the particle, n p is the particle porosity, ρ s , w are the solid phase density and the liquid phase density respectively; The formula for calculating node momentum is: in, represents the liquid phase node momentum, represents the solid phase node momentum, represents the particle-liquid phase velocity, represents the solid phase velocity of the particle, and the subscript i represents the spatial coordinate component.

3. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 2, characterized in that: In step 2, the calculation formula of the nodal force of the liquid phase on the background grid node is: in, represents the nodal liquid phase external force, represents the node liquid internal force, represents the solid-liquid drag force at the node, b ip represents the mass gravity component, represents the pressure boundary, represents the boundary pressure, S represents the boundary area, the subscript j represents the spatial coordinate component, m represents the Kronecker function vector, and p w represents the pore water pressure, V wp represents the volume of liquid phase, k w is the permeability coefficient, γ w is the liquid phase density; The momentum equation for the liquid phase is: in, Represents the liquid phase acceleration of the background grid node.

4. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 3, characterized in that: In step 3, the calculation formula of the node force of the hybrid body on the background grid node is: in, represents the internal force of the node mixture, represents the external force at the node, τ i represents the boundary force, Represents the traction boundary, m p represents the total mass of the particle, σ represents the total stress, V p represents the total volume of particles; The momentum equation of the mixture is: in, represents the solid phase acceleration of the background grid node.

5. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 4, characterized in that: In step 4, the particle velocity update formula is: Where Δt is the time step, N g is the total number of nodes.

6. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 5, characterized in that: In step 5, the particle momentum is interpolated back to the background grid node using the following formula: The superscript k+1 represents the calculation step.

7. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 6, characterized in that: In step 6, the strain increment Δε sijp and the curl increment Δω sijp , and the strain increment Δε of the liquid phase wijp The calculation formula is: in, represents the nodal solid phase velocity, represents the nodal liquid velocity, and the subscripts i and j represent the spatial coordinate components.

8. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 7, characterized in that: In step 7, the updated formula for the pore water pressure p of the particle is: Where n represents porosity, t represents time, v and w represent solid phase velocity and liquid phase true velocity respectively. Represents the gradient operator.

9. The solid-liquid two-phase flow numerical simulation method based on the material point method according to claim 8, characterized in that: In step 8, the volume V of the particle p , porosity n p , permeability coefficient k w The update formula is: Among them, x ip is the particle position coordinate, Δε iip is the particle volume strain increment, and C1 is an empirical parameter.

Images