Mixed representation and concurrent solution method for soil-rock mixture

By combining the material point method and the discrete element method, the problems of microscopic mechanisms and large computational load in solving soil-rock mixtures are solved, achieving efficient and accurate reproduction of the mechanical properties of soil-rock mixtures, which is suitable for natural disaster early warning in geotechnical engineering.

CN116467959BActive Publication Date: 2026-05-01INST OF ROCK & SOIL MECHANICS CHINESE ACAD OF SCI
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
INST OF ROCK & SOIL MECHANICS CHINESE ACAD OF SCI
Filing Date
2023-03-28
Publication Date
2026-05-01

AI Technical Summary

Technical Problem

Existing methods for solving soil-rock mixtures cannot effectively reproduce microscopic mechanisms and involve large computational loads. They also have low versatility and are difficult to meet the early warning needs of natural disasters such as debris flows, debris flows, and landslides in geotechnical engineering.

Method used

Soil is modeled using the Material Point Method (MPM) and rocks are modeled using the Discrete Element Method (DEM). The interaction between soil and rocks is described by a pre-defined contact model. An MPM-DEM characterization scheme is constructed, the virtual radius and coupled contact force are calculated, equilibrium equations are constructed, and the concurrent solution of material points and discrete elements is achieved by combining the MUSL stress update scheme and the FLIP velocity update scheme.

Benefits of technology

While ensuring computational efficiency, it accurately reproduces the microscopic mechanism and mechanical properties of soil-rock mixtures, reduces computational resource requirements, improves computational efficiency, simplifies code implementation, and is suitable for large-scale engineering problems.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN116467959B_ABST
    Figure CN116467959B_ABST
Patent Text Reader

Abstract

The application discloses a kind of soil and rock mixture mixed representation and concurrent solving method, comprising: after determining the soil and rock block of target soil and rock mixture, soil is modeled with material point method, and the material point corresponding to soil is generated, rock block is modeled with discrete element method, and the discrete unit corresponding to rock block is generated;The interaction between soil and rock block is described with a preset contact model;The virtual radius of each material point is calculated;And determine the material point in contact with discrete unit, marked as virtual particle;The coupling contact force of virtual particle and discrete unit is calculated using discrete element contact model;According to coupling contact force, construct balance equation, and calculate the acceleration on each background grid node in material point method, and the acceleration and angular acceleration on each discrete unit;According to the acceleration on background grid node, the acceleration of discrete unit and the angular acceleration of discrete unit, the real-time speed and displacement of material point and discrete unit are calculated.
Need to check novelty before this filing date? Find Prior Art

Description

Characterization and Concurrent Solution Methods for Soil-Rock Mixtures Technical Field

[0001] This invention relates to the field of geotechnical engineering technology, and in particular to a method for characterizing and solving soil-rock mixtures concurrently. Background Technology

[0002] Soil-rock mixtures are a widespread material in geotechnical engineering. Their heterogeneity and nonlinearity result in highly complex mechanical properties. Accurately reproducing and predicting the mechanical properties of soil-rock mixtures is crucial for early warning of natural disasters such as debris flows, clastic flows, and landslides. Generally, there are two main methods for solving soil-rock mixture problems: continuous methods and discontinuous methods. In continuous methods, such as finite element method (FEM), finite difference method (FD), and smoothed particle hydrodynamics, the soil-rock mixture is treated as a continuous medium, and its mechanical behavior is described by a constitutive model. This can greatly simplify the problem and reduce the computational load, but the interaction between soil and rock cannot be considered, and some microscopic mechanisms cannot be reproduced. In discontinuous methods, such as discrete element method (DEM), each soil particle and rock can be modeled in great detail, which can account for microscopic mechanisms to some extent. However, for large-scale engineering problems, the computational load increases significantly. Furthermore, some microscopic parameters are difficult to calibrate, which is one of the important factors limiting the application of the DEM in practical engineering. Summary of the Invention

[0003] The main objective of this invention is to propose a mixed characterization and concurrent solution method for soil-rock mixtures, aiming to solve the problems of existing soil-rock mixture solution methods being unable to reproduce microscopic mechanisms, having large computational loads, and low versatility.

[0004] To achieve the above objectives, this invention proposes a method for characterizing and concurrently solving soil-rock mixtures, the method comprising:

[0005] After determining the soil and rocks of the target soil-rock mixture, the soil is modeled using the material point method, and material points corresponding to the soil are generated. The rocks are modeled using the discrete element method, and discrete elements corresponding to the rocks are generated.

[0006] The interaction between the soil and the rocks is described using a preset contact model, and an MPM-DEM characterization scheme for the target soil-rock mixture is constructed.

[0007] Calculate the imaginary radius of each of the aforementioned material points;

[0008] Based on the virtual radius, determine the material points that are in contact with the discrete unit and mark them as virtual particles;

[0009] The coupling contact force between the virtual particle and the discrete unit is calculated using a discrete element contact model.

[0010] Based on the coupled contact force, an equilibrium equation is constructed, and the acceleration on each background grid node, as well as the acceleration and angular acceleration on each discrete element, are calculated in the material point method.

[0011] Based on the acceleration of the background grid nodes, the acceleration of the discrete unit, and the angular acceleration of the discrete unit, the real-time velocity and position parameters of the material point and the discrete unit are calculated.

[0012] Optionally, calculating the imaginary radius of each of the said material points includes:

[0013] Using formula Calculate the imaginary radius of each of the aforementioned material points, where, Let P be the imaginary radius of the P-th material point. Let k be the volume of the p-th material point. p Let be the porosity of the p-th material point.

[0014] Optionally, the discrete element method is used to calculate the coupling contact force between the virtual particle and the discrete unit, including:

[0015] Using formula Calculate the resultant normal force, where, Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. p is the normal force exerted by the q-th material point on the P-th discrete unit, and pContact is the neighborhood list of particles in contact with the P-th discrete unit.

[0016] Using formula Calculate the tangential resultant force, where, Let be the resultant tangential force from the corresponding material point acting on the P-th discrete unit. p is the tangential force exerted by the q-th material point on the P-th discrete unit, and pContact is the neighborhood list of particles in contact with the P-th discrete unit.

[0017] The resultant normal force and the resultant tangential force are combined to form the coupled contact force.

[0018] Optionally, based on the coupled contact force, an equilibrium equation is constructed, and the acceleration at each background grid node in the material point method, as well as the acceleration and angular acceleration at each of the discrete elements, are calculated, including:

[0019] Constructing equilibrium equations using the material point method in, The resultant internal force at the i-th node in the material point method. The resultant external force at the i-th node in the material point method is... Let be the mass of the i-th node in the material point method. Let be the acceleration of the i-th node in the material point method;

[0020] Solve the equilibrium equations of the material point method to obtain the accelerations at each background grid node in the material point method.

[0021] Optionally, solving the equilibrium equations of the material point method to obtain the acceleration at each background grid node and the acceleration and angular acceleration at each discrete element in the material point method further includes:

[0022] Using formula Calculate the resultant internal force, where N IP For shape functions, Physical force applied to a point of matter from the outside world. Let be the surface force acting on the i-th node in the material point method. Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. Let be the resultant tangential force from the corresponding material point acting on the Pth discrete unit;

[0023] Using formula Calculate the net external force, where, Let p be the stress at the p-th material point. Let B be the volume of the p-th material point. Ip The gradient of the shape function;

[0024] Calculate the acceleration at each background grid node based on the resultant internal force and the resultant external force.

[0025] Optionally, based on the coupled contact force, constructing equilibrium equations and calculating the acceleration at each background grid node in the material point method, as well as the acceleration and angular acceleration at each of the discrete elements, further includes:

[0026] Constructing equilibrium equations using the discrete element method as well as in, Let be the resultant normal force from other discrete elements acting on the p-th discrete element. Let be the resultant tangential force from other discrete elements acting on the p-th discrete element. Let be the body force acting on the p-th discrete unit. Let be the rotational damping force acting on the p-th discrete element. Let p be the mass of the p-th discrete unit. Let I be the acceleration of the p-th discrete unit. p Let be the moment of inertia of the p-th discrete element. Let be the angular acceleration of the p-th discrete unit. Let n be the particle size of the p-th discrete unit. pq Let δ be the unit normal vector at the contact point. n,pq Let p be the normal stacking amount between the p discrete units and the q-th discrete unit (or material point). Let be the force exerted on the p-th discrete element by the q-th discrete element;

[0027] Solve the equilibrium equations of the discrete element method to obtain the acceleration and angular acceleration on each discrete element.

[0028] Optionally, the real-time parameters of the material point include the material point velocity, and the real-time parameters of the discrete unit include the discrete unit velocity;

[0029] Based on the acceleration of the background grid nodes, the acceleration of the discrete element, and the angular acceleration of the discrete element, the real-time parameters of the material point and the discrete element are calculated, including:

[0030] The velocity of the material points is updated in real time using the FLIP velocity update format, wherein the FLIP velocity update format is as follows: as well as Let be the neighborhood list of the adjacent nodes of the p-th material point, and let Δt be the time step of the material point. For the p-th matter point at t n+1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time, For the p-th matter point at t n+1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time;

[0031] The velocity of the discrete unit is updated in real time using a preset velocity update format, wherein the preset velocity update format is as follows: as well as For the p-th matter point at t n+1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time, For the p-th matter point at t n+1 / 2 acceleration at any moment For the p-th matter point at t n-1 / 2 Acceleration at any moment.

[0032] Optionally, the real-time parameters of the material point include the position of the material point, and the real-time parameters of the discrete unit include the position of the discrete unit;

[0033] Based on the acceleration of the background grid nodes, the acceleration of the discrete element, and the angular acceleration of the discrete element, the real-time parameters of the material point and the discrete element are calculated, including:

[0034] use The update format updates the position of the material point in real time, wherein, For the p-th matter point at t n+1 The coordinates of the matter point at time [time]. For the p-th matter point at t n-1 The coordinates of the matter point at that moment;

[0035] The positions of the discrete units are updated in real time using a preset position update format, wherein the preset position update format is as follows: as well as For the p-th discrete unit at t n+1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n-1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n+1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n-1 The coordinates of the discrete unit at time t.

[0036] Optionally, based on the solution results, the stress, velocity, and position of the material points, as well as the velocity and position of the discrete elements, are updated according to a preset update format, further including:

[0037] Map the momentum parameters of the material points onto the background mesh;

[0038] Based on the momentum parameters, the velocities at the grid nodes are solved using a preset equation, wherein the preset equation is: as well as

[0039] The strain at the material point is updated using a strain update format, wherein the strain update format is set to...

[0040] The spinor at the material point is updated using a spinor update format, wherein the spinor update format is set as follows:

[0041] Based on the strain and the spinor, the stress at the material point is updated using a preset constitutive model.

[0042] Optionally, before the steps of modeling the soil and rocks of the target soil-rock mixture using the material point method to obtain the material points corresponding to the soil and modeling the rocks using the discrete element method to obtain the discrete elements corresponding to the rocks, the method further includes:

[0043] Obtain the boundary particle size of the target soil-rock mixture;

[0044] Based on the defined boundary particle size, the soil and rocks of the target soil-rock mixture are determined.

[0045] In the technical solution of this invention, particles smaller than the boundary particle size are considered as soil, modeled using the Material Point Method (MPM), and their mechanical behavior is described using macroscopic constitutive models. The particle size of the material points can be appropriately increased to reduce the total number of particles and save computational resources. Particles larger than the boundary particle size are considered as rocks, modeled using the Discrete Element Method (DEM) to track the movement and contact behavior of individual rocks. Each material point is assigned a virtual radius for contact detection, and material points in contact with DEM elements are treated as virtual particles to calculate the coupling contact force. The coupling contact force is treated as an external boundary condition, applied to both the material points and the discrete elements, and their respective equilibrium equations are constructed. The equilibrium equations of the material points and the discrete elements are solved, and the stress and velocity of the material points, as well as their positions, are updated using the MUSL stress update scheme and the FLIP velocity update scheme, respectively. For the discrete elements, their velocity and position can be directly updated. It should be noted that the Material Point Method is a numerical method that combines Lagrange and Eulerian descriptions. Compared to mesh-based numerical methods, the material point method (MPM) overcomes problems such as mesh distortion; compared to other meshless methods, the MPM also boasts excellent computational efficiency. The discrete element method (DEM) has significant advantages in handling multi-body contact, making it highly suitable for dealing with rocks in soil-rock mixtures and the contact between soil and rocks. This invention couples the MPM and DEM, using the MPM to simulate soil and the DEM to simulate rocks, fully leveraging the advantages of both. Given the continuous nature of the MPM, macroscopic constitutive models can be used to describe the mechanical response of the soil, without needing to consider the contact behavior between individual soil particles. While maintaining the overall mechanical properties, the size of the material points can be appropriately increased to reduce the total number of particles in the model and improve computational efficiency; treating each rock as a separate DEM element without fine subdivision also reduces the total number of particles in the model and improves computational efficiency. Simultaneously, the contact between rocks and between rocks and soil is preserved, which can reflect the microscopic mechanisms to some extent. Finally, in the material point-discrete element coupling framework given in this invention, the calculations of the material point method and the discrete element method can be performed concurrently. Apart from calculating the coupling contact force, the other steps of the material point method and the discrete element method do not affect each other, which simplifies the code implementation difficulty and improves the calculation efficiency. Attached Figure Description

[0046] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.

[0047] Figure 1 is a flowchart of the first embodiment of the soil-rock mixture characterization and concurrent solution method provided by the present invention;

[0048] Figure 2 shows the MPM-DEM characterization scheme of the soil-rock mixture provided by the present invention;

[0049] Figure 3 shows the contact force coupling scheme based on virtual particle contact detection technology provided by the present invention;

[0050] Figure 4 shows the gradation curve of the binary mixture provided by the present invention;

[0051] Figure 5 is a modeling diagram of a binary mixture sample with a particle size ratio of 1:10 provided by the present invention;

[0052] Figure 6 is a graph comparing the calculation results of the binary mixture sample with a particle size ratio of 1:10 provided by the present invention with those obtained by the discrete element method.

[0053] Figure 7 is a modeling diagram of a binary mixture sample with a particle size ratio of 1:5 provided by the present invention;

[0054] Figure 8 is a graph comparing the calculation results of the sample with a particle size ratio of 1:5 and the results with a particle size ratio of 1:10 provided by the present invention.

[0055] Figure 9 is a modeling diagram of the specimen for the indoor medium-sized triaxial test provided by the present invention;

[0056] Figure 10 is a graph comparing the simulation results of an indoor medium-sized triaxial test under a confining pressure of 200 kPa provided by the present invention with the actual results;

[0057] Figure 11 is a graph comparing the simulation results of an indoor medium-sized triaxial test under a confining pressure of 400 kPa provided by the present invention with the actual results;

[0058] Figure 12 is a graph comparing the simulation results of an indoor medium-sized triaxial test under 800 kPa confining pressure provided by the present invention with the actual results.

[0059] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0061] It should be noted that if the embodiments of the present invention involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of the components in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0062] Furthermore, if the embodiments of this invention involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the meaning of "and / or" throughout the text includes three parallel solutions; for example, "A and / or B" includes solution A, solution B, or a solution where both A and B are satisfied simultaneously. Furthermore, the technical solutions of the various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by this invention.

[0063] Soil-rock mixtures are a widespread material in geotechnical engineering. Their heterogeneity and nonlinearity result in highly complex mechanical properties. Accurately reproducing and predicting the mechanical properties of soil-rock mixtures is crucial for early warning of natural disasters such as debris flows, clastic flows, and landslides. Generally, there are two main methods for solving soil-rock mixture problems: continuous methods and discontinuous methods. In continuous methods, such as finite element method (FEM), finite difference method (FD), and smoothed particle hydrodynamics, the soil-rock mixture is treated as a continuous medium, and its mechanical behavior is described by a constitutive model. This can greatly simplify the problem and reduce the computational load, but the interaction between soil and rock cannot be considered, and some microscopic mechanisms cannot be reproduced. In discontinuous methods, such as discrete element method (DEM), each soil particle and rock can be modeled in great detail, which can account for microscopic mechanisms to some extent. However, for large-scale engineering problems, the computational load increases significantly. Furthermore, some microscopic parameters are difficult to calibrate, which is one of the important factors limiting the application of the DEM in practical engineering.

[0064] In view of this, the present invention proposes a method for characterizing and solving soil-rock mixtures concurrently. Figures 1 to 12 show specific embodiments of the method for characterizing soil-rock mixtures concurrently.

[0065] Please refer to Figure 1. The method for characterizing and solving the soil-rock mixture includes:

[0066] S30: After determining the soil and rocks of the target soil-rock mixture, the soil is modeled using the material point method and material points corresponding to the soil are generated; the rocks are modeled using the discrete element method and discrete elements corresponding to the rocks are generated.

[0067] S40: Describe the interaction between the soil and the rocks using a preset contact model, and construct an MPM-DEM characterization scheme for the target soil-rock mixture;

[0068] S50: Calculate the imaginary radius of each of the material points;

[0069] S60: Based on the virtual radius, determine the material points that are in contact with the discrete unit and mark them as virtual particles;

[0070] S70: The coupling contact force between the virtual particle and the discrete unit is calculated using the discrete element contact model;

[0071] S80: Based on the coupled contact force, construct the equilibrium equation and calculate the acceleration on each background grid node in the material point method, as well as the acceleration and angular acceleration on each discrete element;

[0072] S90: Calculate the real-time velocity and position parameters of the material point and the discrete unit based on the acceleration of the background grid node, the acceleration of the discrete unit, and the angular acceleration of the discrete unit.

[0073] In the technical solution of this invention, particles smaller than the boundary particle size are considered as soil, modeled using the Material Point Method (MPM), and their mechanical behavior is described using macroscopic constitutive models. The particle size of the material points can be appropriately increased to reduce the total number of particles and save computational resources. Particles larger than the boundary particle size are considered as rocks, modeled using the Discrete Element Method (DEM) to track the movement and contact behavior of individual rocks. Each material point is assigned a virtual radius for contact detection, and material points in contact with DEM elements are treated as virtual particles to calculate the coupling contact force. The coupling contact force is treated as an external boundary condition, applied to both the material points and the discrete elements, and their respective equilibrium equations are constructed. The equilibrium equations of the material points and the discrete elements are solved, and the stress and velocity of the material points, as well as their positions, are updated using the MUSL stress update scheme and the FLIP velocity update scheme, respectively. For the discrete elements, their velocity and position can be directly updated. It should be noted that the Material Point Method is a numerical method that combines Lagrange and Eulerian descriptions. Compared to mesh-based numerical methods, the material point method (MPM) overcomes problems such as mesh distortion; compared to other meshless methods, the MPM also boasts excellent computational efficiency. The discrete element method (DEM) has significant advantages in handling multi-body contact, making it highly suitable for dealing with rocks in soil-rock mixtures and the contact between soil and rocks. This invention couples the MPM and DEM, using the MPM to simulate soil and the DEM to simulate rocks, fully leveraging the advantages of both. Given the continuous nature of the MPM, macroscopic constitutive models can be used to describe the mechanical response of the soil, without needing to consider the contact behavior between individual soil particles. While maintaining the overall mechanical properties, the size of the material points can be appropriately increased to reduce the total number of particles in the model and improve computational efficiency; treating each rock as a separate DEM element without fine subdivision also reduces the total number of particles in the model and improves computational efficiency. Simultaneously, the contact between rocks and between rocks and soil is preserved, which can reflect the microscopic mechanisms to some extent. Finally, in the material point-discrete element coupling framework given in this invention, the calculations of the material point method and the discrete element method can be performed concurrently. Apart from calculating the coupling contact force, the other steps of the material point method and the discrete element method do not affect each other, which simplifies the code implementation difficulty and improves the calculation efficiency.

[0074] Specifically, the calculation of the imaginary radius S50 of each of the aforementioned material points includes:

[0075] S501: Using the formula Calculate the imaginary radius of each of the aforementioned material points, where, Let P be the imaginary radius of the P-th material point. Let k be the volume of the p-th material point. p Let be the porosity of the p-th material point.

[0076] In one embodiment, the discrete element method is used to calculate the coupling contact force S70 between the virtual particle and the discrete unit, including:

[0077] S701: Using the formula Calculate the resultant normal force, where, Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. p is the normal force exerted by the q-th material point on the P-th discrete unit, and pContact is the neighborhood list of particles in contact with the P-th discrete unit.

[0078] S702: Using the formula Calculate the tangential resultant force, where, Let be the resultant tangential force from the corresponding material point acting on the P-th discrete unit. p is the tangential force exerted by the q-th material point on the P-th discrete unit, and pContact is the neighborhood list of particles in contact with the P-th discrete unit.

[0079] S703: Combine the resultant normal force and the resultant tangential force to form the coupled contact force.

[0080] In this embodiment, the material point method and the discrete element method interact only when calculating the coupled contact force, thereby performing contact detection and contact force solution. Since their respective calculation processes are relatively independent, a concurrent solution method can be used, whereby the material point method and the discrete element method are calculated separately, and forced synchronization is only required when solving the coupled contact force.

[0081] Simultaneously, based on the coupled contact force, equilibrium equations are constructed, and the accelerations at each background grid node in the material point method, as well as the accelerations and angular accelerations S80 at each discrete element, are calculated, including:

[0082] S801: Constructing equilibrium equations using the material point method in, The resultant internal force at the i-th node in the material point method. The resultant external force at the i-th node in the material point method is... Let be the mass of the i-th node in the material point method. Let be the acceleration of the i-th node in the material point method;

[0083] S802: Solve the equilibrium equations of the material point method to obtain the accelerations at each background grid node in the material point method.

[0084] In this embodiment, the volume of the material point used for integration and the volume used for calculating the contact force are different. The volume used for integration is the total volume of the solid part and the porous part of the material point, while the volume used for calculating the contact force is the volume of the solid part of the material point.

[0085] Specifically, solving the equilibrium equations of the matter point method to obtain the acceleration S802 at each background grid node in the matter point method also includes:

[0086] S8021: Using the formula Calculate the resultant internal force, where N IP For shape functions, Physical force applied to a point of matter from the outside world. Let be the surface force acting on the i-th node in the material point method. Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. Let be the resultant tangential force from the corresponding material point acting on the Pth discrete unit;

[0087] S8022: Using formula Calculate the net external force, where, Let p be the stress at the p-th material point. Let B be the volume of the p-th material point. Ip The gradient of the shape function;

[0088] S8023: Calculate the acceleration at each background grid node based on the resultant internal force and the resultant external force.

[0089] Simultaneously, based on the coupled contact force, an equilibrium equation is constructed, and the acceleration on each background grid node in the material point method, as well as the acceleration and angular acceleration S80 on each discrete element, are calculated. This also includes:

[0090] S801′: Constructing the equilibrium equations using the discrete element method as well as in, Let be the resultant normal force from other discrete elements acting on the p-th discrete element. Let be the resultant tangential force from other discrete elements acting on the p-th discrete element. Let be the body force acting on the p-th discrete unit. Let be the rotational damping force acting on the p-th discrete element. Let p be the mass of the p-th discrete unit. Let I be the acceleration of the p-th discrete unit. p Let be the moment of inertia of the p-th discrete element. Let be the angular acceleration of the p-th discrete unit. Let n be the particle size of the p-th discrete unit. pq Let δ be the unit normal vector at the contact point. n,pq Let p be the normal stacking amount between the p discrete units and the q-th discrete unit (or material point). Let be the force exerted on the p-th discrete element by the q-th discrete element;

[0091] S802′: Solve the equilibrium equations of the material point method to obtain the acceleration and angular acceleration on each of the discrete units.

[0092] In this invention, the real-time parameters of the material point include the material point velocity, and the real-time parameters of the discrete unit include the discrete unit velocity.

[0093] Based on the acceleration of the background grid nodes, the acceleration of the discrete element, and the angular acceleration of the discrete element, the real-time parameter S90 of the material point and the discrete element is calculated, including:

[0094] S901: The velocity of the material points is updated in real time using the FLIP velocity update format, wherein the FLIP velocity update format is as follows: as well as Let be the neighborhood list of the adjacent nodes of the p-th material point, and let Δt be the time step of the material point. For the p-th matter point at t n+1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time, For the p-th matter point at t n +1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time;

[0095] S902: The velocity of the discrete unit is updated in real time using a preset velocity update format, wherein the preset velocity update format is as follows: as well as For the p-th matter point at t n+1 / 2 The speed of time, For the p-th matter point at t n-1 / 2 The speed of time, For the p-th matter point at t n+1 / 2 acceleration at any moment For the p-th matter point at t n-1 / 2 Acceleration at any moment.

[0096] Meanwhile, the real-time parameters of the material point include the position of the material point, and the real-time parameters of the discrete unit include the position of the discrete unit.

[0097] Based on the acceleration of the background grid nodes, the acceleration of the discrete element, and the angular acceleration of the discrete element, the real-time parameter S90 of the material point and the discrete element is calculated, including:

[0098] S901′: Adopted The update format updates the position of the material point in real time, wherein, For the p-th matter point at t n+1 The coordinates of the matter point at time [time]. For the p-th matter point at t n-1 The coordinates of the matter point at that moment;

[0099] S902′: The position of the discrete unit is updated in real time using a preset position update format, wherein the preset position update format is as follows: as well as For the p-th discrete unit at t n+1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n-1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n+1 Discrete cell coordinates at time 10:00 For the p-th discrete unit at t n-1 The coordinates of the discrete unit at time t.

[0100] Furthermore, based on the solution results, the stress, velocity, and position of the material points, as well as the velocity and position of the discrete elements S90, are updated according to a preset update format, and the process also includes:

[0101] S901″: Maps the momentum parameters of the material point onto the background mesh;

[0102] S902″: Based on the momentum parameters, solve for the velocities at the grid nodes using a preset equation, wherein the preset equation is: as well as

[0103] S903″: Update the strain at the material point using a strain update format, wherein the strain update format is set to...

[0104] S904″: Update the spinor at the material point using a spinor update format, wherein the spinor update format is set to...

[0105] S905″: Based on the strain and the spinor, update the stress at the material point using a preset constitutive model.

[0106] In this invention, before step S30, after determining the soil and rocks of the target soil-rock mixture, the soil is modeled using the material point method to obtain the material points corresponding to the soil, and the rocks are modeled using the discrete element method to obtain the discrete elements corresponding to the rocks, the invention further includes:

[0107] S10: Obtain the boundary particle size of the target soil-rock mixture;

[0108] S20: Determine the soil and rocks of the target soil-rock mixture based on the defined boundary particle size.

[0109] It should be noted that the verification examples provided in this invention mainly use an ideal elastoplastic constitutive model based on the Drucker-Prager yield criterion.

[0110] The following is a specific experiment to verify this:

[0111] Example 1:

[0112] The experimental sample was simplified to contain only two particle sizes—coarse and fine particles—in the same system; this is also known as a binary mixture. Its gradation curve is shown in Figure 4, with a particle size ratio of 1:10. The sample is shown in Figure 5; it is a cubic sample with a side length of 0.075 m.

[0113] First, undrained triaxial tests were conducted on pure fine particle samples under three confining pressures of 50 kPa, 100 kPa, and 150 kPa to calibrate the macroscopic parameters of the pure fine particles for simulation using the material point-discrete element method. The final simulation parameters are shown in Table 1.

[0114] Table 1 Macroscopic parameters of pure fine particles

[0115]

[0116] To simulate the strain hardening and softening phenomena in triaxial tests, a cubic hardening function and an exponential softening function are introduced, as shown in the following formulas:

[0117]

[0118] In the formula: ε p It is the cumulative equivalent plastic strain at the current moment; The equivalent plastic strain at which the peak value is reached is taken as 0.022 here; and These are the residual and peak internal friction angles. is the internal friction angle of the elastic segment, which is assumed to be equal to ; a, b, c, d are the coefficients of the hardening function, which are taken as 1075.64, -407.93, 16.39, and 0.32 respectively; H is the coefficient related to the softening rate, which is taken as 50.0.

[0119] Finally, the above parameters were applied to soil represented by the material point method, and triaxial tests were conducted under a confining pressure of 50 kPa with stone contents of 0%, 20%, and 40%. The stress-strain curves of each specimen are shown in Figure 6. It can be seen that under various stone content conditions, the results of the material point-discrete element method are very close to those of the discrete element method, accurately reflecting both hardening and softening phenomena. Moreover, the trend that the specimen reaches its peak strength earlier and at a higher peak strength with increasing stone content was also accurately reproduced.

[0120] Using the aforementioned macroscopic parameters, a single material point was used to replace multiple soil particles, and samples with a particle size ratio of 1:5 were established, containing three stone contents: 0%, 20%, and 40%, as shown in Figure 7. These samples were then subjected to triaxial tests under a confining pressure of 50 kPa, and the final results are shown in Figure 8.

[0121] It can be seen that even with a reduction in the number of material particles, the mechanical properties under a 1:10 particle size ratio can still be accurately reflected. Furthermore, by comparing the total number of particles in the samples, it was found that the number of particles in the 1:5 sample is only 1 / 8 of that in the 1:10 sample, indicating that at least 87.5% of the computational effort can be saved.

[0122] Example 2:

[0123] A medium-sized triaxial test was simulated using the material point-discrete element method. The samples were cuboids measuring 0.1m x 0.1m x 0.2m, with a stone content of 50%, but with three different particle size ratios. The model is shown in Figure 9. The gradation is shown in Table 2.

[0124] Table 2 Grading of Medium-Sized Triaxial Specimens

[0125] Particle size (m) gradation Grade 1 Grade 2 3 0.016–0.020 50% 0% 0% 0.010–0.016 0% 500% 0.005–0.010 0% 50% 0.5e-3–1.0e-3 0.6% 0.6% 0.6% 0.25e-3–0.5e-3 1.3% 1.3% 1.3% 0.075e-3–0.25e-3 14.1% 14.1% 14.1% <0.075e-3 34.0% 34.0% surface

[0126] For each gradation sample, simulations were performed using the material point-discrete element method under confining pressures of 200 kPa, 400 kPa, and 800 kPa. The density, Young's modulus, and Poisson's ratio of the stones were set to 2700 kg / m³, 50 GPa, and 0.2, respectively. Since the microscopic parameters of the soil particles were difficult to obtain directly, they were assumed to be consistent with those of the stones. The macroscopic parameters of the soil are shown in Table 3. To simulate the hardening phenomenon, a linear hardening function was used. The friction coefficients between the soil and stones, and between stones themselves, were both set to 0.3, and the bond between the soil and stones was ignored.

[0127]

[0128] The final results are shown in Figure 10-12. It can be seen that the simulated results are very close to the stress-strain curves of the actual triaxial test, and the continuous hardening phenomenon of the specimen is well reflected. At the same time, it can be seen that the particle size of the stones significantly affects the mechanical properties of the specimen. Especially under high confining pressure conditions, the smaller the particle size of the stones, the more pronounced the hardening phenomenon of the specimen, and the higher the peak strength of the specimen.

[0129] As demonstrated above, the material point method and the discrete element method (DEM) are coupled, fully leveraging their respective advantages. Using the material point method to discretize the soil, taking advantage of its continuous medium characteristics, eliminates the need to focus on the contact between soil particles, instead employing a macroscopic constitutive model to describe the soil's mechanical response. Simultaneously, without affecting the overall mechanical properties of the mixture, the size of the material points can be appropriately increased to reduce the number of particles in the model, saving computational resources. Using the DEM to model the stones accurately captures the motion and contact behavior of each stone. Furthermore, the interactions between stones and between stones and the soil can be considered, preserving, to some extent, the ability to reflect microscopic mechanisms. Using the material point method for soil modeling avoids the influence of mesh distortion, facilitating the simulation of large deformation problems. Moreover, the material point method boasts high computational efficiency among meshless methods, ensuring the efficiency of the coupled method. In this invention, the material point method and the DEM only interact when calculating the coupled contact force; their respective computational processes are independent, thus allowing for a concurrent solution method, which is convenient and significantly improves computational efficiency. Furthermore, when using this method to solve soil-rock mixture problems, compared to the discrete element method, the parameters of individual soil particles do not need to be precisely measured; only an approximate value is required, significantly reducing the modeling difficulty. In other words, the method provided by this invention is convenient, fast, and highly accurate, simplifying the workflow and facilitating the handling of soil-rock mixture problems with complex mechanical properties. It can better guide engineering practice and has broad application potential.

[0130] The above description is only a preferred embodiment of the present invention and does not limit the patent scope of the present invention. All equivalent structural transformations made under the concept of the present invention using the contents of the present invention specification and drawings, or direct / indirect applications in other related technical fields, are included within the patent protection scope of the present invention.

Claims

1. A method for characterizing and concurrently solving soil-rock mixtures, characterized in that, The method for characterizing and solving the soil-rock mixture includes: after determining the soil and rocks of the target soil-rock mixture, modeling the soil using the material point method and generating material points corresponding to the soil; modeling the rocks using the discrete element method and generating discrete elements corresponding to the rocks; describing the interaction between the soil and the rocks using a preset contact model and constructing an MPM-DEM characterization scheme for the target soil-rock mixture; calculating the virtual radius of each material point; determining the material points in contact with the discrete elements based on the virtual radius and marking them as virtual particles; calculating the coupling contact force between the virtual particles and the discrete elements using the discrete element contact model; constructing equilibrium equations based on the coupling contact force and calculating the acceleration on each background grid node in the material point method, as well as the acceleration and angular acceleration on each discrete element; calculating the real-time velocity and position parameters of the material points and the discrete elements based on the acceleration on the background grid nodes, the acceleration of the discrete elements, and the angular acceleration of the discrete elements.

2. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, Calculate the imaginary radius of each of the aforementioned material points, including: using the formula Calculate the imaginary radius of each of the aforementioned material points, where, Let P be the imaginary radius of the P-th material point. Let p be the volume of the p-th material point. Let be the porosity of the p-th material point.

3. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, The discrete element method is used to calculate the coupling contact force between the virtual particle and the discrete unit, including: using the formula... Calculate the resultant normal force, where, Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. Let be the normal force exerted by the q-th material point on the P-th discrete element. This is a list of the neighborhoods of particles in contact with the Pth discrete unit; the formula is used. Calculate the tangential resultant force, where, Let be the resultant tangential force from the corresponding material point acting on the P-th discrete unit. Let be the tangential force exerted by the q-th material point on the P-th discrete element. A neighborhood list of particles in contact with the Pth discrete unit; the resultant normal force and the resultant tangential force are combined to form the coupling contact force.

4. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, Based on the coupled contact force, an equilibrium equation is constructed, and the accelerations at each background grid node, as well as the accelerations and angular accelerations at each discrete element, are calculated in the material point method, including: constructing the material point method equilibrium equation. ,in, The resultant internal force at the i-th node in the material point method. The resultant external force at the i-th node in the material point method is... Let be the mass of the i-th node in the material point method. Let be the acceleration of the I-th node in the material point method; solve the equilibrium equation of the material point method to obtain the acceleration at each background grid node in the material point method.

5. The method for characterizing and concurrently solving soil-rock mixtures according to claim 4, characterized in that, Solving the equilibrium equations of the material point method to obtain the acceleration at each background grid node, as well as the acceleration and angular acceleration at each discrete element, also includes: using the formula Calculate the resultant internal force, where, For shape functions, Physical force applied to a point of matter from the outside world. Let be the surface force acting on the i-th node in the material point method. Let be the resultant normal force from the corresponding material point acting on the P-th discrete unit. Let P be the resultant tangential force from the corresponding material point acting on the Pth discrete element; using the formula Calculate the net external force, where, Let p be the stress at the p-th material point. Let p be the volume of the p-th material point. The gradient of the shape function is used; the acceleration at each background grid node is calculated based on the resultant internal force and the resultant external force.

6. The method for characterizing and concurrently solving soil-rock mixtures according to claim 4, characterized in that, Based on the coupled contact force, an equilibrium equation is constructed, and the accelerations at each background grid node in the material point method, as well as the accelerations and angular accelerations at each discrete element, are calculated. The method also includes: constructing the equilibrium equations for the discrete element method. as well as ,in, Let be the resultant normal force from other discrete elements acting on the p-th discrete element. Let be the resultant tangential force from other discrete elements acting on the p-th discrete element. Let be the body force acting on the p-th discrete unit. Let be the rotational damping force acting on the p-th discrete element. Let p be the mass of the p-th discrete unit. Let p be the acceleration of the p-th discrete unit. Let be the moment of inertia of the p-th discrete element. Let be the angular acceleration of the p-th discrete unit. Let p be the particle size of the p-th discrete unit. Let be the unit normal vector at the contact point. Let p be the normal stacking amount between the p discrete units and the q-th discrete unit or material point. Let p be the force exerted on the p-th discrete element by the q-th discrete element; solve the equilibrium equations of the discrete element method to obtain the acceleration and angular acceleration on each discrete element.

7. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, The real-time parameters of the material point include the material point velocity, and the real-time parameters of the discrete element include the discrete element velocity. The real-time parameters of the material point and the discrete element are calculated based on the acceleration of the background mesh nodes, the acceleration of the discrete element, and the angular acceleration of the discrete element. This includes updating the material point velocity in real-time using the FLIP velocity update format, wherein the FLIP velocity update format is... as well as , Let p be the list of neighboring nodes of the p-th material point. The time step of the material point. For the p-th matter point in The speed of time, For the p-th matter point in The speed of time, For the p-th matter point in The speed of time, For the p-th matter point in The velocity at any given time; the velocity of the discrete unit is updated in real time using a preset velocity update format, wherein the preset velocity update format is as follows: as well as , For the p-th matter point in The speed of time, For the p-th matter point in The speed of time, For the p-th matter point in acceleration at any moment For the p-th matter point in Acceleration at any moment.

8. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, The real-time parameters of the material point include the material point position, and the real-time parameters of the discrete unit include the discrete unit position; the real-time parameters of the material point and the discrete unit are calculated based on the acceleration on the background mesh nodes, the acceleration of the discrete unit, and the angular acceleration of the discrete unit, including: using... The update format updates the position of the material point in real time, wherein, For the p-th matter point in The coordinates of the matter point at time [time]. For the p-th matter point in The coordinates of the material point at time t; the position of the discrete unit is updated in real time using a preset position update format, wherein the preset position update format is as follows: as well as , For the p-th discrete unit in Discrete cell coordinates at time 10:00 For the p-th discrete unit in Discrete cell coordinates at time 10:00 For the p-th discrete unit in Discrete cell coordinates at time 10:00 For the p-th discrete unit in The coordinates of the discrete unit at time t.

9. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, Based on the solution results, the stress, velocity, and position of the material point, as well as the velocity and position of the discrete element, are updated according to a preset update format. The process also includes: mapping the momentum parameters of the material point onto the background mesh; and solving for the velocity at the mesh nodes using a preset equation based on the momentum parameters, wherein the preset equation is... as well as The strain at a material point is updated using a strain update format, wherein the strain update format is set to... Update the spinor at the material point using a spinor update format, wherein the spinor update format is set to... Based on the strain and the spinor, the stress at the material point is updated using a preset constitutive model.

10. The method for characterizing and concurrently solving soil-rock mixtures according to claim 1, characterized in that, After determining the soil and rocks of the target soil-rock mixture, before the steps of modeling the soil using the material point method to obtain the material points corresponding to the soil and modeling the rocks using the discrete element method to obtain the discrete elements corresponding to the rocks, the method further includes: obtaining the boundary particle size of the target soil-rock mixture; and determining the soil and rocks of the target soil-rock mixture based on the boundary particle size.