Discrete element method for multi-scale random mechanical properties of soil-rock mixture

By generating realistic rock block morphology through CT scanning and multi-sphere cluster algorithm, and combining Halton sequence and dual threshold correction, the error problem of soil-rock mixture mechanical simulation is solved, achieving efficient and accurate mechanical parameter statistics, and improving design accuracy and engineering application reliability.

CN122113542APending Publication Date: 2026-05-29THREE GORGES JINSHAJIANG CHUANYUN HYDROPOWER DEV CO LTD
View PDF 2 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
THREE GORGES JINSHAJIANG CHUANYUN HYDROPOWER DEV CO LTD
Filing Date
2026-02-04
Publication Date
2026-05-29

Smart Images

  • Figure CN122113542A_ABST
    Figure CN122113542A_ABST
Patent Text Reader

Abstract

The application discloses a discrete element analysis method for multi-scale random mechanical characteristics of soil and rock mixture, and comprises the following steps: obtaining real morphological characteristics of rock blocks; establishing a discrete element model based on the real morphological characteristics; generating a random sample space configuration; correcting particle contact relations; sequentially performing initial mechanical balance and triaxial shear loading on the corrected sample, obtaining stress-strain whole-process data, and extracting mechanical parameters; and based on the mechanical parameters, using Bootstrap resampling and taking a set variation coefficient as a convergence basis, obtaining a minimum simulation number, and constructing a quantitative relation table. The application solves the problems of morphological distortion, uneven distribution, unreasonable contact and arbitrary statistics existing in traditional analysis methods, improves simulation accuracy and engineering applicability, and can directly provide reliable parameters for slope, roadbed and dam body design.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology in geotechnical engineering, and in particular to a discrete element analysis method for multi-scale stochastic mechanical properties of soil-rock mixtures. Background Technology

[0002] In geotechnical engineering practice, soil-rock mixtures (such as landslide deposits, roadbed fill materials, and rockfill dam materials) are widely present. They are typical heterogeneous materials composed of soil matrix and randomly distributed rock blocks. Their mechanical properties directly determine the accuracy of analysis and design schemes for slope stability, roadbed settlement, and dam safety. However, soil-rock mixtures are characterized by a large range of rock block sizes, highly irregular shapes, and strong spatial randomness. In-situ testing is costly, produces highly discrete data, and makes it difficult to obtain repeatable mechanical parameters.

[0003] The Discrete Element Method (DEM) is currently the mainstream numerical method for analyzing the mechanical behavior of soil-rock mixtures in geotechnical engineering. Existing techniques typically use ideal spheres or simple ellipsoids to approximate rock particles and generate samples through random drop algorithms. However, ideal spheres / ellipsoids cannot reflect the edges, concavity, and aspect ratio of natural rock blocks, leading to significant deviations between simulation results and actual mechanical behavior. Furthermore, traditional random drop methods are prone to particle aggregation or voids, resulting in large initial contact force artifacts and insufficient sample representativeness. For example, a DEM simulation method for gravelly soil layers considering the randomness of particle shape disclosed in patent number CN201711063339.6 includes steps such as generating gradation curves, establishing a library of random convex polygon shapes, disk filling and covering, and generating polygonal particles. It focuses on the randomness of particle shape, but it is limited to convex polygonal particle shapes, while actual soil-rock mixtures contain a large number of non-convex stones. Moreover, the disk filling method has low computational efficiency, limiting its application in large-scale engineering projects. It does not consider the differences in cementation characteristics at the soil-rock contact surface and does not involve multi-physics coupling effects. The numerical simulation method for the failure law and mechanical characteristics of soil-rock mixtures disclosed in patent number CN202410570323.8 uses the meshless particle method and generates rock aggregates by employing random angle increments and random length sequences, belonging to the stochastic modeling method. However, its contact surface mechanical model is overly simplified, neglecting complex behaviors such as soil-rock embedding; it fails to address the issue of differential distribution of rock contact surfaces in the 3D model construction and lacks a micro-to-macro parameter transfer mechanism; it does not consider the coupling effects of multiple stress fields (such as the combination of rainfall-earthquake-low temperature); and the SPH method has high computational cost, limiting its engineering applications.

[0004] In terms of parameter statistics, the current practice is to select the sample size and number of simulations based on experience, lacking quantitative criteria for the convergence of mechanical parameters. Engineers find it difficult to determine the minimum reliable sample size and the minimum number of calculations, resulting in arbitrary parameter selection and overly conservative or dangerous design results. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a discrete element method for analyzing the multi-scale stochastic mechanical properties of soil-rock mixtures. By designing an analysis process involving CT scanning, multi-sphere clusters, Halton sequences, dual threshold correction, Bootstrap resampling, and quantification of relationships, this method simultaneously achieves the restoration of realistic rock block morphology, uniform spatial distribution of particles, quantitative determination of the statistical stability of mechanical parameters, and quantification of the relationship between engineering dimensions and computational costs within the discrete element framework. This solves the problems of morphological distortion, uneven distribution, and insufficient statistical reliability in existing technologies. It allows engineering projects such as slope stability and roadbed design to directly utilize the parameter packages output by the simulation without additional empirical correction, thus improving design accuracy and efficiency.

[0006] This invention provides a discrete element analysis method for the multi-scale stochastic mechanical properties of soil-rock mixtures, the specific technical solution of which is as follows: S1: Obtain the true morphological characteristics of the rock block through CT scanning technology; S2: Based on the real morphological features, a discrete element model considering the three-dimensional morphological features of particles is established using the multi-sphere cluster algorithm; S3: Generate random sample spatial configuration; S4: A dual-threshold contact correction mechanism is used to correct the particle contact relationship; S5: Initial mechanical equilibrium and triaxial shear loading are performed on the modified specimen in sequence to obtain stress-strain data for the entire process and extract mechanical parameters; S6: Based on the mechanical parameters, Bootstrap resampling is used and the set coefficient of variation is used as the convergence criterion to obtain the minimum number of simulations and construct a quantitative relationship table of sample size-stone content-minimum number of simulations.

[0007] Furthermore, the discrete element model is constructed as follows: Set a threshold for the soil-rock boundary; For soil matrix with particle size smaller than the soil-rock boundary threshold, spheres with a particle size between 1 mm and the soil-rock boundary threshold are used, and rolling resistance torque is set to compensate for morphological simplification errors for simulation. For rock blocks with a particle size not less than the soil-rock boundary threshold, a non-convex particle model corresponding to the real morphology is generated using a multi-sphere rigid cluster algorithm based on the actual morphological characteristics of the rock block, and a rock particle cluster library is established.

[0008] Furthermore, the generation of the random sample spatial configuration is specifically as follows: Based on the target stone content and gradation curve, a set of pure spheres is generated; Replace the occupiers with a particle size not less than the soil-rock boundary threshold with rock particle clusters randomly selected from the rock particle cluster library; The spatial coordinates of the centroid of rock grain clusters were determined using the Halton low-difference sequence algorithm.

[0009] Furthermore, when replacing the occupier spheres with rock particle clusters, a random rotation angle uniformly distributed in the range of [0, 2π] is applied to the clusters.

[0010] Furthermore, the modification of the particle contact relationship is as follows: Calculate the overlap ratio between soil particles and rock particles; Set deletion threshold With retention threshold Among them, the deletion threshold =1, retain threshold =0.01; The overlap rate of the deletion is not less than the deletion threshold. Soil particles, for retention threshold ≤overlap rate<deletion threshold The soil particles are proportionally reduced in radius; Generate a binary image of a cross-sectional CT slice and verify whether geometric penetration exists. If penetration exists, continue to reduce the radius proportionally until the binary image shows no geometric penetration, and proceed to the next step. Otherwise, proceed directly to the next step.

[0011] Furthermore, the initial mechanical equilibrium is as follows: A preset confining pressure is applied to the six boundary walls, and the wall pose loading servo pressure is dynamically adjusted through a force-displacement hybrid servo algorithm. The average contact force of the system is calculated at every first time step, and the real-time pressure of the boundary wall is obtained; When the average contact force fluctuation and the difference between the boundary wall pressure and the target value are both less than the set threshold, and the second time step is continuously set, the initial static equilibrium is reached.

[0012] Furthermore, the triaxial shear loading is specifically as follows: The lateral boundary walls are fixed, and the top wall is pressed down at a constant rate. Real-time acquisition of stress-strain data; calculation of deviatoric stress at each set axial strain step size; and simultaneous calculation of volumetric strain until the axial strain reaches the set threshold.

[0013] Furthermore, the mechanical parameters are extracted as follows: Based on the stress-strain process data, the first extreme point is extracted as the peak strength, the average stress within the set range of axial strain is taken as the residual strength, and the dilatation angle and initial elastic modulus are calculated based on the volumetric strain evolution.

[0014] Furthermore, the quantification relationship table is constructed as follows: Using Bootstrap, resampled data with replacement was used to generate several sets of resampled datasets. The sample size of each group is gradually increased from the initial value to the same as the original data, and the coefficient of variation of the key parameters is calculated. Using the set coefficient of variation threshold as the convergence condition, the minimum number of simulations was determined, and a quantitative relationship table of stone content, sample size, and minimum number of simulations was established.

[0015] The beneficial effects of this invention are as follows: 1. This invention establishes a discrete element model that considers the true morphological characteristics of rock blocks by using CT scanning and multi-spherical rigid clusters. Natural edges and concavities are restored in one step, and the morphological characteristics of rock blocks such as aspect ratio, flatness ratio and sphericity are accurately characterized. This significantly improves the accuracy of microstructure modeling and breaks through the limitations of traditional simplified spherical particles.

[0016] 2. This invention uses Halton low-difference sequences to optimize the spatial distribution of particles, and with the help of dual-threshold contact correction, generates random sample spatial configurations, achieving spatial uniform distribution and anisotropic arrangement of rock particles, thus ensuring the quality of sample generation.

[0017] 3. This invention introduces Bootstrap resampling technology to evaluate the statistical stability of mechanical parameters, determines the minimum number of simulations through coefficient of variation convergence analysis, establishes a quantitative relationship between sample size and the minimum number of simulations, solves the statistical reliability problem of traditional discrete element simulation, and reduces computational cost while ensuring the reliability of the results. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the method flow of the present invention.

[0019] Figure 2 It is a schematic diagram of the gradation curve of soil-rock mixture and irregular rock particle clusters.

[0020] Figure 3 This is a schematic diagram of the spatial configuration of a random sample.

[0021] Figure 4 This is a schematic diagram of particle correction.

[0022] Figure 5 This is a schematic diagram of CT numerical profile.

[0023] Figure 6 This is a schematic diagram for calculating mechanical parameters. Detailed Implementation

[0024] The technical solutions in the embodiments of the present invention are clearly and completely described in the following description. 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.

[0025] In the description of the embodiments of the present invention, it should be noted that the indicated orientation or positional relationship is based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship in which the product of the invention is conventionally placed during use, or the orientation or positional relationship in which those skilled in the art conventionally understand it during use. This is only for the convenience of describing the present invention and simplifying the description, and is not intended to indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, it should not be construed as a limitation of the present invention. Furthermore, the terms "first" and "second" are only used to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0026] In the description of the embodiments of the present invention, it should also be noted that, unless otherwise explicitly specified and limited, the terms "set" and "connection" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a direct connection or an indirect connection through an intermediate medium. Those skilled in the art can understand the specific meaning of the above terms in the present invention based on the specific circumstances.

[0027] Example 1 Embodiment 1 of the present invention discloses a discrete element analysis method for multi-scale stochastic mechanical properties of soil-rock mixtures, such as... Figure 1 As shown, the details are as follows: S1: A three-dimensional morphological database was established by acquiring the true morphological features of rock blocks through CT scanning; The actual morphological features include the aspect ratio (the ratio of the central axis of the rock block to the long axis) and the planarity ratio (the ratio of the short axis of the rock block to the central axis). Simultaneously calculate the sphericity index; sphericity The calculation formula is: ; in, Let be the volume of the rock particles. It represents the surface area of ​​the rock particles.

[0028] S2: Based on the aforementioned real morphological features, a discrete element model considering the three-dimensional morphological features of particles is established using the multi-sphere cluster algorithm. The specific process is as follows: A soil-rock boundary threshold is set. In this embodiment, a particle size of 5 mm is used as the soil-rock boundary threshold. For soil matrix with a particle size of less than 5 mm, 4-5 mm spheres are used for simulation, and a rolling resistance coefficient is set to compensate for the error caused by morphological simplification. In this embodiment, the rolling resistance coefficient is set to 0.2. For rock blocks with a particle size of 5-50 mm, based on a three-dimensional morphology database, a non-convex particle model corresponding to the actual morphology is generated using a multi-sphere rigid cluster algorithm, thus establishing a rock particle cluster library, such as... Figure 2 As shown in this embodiment, the cluster consists of 100 small balls with fixed relative positions.

[0029] In this embodiment, the rock particle model is stored in a six-column structured text format, consisting of cluster number, member sphere number, coordinates, and radius, and is based on... The value classifies rock particles into three categories: spheroids ( >0.9), sub-sphere (0.7< ≤0.9) and prismatic ( ≤0.7) S3: Generate random sample spatial configuration, the specific process is as follows: Based on the target stone content and gradation curve, the particle number distribution of each particle size group is determined. First, a set of pure spheres is generated according to the target stone content to ensure that the total number of particles meets the volume constraint condition. Subsequently, the occupier spheres with a particle size greater than 5 mm were replaced with rock particle clusters randomly selected from the rock particle cluster library. As a preferred embodiment, during the replacement process, a random rotation angle of [0, 2π] uniformly distributed was applied to the clusters to simulate the random spatial orientation of rock particles during the actual filling process. Finally, the spatial coordinates of the centroids of rock particle clusters were determined using the Halton low-difference sequence algorithm, realizing the random distribution of the spatial positions of rock particles. The constructed random soil-rock mixture sample is shown below. Figure 3 As shown, the red particles are soil, and the green particles are irregular rock fragments, which are evenly and randomly distributed in space.

[0030] S4: A dual-threshold contact correction mechanism is used to correct the particle contact relationship, as detailed below: Calculate the overlap rate between soil particles (simulated spheres of 4-5 mm) and rock particles (non-convex particle model); Set deletion threshold With retention threshold Among them, the deletion threshold =1, retain threshold =0.01; Delete blocks with a high overlap rate. Soil particles (ineffective soil particles completely embedded in the rock) at the retention threshold ≤overlap rate<deletion threshold The soil particles are proportionally reduced in radius to optimize the soil-rock contact state; By generating CT slices in the cross-section, binary images are used to verify whether geometric penetration exists. If penetration exists, the radius is continuously reduced proportionally until the binary image shows no geometric penetration, and the next step is executed. Otherwise, the next step is executed directly.

[0031] like Figure 4 As shown, the soil-rock mixture samples obtained in the end do not have any physical overlap on any cross section.

[0032] S5: Perform initial mechanical equilibrium and triaxial shear loading on the modified specimen in sequence to obtain stress-strain data (stress-strain response curve) and extract mechanical parameters; The initial mechanical equilibrium is as follows: A preset confining pressure (100-500kPa) is applied to the six boundary walls, and the wall pose loading servo pressure is dynamically adjusted by a force-displacement hybrid servo algorithm. In this embodiment, the iteration step size Δt=1×10-6s is set for dynamic compression balancing. The average contact force of the system is calculated every 1000 steps, and the real-time pressure of the boundary wall is obtained. When the difference between the average contact force fluctuation and the boundary wall pressure and the target value is less than 0.1% and is maintained for 500 time steps, the initial static equilibrium is determined to be reached.

[0033] The triaxial shear loading is specifically as follows: The lateral boundary wall is fixed, and the top wall is pressed down at a constant rate of 0.1% of the initial height of the specimen per second to achieve strain-controlled loading. Real-time acquisition of stress-strain data; calculation of deviatoric stress at 0.1% axial strain increments. Simultaneous calculation of volumetric strain (volume change) Divide by the total volume of the sample (until the axial strain reaches 40%).

[0034] like Figure 5 As shown, the first extreme point is extracted from the stress-strain response curve as the peak intensity. The average stress when the axial strain reaches 30%-40% is taken as the residual strength. Calculation of dilatation angle based on volumetric strain evolution , For shear strain, the elastic modulus under the reference confining pressure is determined by the slope of the tangent in the initial loading segment. .

[0035] S6: Based on the aforementioned mechanical parameters, Bootstrap resampling with replacement is used to generate 1000 sets of resampled datasets; The sample size of each group is [number]. Random simulation. The value increases from 3 to the same as the original data; For each set of resampled data, calculate the mean estimate of the key parameters and the corresponding variance estimate; Based on variance estimator and mean estimator Calculate the coefficient of variation ; Using a coefficient of variation < 2.5% as the convergence criterion, the minimum sample size was obtained. As the minimum number of simulations, a quantitative relationship table of sample size, stone content, and minimum number of simulations was constructed.

[0036] As shown in Table 1 below, Table 1 presents the minimum sample size obtained based on resampling analysis under different combinations of stone content and sample size.

[0037] Table 1: Relationship of minimum sample size obtained from resampling analysis under different combinations of stone content and sample size.

[0038] This invention is not limited to the specific embodiments described above. The invention extends to any new feature or combination disclosed in this specification, as well as any new method or process step or combination disclosed herein.

Claims

1. A discrete element analysis method for multi-scale stochastic mechanical properties of soil-rock mixtures, characterized in that, include: S1: Obtain the true morphological characteristics of the rock block; S2: Based on the aforementioned real morphological features, establish a discrete element model; S3: Generate random sample spatial configuration; S4: Correct the particle contact relationship; S5: Initial mechanical equilibrium and triaxial shear loading are performed on the modified specimen in sequence to obtain stress-strain data for the entire process and extract mechanical parameters; S6: Based on the mechanical parameters, Bootstrap resampling is used, and the set coefficient of variation is used as the convergence criterion to obtain the minimum number of simulations and construct a quantization relationship table.

2. The discrete element analysis method according to claim 1, characterized in that, The discrete element model is constructed as follows: Set a threshold for the soil-rock boundary; For soil matrix with particle size smaller than the soil-rock boundary threshold, spheres with a particle size ranging from 1 mm smaller than the soil-rock boundary threshold to the soil-rock boundary threshold were used, and rolling resistance torque was set for simulation. For rock blocks with a particle size not less than the soil-rock boundary threshold, a non-convex particle model corresponding to the actual shape of the rock block is generated based on the actual shape characteristics of the rock block, and a rock particle cluster library is established.

3. The discrete element analysis method according to claim 2, characterized in that, The specific details of generating the random sample spatial configuration are as follows: Based on the target stone content and gradation curve, a set of pure spheres is generated; Replace the occupiers with a particle size not less than the soil-rock boundary threshold with rock particle clusters randomly selected from the rock particle cluster library; The spatial coordinates of the centroid of rock grain clusters were determined using the Halton low-difference sequence algorithm.

4. The discrete element analysis method according to claim 3, characterized in that, When replacing the occupier spheres with rock particle clusters, apply a random rotation angle uniformly distributed in the range of [0, 2π] to the clusters.

5. The discrete element analysis method for multi-scale stochastic mechanical properties of soil-rock mixtures according to claim 1, characterized in that, The correction of the particle contact relationship is as follows: Calculate the overlap ratio between soil particles and rock particles; Set deletion threshold With retention threshold Among them, the deletion threshold =1, retain threshold =0.01; The overlap rate of the deleted text is not less than the deletion threshold. Soil particles, for retention threshold ≤overlap rate<deletion threshold The soil particles are proportionally reduced in radius; Generate a binary image of a cross-sectional CT slice and verify whether geometric penetration exists. If penetration exists, continue to reduce the radius proportionally until the binary image shows no geometric penetration, and proceed to the next step. Otherwise, proceed directly to the next step.

6. The discrete element analysis method according to claim 1, characterized in that, The initial mechanical equilibrium is as follows: A preset confining pressure is applied to the six boundary walls, and the wall pose loading servo pressure is dynamically adjusted through a force-displacement hybrid servo algorithm. The average contact force of the system is calculated at every first time step, and the real-time pressure of the boundary wall is obtained; When the average contact force fluctuation and the difference between the boundary wall pressure and the target value are both less than the set threshold, and the second time step is continuously set, the initial static equilibrium is reached.

7. The discrete element analysis method according to claim 1, characterized in that, The triaxial shear loading is specifically as follows: The lateral boundary walls are fixed, and the top wall is pressed down at a constant rate. Real-time acquisition of stress-strain data; calculation of deviatoric stress at each set axial strain step size; and simultaneous calculation of volumetric strain until the axial strain reaches the set threshold.

8. The discrete element analysis method according to claim 1, characterized in that, The mechanical parameters were extracted as follows: Based on the stress-strain process data, the first extreme point is extracted as the peak strength, the average stress within the set range of axial strain is taken as the residual strength, and the dilatation angle and initial elastic modulus are calculated based on the volumetric strain evolution.

9. The discrete element analysis method according to claim 1, characterized in that, The quantitative relationship table is constructed as follows: Several sets of resampled datasets are generated by resampling using Bootstrap; The sample size of each group is gradually increased from the initial value to the same as the original data, and the coefficient of variation of the key parameters is calculated. Using the set coefficient of variation threshold as the convergence condition, the minimum number of simulations was determined, and a quantitative relationship table of stone content, sample size, and minimum number of simulations was established.