Rockfill material space variability discrete element simulation method based on three-dimensional random field

Abstract

The invention discloses a three-dimensional random field-based rockfill material space variability discrete element simulation method. The method specifically comprises the following steps of: 1, constructing a triaxial test discrete element model and determining related parameters; step 2, generating a three-dimensional random field based on Karhunen-Loeve expansion; and step 3, random interpolation and assignment of particle contact parameters. Aiming at the technical defect that the real mechanical behavior of the rockfill material cannot be accurately reflected due to the fact that homogeneous hypothesis, two-dimensional modeling and continuous medium analysis are generally adopted in the prior art, the method comprises the following steps: constructing a three-dimensional random field of the rockfill material and introducing a discrete element particle system; more reasonable and more practical description of rockfill material space variability, macroscopic mechanical characteristics and mesoscopic mechanism deformation is realized.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of geotechnical engineering and computational mechanics, specifically involving a discrete element simulation method for the spatial variability of rockfill based on three-dimensional random fields. Background Technology

[0002] Rockfill is a heterogeneous granular material composed of particles of different sizes, shapes, and compositions, possessing unique engineering characteristics such as high strength, good permeability, and loose structure. Due to its superior engineering performance, rockfill is widely used in major infrastructure projects such as rockfill dams, roadbeds, high embankments, and slope reinforcement. In these engineering structures, the mechanical stability of the rockfill directly affects the safety and long-term service performance of the project; therefore, accurate modeling and simulation of its mechanical behavior are of great significance.

[0003] Rockfill exhibits significant spatial variability at the macroscopic scale, primarily due to the randomness and complexity of microscopic factors such as particle morphology, mechanical properties, and particle size distribution. Traditional modeling methods generally assume that rockfill is a spatially homogeneous material, neglecting the randomness and differences in its particle composition, mechanical properties, and spatial distribution. While this simplification reduces modeling complexity, it fails to accurately reflect the heterogeneity of local areas within the rockfill, leading to discrepancies between simulation results and actual engineering behavior. For example, the study in "Z. Qiu, W. Ling, L. Tian, ​​Effects of rock content and spatial distribution on the mechanical properties of soil–rock mixtures, ScientificReports 14(2024)80812" shows that the rock content and spatial distribution in soil-rock mixtures significantly affect their mechanical response, and ignoring material heterogeneity reduces the accuracy of simulation results.

[0004] Furthermore, many studies still employ two-dimensional numerical models to analyze rockfill structures. Two-dimensional models have inherent limitations in their geometric representation capabilities and description of structural deformation and stress differences along the thickness direction, making it difficult to realistically simulate the impact of the three-dimensional granular structure of rockfill on shear zone formation, seepage channel evolution, and local failure mechanisms. For example, the paper "D. Liu, L. Sun, J. Zhang, Mesoscopic analysis of the compaction characteristics of rockfill by 3D discrete element model, International Journal of Geomechanics 24(2023)04024154" points out that compared to two-dimensional simulations, three-dimensional models can more accurately characterize the evolution of rockfill granular structure and mechanical behavior, while two-dimensional models exhibit significant prediction biases.

[0005] In addition, most current studies still use the finite element method (FEM) to model and analyze riprap as a continuous medium. This method can describe the overall stress field at the macroscopic scale, but it is difficult to accurately characterize microscopic mechanical behaviors such as particle contact, force chain transmission, rolling and breakage. For example, "C. Zhang, X. Chen, Y. Liu, Discrete element methodsimulation of granular materials: particle breakage and shape effects, Granular Materials 2(2024)100035" points out that the finite element method (FEM) is limited in simulating particle breakage and shape effects compared to the discrete element method (DEM). Summary of the Invention

[0006] The purpose of this invention is to provide a discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field. This method addresses the shortcomings of existing technologies that commonly employ homogeneous assumptions, two-dimensional modeling, and continuous medium analysis, which fail to accurately reflect the true mechanical behavior of riprap. By constructing a three-dimensional random field for riprap and introducing a discrete element particle system, this invention achieves a more reasonable and practical description of the spatial variability, macroscopic mechanical characteristics, and microscopic deformation mechanisms of riprap.

[0007] The technical solution adopted in this invention is a discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field, specifically: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 3: Random interpolation and assignment of particle contact parameters.

[0008] The invention is further characterized in that: Step 1 is as follows: In the Discrete Element Method (PFC) software, the geometry, dimensions, and contact properties of the triaxial test riprap were defined, and the interparticle contact forces were calculated using the Discrete Element Method. Mechanical boundary conditions were applied at each contact point, ensuring the transfer of contact forces between particles, resulting in the triaxial test Discrete Element model. Then, the particle information of the riprap triaxial test was exported as a Ball-info.txt file using PFC's Python function area; simultaneously, the expected values ​​of the randomized normal and tangential stiffness required for the triaxial test were determined. with standard deviation .

[0009] Step 2 is as follows: Step 2.1: Constructing the random field matrix; Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; Step 2.3: Log-normal random field mapping.

[0010] Step 2.1 specifically involves: First, using Matlab, a uniform spatial resolution mesh is constructed within the cylindrical shape of the triaxial experimental discrete element model established in step 1. Then, four resolution mesh elements are extended to the outside of the mesh to form a cuboid random field mesh that completely encloses the triaxial experimental discrete element model in terms of length, width, and height. The expression for constructing the mesh is as follows: (1) (2) In the formula: L x It is the total length of the physical domain in the x-direction. L y It is the total length of the physical domain in the y-direction. L z It is the total length of the physical domain in the z-direction; n x It is the number of mesh cells in the x-direction. n y It represents the number of mesh cells in the y-direction. n z It represents the number of mesh cells in the z-direction; This is the total number of nodes; Next, the Matérn covariance function is used. To establish a spatial correlation model of random fields, and through the correlation length λ and variance σ²The Matérn covariance function is used to control the degree of randomness in spatial correlation models. The expression is as follows: (3) In the formula: d The distance between two nodes in space; λ is the smoothness parameter; λ is the correlation length, which is an indicator of the extent of spatial correlation influence. σ² The variance is calculated from the standard deviation. For Bessel functions; It is a gamma function.

[0011] Then based on the number of generated nodes and its coordinates and Matérn covariance function Calculate the covariance values ​​between nodes in the spatial correlation model. C ij Construct a covariance matrix C with a total dimension of N×N: (4) In the formula: C ij This represents the element in the i-th row and j-th column of the covariance matrix; It is the distance between the i-th and j-th nodes in the space.

[0012] Step 2.2 specifically involves: First, the covariance matrix C obtained by formula (4) is decomposed into eigenvalues ​​to obtain the eigenvalues. and eigenvectors By selecting the number of terms to retain based on energy, the generated random field is ensured to have sufficient accuracy, thus obtaining the Gaussian random field. The calculation formula is as follows: (5) In the formula: , Let i be the i-th eigenvalue and eigenvector; They are independent standard normal random variables; N KL Select the number of terms for energy; Then, the LHS method is used to generate the basic variable set. And through correlation coefficient Construct relevant random variables Specifically: (6) In the formula: The correlation coefficient is set by the user. Given a set of independent standard normal vectors, and They are independent of each other; , It conforms to the standard normal distribution.

[0013] Step 2.3 specifically involves: The Gaussian random field obtained in step 2.2 Based on the set random normal and tangential stiffness expectations 2.8×10 6 N / m, standard deviation A parameterized transformation is performed on the expected value in the range of 0.1-0.5, mapping it to a log-normal field. To match the distribution of material parameters, specifically: (7) (8) (9) Where: control parameters of the log-normal distribution and Based on target expected value and standard deviation The calculated log-normal distribution parameters, It is a Gaussian random field The average value, It is a Gaussian random field The standard deviation.

[0014] Step 3 specifically involves: Step 3.1: Use a suitable mesoscopic contact parameter interpolation mapping method; Step 3.2: Complete the random assignment of contact stiffness and the export of the random discrete element triaxial test model.

[0015] Step 3.1 specifically involves: Step 3.1.1: Loading and reconstructing the physical domain of the random field; In PFC software, the grid point data of the three-dimensional random field file random.txt is imported, and the one-dimensional field value vector is reconstructed into a three-dimensional tensor form to form a discrete representation of the continuous physical domain. Step 3.1.2: Particle positioning and interpolation; In the PFC software, the fish language is used to traverse all particles in the triaxial test of the rockfill. Based on the three-dimensional discrete mesh generated in step 3.1.1, the field values ​​of the mesh cell and its eight adjacent vertices are determined according to the coordinates of the particle's center. For particles completely contained within the mesh, the normal and tangential random parameter values ​​of the corresponding particles are calculated using the trilinear interpolation formula. For particles located at the boundary and only partially covering the mesh cell, trilinear interpolation combined with nearest neighbor extrapolation is used to avoid parameter loss and ensure numerical stability in the boundary region. Step 3.1.3: Denormalization mapping and structured file output; The random field values ​​obtained by interpolation in step 3.1.2 Since it is a dimensionless quantity, it needs to be mapped to the actual triaxial test contact stiffness range, i.e., 10, through a logarithmic scaling normalization strategy. 6 -10 8 The following is the mapping formula; (10) In the formula: For the i-th particle or particle contact pair, the actual contact stiffness is given. This represents the original value of the i-th particle; These represent the minimum and maximum values ​​of the field in the current simulation, respectively; ε is a small constant to prevent the denominator from being zero. Then, the particle ID, coordinates, radius, group, and corresponding random parameter values ​​are output as a structured file assigned-values.txt.

[0016] Step 3.2 specifically involves: After obtaining the structured file assigned-values.txt of the particles, the random normal and tangential stiffness contact parameters are read in the PFC software and mapped to each particle and particle contact in the discrete element model to form a loop. To ensure the continuity of parameter distribution and numerical stability, the contact stiffness is taken from the arithmetic mean of the stiffness of adjacent particles. After the assignment is completed, model samples containing spatial variability are exported.

[0017] The beneficial effects of this invention are: (1) The method of the present invention constructs an integrated modeling process from "three-dimensional random field to discrete element", which has the characteristics of portability, modularity and controllable parameters, and can be applied to different types of rockfill and various engineering scenarios. This method supports multi-parameter sensitivity analysis, reliability assessment, and automated batch calculation of examples, significantly improving the model's engineering adaptability and application efficiency. Compared to the homogeneous discrete element model, the stress-strain curve obtained by this invention exhibits mechanical properties that better match experimental observations: peak intensity is reduced by approximately 15%, peak strain is delayed by approximately 2.0%–3.5%, and the post-peak softening process is smoother. Furthermore, the curve dispersion increases with increasing confining pressure, consistent with the random fluctuations in the mechanical properties of riprap as external constraints change in actual engineering. This effectively compensates for the shortcomings of homogeneous models, which have overly concentrated responses and are difficult to reflect experimental dispersion. This invention can capture, at the microscopic level, the characteristics of more dispersed force chain structures, more random fracture locations, and a slower rate of contact force anisotropy enhancement in riprap. It quantitatively reveals that spatial variability weakens the stability of the dominant force chain path, leading to increased randomness in the macroscopic mechanical response. From a mechanistic perspective, it establishes a causal chain of "parameter spatial variability—microscopic response dispersion—macroscopic mechanical dispersion," providing theoretical support for explaining the dispersion of the engineering properties of riprap.

[0018] (2) The method of this invention addresses the shortcomings of existing technologies that commonly employ homogeneous assumptions, two-dimensional modeling, and continuous medium analysis, which fail to accurately reflect the true mechanical behavior of riprap. By constructing a multi-scale spatial random field for riprap and introducing a three-dimensional discrete element particle system, this invention achieves a more reasonable and practical description of the spatial variability, macroscopic mechanical characteristics, and microscopic deformation mechanisms of riprap. This method can characterize particle contact, force chain evolution, and breakage behavior at the microscopic level, and reflect the stress-strain curve of riprap at the macroscopic scale to demonstrate the influence of spatial variability on the macroscopic and microscopic mechanical responses of riprap. This provides a more reliable and applicable new modeling method for the analysis of riprap material properties, evaluation of engineering performance, and research on structural safety. Attached Figure Description

[0019] Figure 1 This is a flowchart of the discrete element simulation method for spatial variability of rockfill based on KL expansion, as described in this invention. Figure 2 This is a diagram showing the number of randomly expanded items in an embodiment of the present invention; Figure 3 This is a schematic diagram of the trilinear interpolation method used to completely cover the particles in the embodiment of the present invention; Figure 4 This is a schematic diagram of the mesh portion covering particles in this embodiment of the invention using trilinear interpolation combined with boundary proximity extrapolation. Figure 5 This is a schematic diagram of the boundary particles using trilinear interpolation combined with the boundary nearest extrapolation method in an embodiment of the present invention; Figure 6 This is a simulation result of a random discrete element method based on a triaxial test of riprap in an embodiment of the present invention; Figure 7 This is a graph showing the random and deterministic stress-strain curves under a confining pressure of 400 kPa in an embodiment of the present invention. Figure 8 This is a graph showing the random and deterministic stress-strain curves under a confining pressure of 800 kPa in an embodiment of the present invention. Figure 9 This is a graph showing the random and deterministic stress-strain curves under a confining pressure of 1200 kPa in an embodiment of the present invention. Figure 10 This is a graph showing the random and deterministic stress-strain curves under a confining pressure of 1600 kPa in an embodiment of the present invention. Figure 11 This is a statistical diagram of normal and tangential fractures of rockfill material determined and randomized under a confining pressure of 400 kPa in an embodiment of the present invention. Figure 12 This is a statistical diagram of normal and tangential fractures of rockfill material determined and randomized under a confining pressure of 1600 kPa in an embodiment of the present invention. Figure 13 This is a two-dimensional polar coordinate normal contact force distribution diagram in an embodiment of the present invention; Figure 14 This is a two-dimensional polar coordinate tangential contact force distribution diagram in an embodiment of the present invention; Figure 15 This is a graph showing the variation of the average normal contact force of the riprap material with strain in an embodiment of the present invention; Figure 16 This is a graph showing the variation of the normal Fourier fitting coefficients of the riprap material with strain in an embodiment of the present invention. Figure 17 This is a graph showing the variation of average tangential contact force of riprap with strain in an embodiment of the present invention; Figure 18 This is a graph showing the variation of the tangential Fourier fitting coefficients of the riprap material with strain in an embodiment of the present invention. Detailed Implementation

[0020] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0021] This invention provides a discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields, such as... Figure 1 As shown, the specific steps include: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 1 is as follows: In the Discrete Element Method (PFC) software, the geometry, dimensions, and contact properties of the triaxial test riprap were defined, and the interparticle contact forces were calculated using the Discrete Element Method. Mechanical boundary conditions were applied at each contact point, ensuring the transfer of contact forces between particles, resulting in the triaxial test Discrete Element model. Then, the particle information of the riprap triaxial test was exported as a Ball-info.txt file using PFC's Python function area; simultaneously, the expected values ​​of the randomized normal and tangential stiffness required for the triaxial test were determined. with standard deviation .

[0022] Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; A three-dimensional random field is generated, and eigenvalue decomposition is performed using Karhunen-Loeve expansion. This step introduces the spatial variability of the riprap into the simulation. Step 2 is implemented as follows: Step 2.1: Constructing the random field matrix; To implement the above method, firstly, a uniform spatial resolution mesh is constructed using Matlab within the cylindrical size of the triaxial experimental discrete element model established in step 1. Then, four resolution mesh elements are extended to the outside of the mesh to form a cuboid random field mesh whose length, width, and height can completely enclose the triaxial experimental discrete element model. The expression for constructing the mesh is as follows: (1) (2) In the formula: L x It is the total length of the physical domain in the x-direction. L y It is the total length of the physical domain in the y-direction. L z It is the total length of the physical domain in the z-direction; n x It is the number of mesh cells in the x-direction. n y It represents the number of mesh cells in the y-direction. n z It represents the number of mesh cells in the z-direction; This is the total number of nodes; Next, the Matérn covariance function is used. To establish a spatial correlation model of random fields, and through the correlation length λ and variance σ² The Matérn covariance function is used to control the degree of randomness in spatial correlation models. The expression is as follows: (3) In the formula: dThe distance between two nodes in space; λ is the smoothness parameter; λ is the correlation length, which is an indicator of the extent of spatial correlation influence. σ² The variance is calculated from the standard deviation. For Bessel functions; It is a gamma function; Then based on the number of generated nodes and its coordinates and Matérn covariance function Calculate the covariance values ​​between nodes in the spatial correlation model. C ij Construct a covariance matrix C with a total dimension of N×N: (4) In the formula: C ij This represents the element in the i-th row and j-th column of the covariance matrix; It is the distance between the i-th and j-th nodes in the space; Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; To improve computational efficiency, the Karhunen-Loeve expansion is used to perform eigenvalue decomposition on the covariance matrix, approximating the random field as a finite-mode form. Specifically, the covariance matrix C obtained from formula (4) is decomposed into eigenvalues ​​to obtain the eigenvalues. and eigenvectors By selecting the number of terms to retain based on energy, the generated random field is ensured to have sufficient accuracy, thus obtaining the Gaussian random field. The calculation formula is as follows: (5) In the formula: , Let i be the i-th eigenvalue and eigenvector; They are independent standard normal random variables; N KL Select the number of terms for energy; Based on the above steps, the principal component coefficients of the Gaussian random field are determined by the standard normal variables. Therefore, Latin hypercube sampling (LHS) is used to generate independent Gaussian samples and control the correlation coefficients. Specifically, the LHS method is used to generate the basic variable set. And can be obtained through correlation coefficient Construct relevant random variables Specifically: (6) In the formula: The correlation coefficient is set by the user. Given a set of independent standard normal vectors, and They are independent of each other; , It conforms to a standard normal distribution; Step 2.3: Mapping to a log-normal random field; Subsequently, the Gaussian random field obtained in step 2.2 Based on the set random normal and tangential stiffness expectations 2.8×10 6 N / m, standard deviation A parameterized transformation is performed on the expected value in the range of 0.1-0.5, mapping it to a log-normal field. To match the distribution of material parameters, specifically: (7) (8) (9) Where: control parameters of the log-normal distribution and Based on target expected value and standard deviation The calculated log-normal distribution parameters, It is a Gaussian random field The average value, It is a Gaussian random field The standard deviation of the standard deviation ensures that the generated random field of material parameters meets the expected distribution and fluctuation characteristics. This yields the three-dimensional random field file random.txt.

[0023] Step 3: Random interpolation and assignment of particle contact parameters; Based on the generated three-dimensional random field, the particle contact parameters of the riprap in the triaxial test are randomly assigned values. This step transforms external randomness into spatial variability within the discrete element model of the riprap. Step 3 is implemented as follows: Step 3.1: Use a suitable mesoscopic contact parameter interpolation mapping method; To construct an efficient and reliable mapping method that reasonably assigns the physical domain parameters of continuous random fields to the particle system, this study uses a trilinear interpolation method to achieve interpolation mapping.

[0024] Step 3.1.1: Loading and reconstructing the physical domain of the random field; In PFC software, the grid point data of the obtained three-dimensional random field file random.txt is imported, and the one-dimensional field value vector is reconstructed into a three-dimensional tensor form to form a discrete representation of the continuous physical domain.

[0025] Step 3.1.2: Particle positioning and interpolation; In the PFC software, the FILE language was used to traverse all particles in the triaxial test of the rockfill. Based on the three-dimensional discrete mesh generated in step 3.1.1, the field values ​​of the mesh cell and its eight adjacent vertices were determined according to the particle's sphere center coordinates. For particles completely contained within the mesh, the normal and tangential random parameter values ​​were calculated using trilinear interpolation. For particles located at the boundary and only partially covering the mesh cell, a combination of trilinear interpolation and nearest neighbor extrapolation was used to avoid parameter loss and ensure numerical stability in the boundary region.

[0026] Step 3.1.3: Denormalization mapping and structured file output; The random field values ​​obtained by interpolation in step 3.1.2 Since it is a dimensionless quantity, it needs to be mapped to the actual triaxial test contact stiffness range, i.e., 10, through a logarithmic scaling normalization strategy. 6 -10 8 The following is the mapping formula.

[0027] (10) In the formula: This represents the actual contact stiffness corresponding to the i-th particle (or particle contact pair); This represents the original value of the i-th particle; These represent the minimum and maximum values ​​of the field in the current simulation, respectively; ε is a small constant to prevent the denominator from being zero. Then, the particle ID, coordinates, radius, group, and corresponding random parameter values ​​are output as a structured file assigned-values.txt.

[0028] Step 3.2: Complete the random assignment of contact stiffness and the export of the random discrete element triaxial test model; After obtaining the structured file `assigned-values.txt` containing the randomized values ​​of the particles, the randomized normal and tangential stiffness contact parameters are read from the PFC software and mapped to each particle and particle contact in the discrete element model to form a loop. To ensure the continuity and numerical stability of the parameter distribution, the contact stiffness is taken as the arithmetic mean of the stiffness of adjacent particles. After the assignment is completed, model samples containing spatial variability are exported, laying the foundation for subsequent mechanical response analysis under multiple random field conditions. Simultaneously, the above process is encapsulated into a modular program to improve the computational efficiency and methodological versatility of discrete element random field modeling.

[0029] Example 1 In this embodiment, the expected values ​​of the random normal and tangential stiffness required for the triaxial test are determined according to the article "Shao Lei, Chi Shichun, Zhang Yong. Study on triaxial shear test of rockfill based on particle flow [J]. Rock and Soil Mechanics, 2013, 34 (3): 711-720." 2.8×106 N / m, standard deviation Set the value to 0.1-0.5 times the expected value. Establish a triaxial experimental discrete element model with a diameter of 300mm and a height of 650mm, setting the spatial resolution to 0.03. The determined cuboid random field mesh is 720mm in diameter and 1420mm in height. Then, following step 2 above, use the Matérn covariance function with a correlation length λ of 0.1 to form a covariance matrix C of total dimension N×N. Perform eigenvalue decomposition based on the Karhunen-Loeve expansion on the covariance matrix C. In this example, the number of expansion terms is selected to be the first 500 terms to achieve an energy of 98%. Figure 2 As shown. Then, the Latin hypercube sampling (LHS) method is used to generate independent Gaussian samples and the correlation coefficient is set. The value is 0.5. Then, following the steps, a log-normal random field mapping is performed to obtain the three-dimensional random field file random.txt. Afterwards, following step 3, random interpolation and assignment of particle contact parameters are performed to obtain a spatial variability model of the random discrete element method for the triaxial test of the riprap, as shown below. Figure 6 As shown. In step 3.1.2: particle positioning and interpolation... Figure 3 The mesh completely covers the particles using trilinear interpolation, such as... Figure 3 As shown, the mesh portion covering the particles uses trilinear interpolation combined with boundary proximity extrapolation, as shown in the figure. Figure 4 As shown, the boundary particles are calculated using trilinear interpolation combined with boundary proximity extrapolation. Figure 5 As shown.

[0030] (1) Analysis of macroscopic mechanical properties of riprap considering spatial variability The aforementioned simulation method can effectively demonstrate the spatial randomness and variability of riprap. Based on this, this study will analyze the changes in the macroscopic mechanics of riprap considering spatial variability to demonstrate the necessity of simulating the spatial variability of riprap. The stress-strain curves from triaxial tests, as characteristic curves characterizing the entire process of material mechanics response, can reflect the deformation, strength, and failure characteristics of riprap. The results of the stress-strain curves from triaxial tests of riprap under different confining pressures considering spatial variability are shown below. Figures 7 to 10 The confining pressures shown are 400 kPa, 800 kPa, 1200 kPa, and 1600 kPa. (The rest of the text appears to be a fragment and doesn't translate accurately.) Figures 7 to 10It can be seen that: ① The calculation results of the 20-times random discrete element model have certain dispersion characteristics within a small range, and are concentrated around its mean. Considering the spatial variability of the riprap, the stress-strain curve of the riprap is significantly discrete, which proves from the side that the spatial variability discrete element simulation method of the riprap can effectively show the randomness of the mechanical properties of the riprap; ② Compared with the deterministic discrete element model, the peak strength of the sample calculated by the random field is reduced by about 15%, and the peak strain is delayed by about 2-3.5% overall. The delay increases with the increase of confining pressure, and the stress-strain curve shows a more obvious "hyperbolic" feature and a smoother post-peak softening behavior, which is closer to the results of indoor tests; ③ With the increase of confining pressure, the dispersion of the calculation results of the 20-times random discrete element model increases significantly: under low confining pressure, the width of the gray shading band is about 0.39 MPa, while under high confining pressure, it expands to 0.9 MPa, and the dispersion range doubles. This increased dispersion indicates a more significant difference in stress paths under high confining pressure, reflecting the amplified effect of the randomness of the internal structure of the riprap on its mechanical characteristics as confining pressure increases. The above analysis shows that constructing a random field in the discrete element model can more accurately reflect the mechanical characteristics of the riprap. Therefore, using a random field in discrete element simulation is essential, and the discrete element simulation method for the spatial variability of riprap proposed in this study can effectively achieve this goal.

[0031] (2) Analysis of the microscopic deformation mechanism of riprap considering spatial variability While the stress-strain curves of triaxial tests can reflect the macroscopic differences between deterministic and stochastic discrete element models (DEMs), they are insufficient to explain the microscopic differences. The following analysis examines the fracture characteristics and contact stress distribution of the stochastic DEM at the microscopic level to demonstrate the influence mechanism of material spatial variability on triaxial tests of riprap. In the deterministic model, strong chains and breakage are continuous and concentrated, while in the stochastic model they are more dispersed. Due to the influence of factors such as shape, size, and strength in actual engineering, the distribution of strong chains and breakage is often uneven; therefore, the stochastic model is more realistic and exhibits better performance. Fracture occurs when the normal or tangential force exceeds the bond strength during triaxial loading. The fracture process is analyzed by statistically analyzing the number of bond failures. The statistical analysis of different types of bond failures in riprap under lower confining pressures of 400 kPa and higher confining pressures of 1600 kPa is shown below. Figures 11 to 12 As shown in the figure, the results indicate that: ① Under the same confining pressure, the number of failures in the normal, tangential, and total directions all increase steadily with axial strain, with tangential failures significantly exceeding normal failures; ② Increasing the confining pressure significantly promotes failure development; ③ The failure rate of the deterministic model is generally higher than that of the stochastic model, especially in tangential failures. This suggests that deterministic structures are prone to forming continuous force chains, leading to local overloads, while stochastic distributions weaken the concentration effect, making the force transmission system more rational. Figure 13 and Figure 14The figures show the distribution of normal and tangential contact forces in two-dimensional polar coordinates. As can be seen, the normal force contact force distribution gradually evolves from an ellipse to a "peanut shape," while the tangential force contact force distribution shows increased concentration from the peak to the residual stage. Increased confining pressure significantly amplifies the numerical value and distribution range of the contact forces, while spatial variability makes the contact force increase more slowly and the distribution narrower. Quantitative indicators further show that randomness slows down the rapid increase in contact force anisotropy, delaying local instability. The changes in the average normal contact force and normal Fourier fitting coefficient of the riprap are shown below. Figure 15 and Figure 16 As shown, the variations in the average tangential contact force and tangential Fourier fitting coefficients of the riprap are as follows: Figure 17 and Figure 18 As shown. By Figure 15 and Figure 18 It can be observed that: ① The average normal contact force and the average tangential contact force increase with increasing axial strain and then stabilize or slightly decrease in the later stages of strain. The curves for low confining pressure (400 kPa) and high confining pressure (1600 kPa) are similar in shape. The stochastic model shows a slower increase and a lag in the decrease during the residual stage compared to the deterministic model. ② The Fourier fitting coefficients represent the degree of anisotropy of the average contact force. The normal fitting coefficients show an initial increase followed by a decrease with increasing strain. The anisotropy is significantly stronger at high confining pressure than at low confining pressure. The stochastic model shows a smoother anisotropy change compared to the deterministic model, reflecting the moderating and stabilizing effect of stochasticity on the evolution of the force chain direction. The tangential fitting coefficients fluctuate with strain, and the influence of confining pressure and stochasticity is not significant.

[0032] In summary, the method of this invention addresses the shortcomings of existing technologies that commonly employ homogeneous assumptions, two-dimensional modeling, and continuous medium analysis, which fail to accurately reflect the true mechanical behavior of riprap. By constructing a three-dimensional random field for riprap and introducing a discrete element particle system, it achieves a more reasonable and practical description of the spatial variability, macroscopic mechanical characteristics, and microscopic deformation mechanisms of riprap.

[0033] Example 2 The discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields is as follows: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 3: Random interpolation and assignment of particle contact parameters.

[0034] Example 3 The discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields is as follows: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 1 is as follows: In the Discrete Element Method (PFC) software, the geometry, dimensions, and contact properties of the triaxial test riprap were defined, and the interparticle contact forces were calculated using the Discrete Element Method. Mechanical boundary conditions were applied at each contact point, ensuring the transfer of contact forces between particles, resulting in the triaxial test Discrete Element model. Then, the particle information of the riprap triaxial test was exported as a Ball-info.txt file using PFC's Python function area; simultaneously, the expected values ​​of the randomized normal and tangential stiffness required for the triaxial test were determined. with standard deviation .

[0035] Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 3: Random interpolation and assignment of particle contact parameters.

[0036] Example 4 The discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields is as follows: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 2 is as follows: Step 2.1: Constructing the random field matrix; Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; Step 2.3: Log-normal random field mapping.

[0037] Step 3: Random interpolation and assignment of particle contact parameters.

[0038] Example 5 The discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields is as follows: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 2 is as follows: Step 2.1: Constructing the random field matrix; Step 2.1 specifically involves: First, using Matlab, a uniform spatial resolution mesh is constructed within the cylindrical shape of the triaxial experimental discrete element model established in step 1. Then, four resolution mesh elements are extended to the outside of the mesh to form a cuboid random field mesh that completely encloses the triaxial experimental discrete element model in terms of length, width, and height. The expression for constructing the mesh is as follows: (1) (2) In the formula: L x It is the total length of the physical domain in the x-direction. L y It is the total length of the physical domain in the y-direction. L z It is the total length of the physical domain in the z-direction; n x It is the number of mesh cells in the x-direction. n y It represents the number of mesh cells in the y-direction. n z It represents the number of mesh cells in the z-direction; This is the total number of nodes; Next, the Matérn covariance function is used. To establish a spatial correlation model of random fields, and through the correlation length λ and variance σ² The Matérn covariance function is used to control the degree of randomness in spatial correlation models. The expression is as follows: (3) In the formula: d The distance between two nodes in space; λ is the smoothness parameter; λ is the correlation length, which is an indicator of the extent of spatial correlation influence. σ² The variance is calculated from the standard deviation. For Bessel functions; It is a gamma function; Then based on the number of generated nodes and its coordinates and Matérn covariance function Calculate the covariance values ​​between nodes in the spatial correlation model. C ij Construct a covariance matrix C with a total dimension of N×N: (4) In the formula: C ij This represents the element in the i-th row and j-th column of the covariance matrix; It is the distance between the i-th and j-th nodes in the space.

[0039] Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; Step 2.3: Log-normal random field mapping.

[0040] Step 3: Random interpolation and assignment of particle contact parameters.

[0041] Example 6 The discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields is as follows: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 2 is as follows: Step 2.1: Constructing the random field matrix; Step 2.1 specifically involves: First, using Matlab, a uniform spatial resolution mesh is constructed within the cylindrical shape of the triaxial experimental discrete element model established in step 1. Then, four resolution mesh elements are extended to the outside of the mesh to form a cuboid random field mesh that completely encloses the triaxial experimental discrete element model in terms of length, width, and height. The expression for constructing the mesh is as follows: (1) (2) In the formula: L x It is the total length of the physical domain in the x-direction. L y It is the total length of the physical domain in the y-direction. L z It is the total length of the physical domain in the z-direction; n x It is the number of mesh cells in the x-direction. n y It represents the number of mesh cells in the y-direction. n z It represents the number of mesh cells in the z-direction; This is the total number of nodes; Next, the Matérn covariance function is used. To establish a spatial correlation model of random fields, and through the correlation length λ and variance σ² The Matérn covariance function is used to control the degree of randomness in spatial correlation models. The expression is as follows: (3) In the formula: d The distance between two nodes in space; λ is the smoothness parameter; λ is the correlation length, which is an indicator of the extent of spatial correlation influence. σ² The variance is calculated from the standard deviation. For Bessel functions; It is a gamma function; Then based on the number of generated nodes and its coordinates and Matérn covariance function Calculate the covariance values ​​between nodes in the spatial correlation model. C ij Construct a covariance matrix C with a total dimension of N×N: (4) In the formula: C ij This represents the element in the i-th row and j-th column of the covariance matrix; It is the distance between the i-th and j-th nodes in the space.

[0042] Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; Step 2.2 specifically involves: First, the covariance matrix C obtained by formula (4) is decomposed into eigenvalues ​​to obtain the eigenvalues. and eigenvectors By selecting the number of terms to retain based on energy, the generated random field is ensured to have sufficient accuracy, thus obtaining the Gaussian random field. The calculation formula is as follows: (5) In the formula: , Let i be the i-th eigenvalue and eigenvector; They are independent standard normal random variables; N KL Select the number of terms for energy; Then, the LHS method is used to generate the basic variable set. And through correlation coefficient Construct relevant random variables Specifically: (6) In the formula: The correlation coefficient is set by the user. Given a set of independent standard normal vectors, and They are independent of each other; , It conforms to the standard normal distribution.

[0043] Step 2.3: Log-normal random field mapping.

[0044] Step 3: Random interpolation and assignment of particle contact parameters.

Claims

1. A discrete element method for simulating the spatial variability of riprap based on three-dimensional random fields, characterized in that, Specifically: Step 1: Construct a discrete element model for the triaxial experiment and determine the relevant parameters; Step 2: Generation of a 3D random field based on Karhunen-Loeve expansion; Step 3: Random interpolation and assignment of particle contact parameters.

2. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 1, characterized in that, Step 1 is as follows: In the Discrete Element Method (PFC) software, the geometry, dimensions, and contact properties of the triaxial test riprap were defined, and the interparticle contact forces were calculated using the Discrete Element Method. Mechanical boundary conditions were applied at each contact point, ensuring the transfer of contact forces between particles, resulting in the triaxial test Discrete Element model. Then, the particle information of the riprap triaxial test was exported as a Ball-info.txt file using PFC's Python function area. Simultaneously, the expected values ​​of the randomized normal and tangential stiffness required for the triaxial test were determined. with standard deviation .

3. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 1, characterized in that, Step 2 is as follows: Step 2.1: Constructing the random field matrix; Step 2.2: Eigenvalue decomposition based on Karhunen-Loeve expansion; Step 2.3: Log-normal random field mapping.

4. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 3, characterized in that, Step 2.1 specifically involves: First, using Matlab, a uniform spatial resolution mesh is constructed within the cylindrical shape of the triaxial experimental discrete element model established in step 1. Then, four resolution mesh elements are extended to the outside of the mesh to form a cuboid random field mesh that completely encloses the triaxial experimental discrete element model in terms of length, width, and height. The expression for constructing the mesh is as follows: (1) (2) In the formula: L x It is the total length of the physical domain in the x-direction. L y It is the total length of the physical domain in the y-direction. L z It is the total length of the physical domain in the z-direction; n x It represents the number of mesh cells in the x-direction. n y It represents the number of mesh cells in the y-direction. n z It represents the number of mesh cells in the z-direction; This is the total number of nodes; Next, the Matérn covariance function is used. To establish a spatial correlation model of random fields, and through the correlation length λ and variance σ² The Matérn covariance function is used to control the degree of randomness in spatial correlation models. The expression is as follows: (3) In the formula: d The distance between two nodes in space; λ is the smoothness parameter; λ is the correlation length, which is an indicator of the extent of spatial correlation influence. σ² The variance is calculated from the standard deviation. For Bessel functions; It is a gamma function; Then based on the number of generated nodes and its coordinates and Matérn covariance function Calculate the covariance values ​​between nodes in the spatial correlation model. C ij Construct a covariance matrix C with a total dimension of N×N: (4) In the formula: C ij This represents the element in the i-th row and j-th column of the covariance matrix; It is the distance between the i-th and j-th nodes in the space.

5. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 4, characterized in that, Step 2.2 specifically involves: First, the covariance matrix C obtained by formula (4) is decomposed into eigenvalues ​​to obtain the eigenvalues. and eigenvectors By selecting the number of terms to retain based on energy, the generated random field is ensured to have sufficient accuracy, thus obtaining the Gaussian random field. The calculation formula is as follows: (5) In the formula: , Let i be the i-th eigenvalue and eigenvector; They are independent standard normal random variables; N KL Select the number of terms for energy; Then, the LHS method is used to generate the basic variable set. And through correlation coefficient Construct relevant random variables Specifically: (6) In the formula: The correlation coefficient is set by the user. Given a set of independent standard normal vectors, and They are independent of each other; , It conforms to the standard normal distribution.

6. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 5, characterized in that, Step 2.3 specifically involves: The Gaussian random field obtained in step 2.2 Based on the set random normal and tangential stiffness expectations 2.8×10 6 N / m, standard deviation A parameterized transformation is performed on the expected value in the range of 0.1-0.5, mapping it to a log-normal field. To match the distribution of material parameters, specifically: (7) (8) (9) Where: control parameters of the log-normal distribution and Based on target expected value and standard deviation The calculated log-normal distribution parameters, It is a Gaussian random field The average value, It is a Gaussian random field The standard deviation.

7. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 1, characterized in that, Step 3 specifically involves: Step 3.1: Use a suitable mesoscopic contact parameter interpolation mapping method; Step 3.2: Complete the random assignment of contact stiffness and the export of the random discrete element triaxial test model.

8. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 7, characterized in that, Step 3.1 specifically involves: Step 3.1.1: Loading and reconstructing the physical domain of the random field; In PFC software, the grid point data of the three-dimensional random field file random.txt is imported, and the one-dimensional field value vector is reconstructed into a three-dimensional tensor form to form a discrete representation of the continuous physical domain. Step 3.1.2: Particle positioning and interpolation; In the PFC software, the fish language is used to traverse all particles in the triaxial test of the rockfill. Based on the three-dimensional discrete mesh generated in step 3.1.1, the field values ​​of the mesh cell and its eight adjacent vertices are determined according to the coordinates of the particle's center. For particles completely contained within the mesh, the normal and tangential random parameter values ​​of the corresponding particles are calculated using the trilinear interpolation formula. For particles located at the boundary and only partially covering the mesh cell, trilinear interpolation combined with nearest neighbor extrapolation is used to avoid parameter loss and ensure numerical stability in the boundary region. Step 3.1.3: Denormalization mapping and structured file output; The random field values ​​obtained by interpolation in step 3.1.2 Since it is a dimensionless quantity, it needs to be mapped to the actual triaxial test contact stiffness range, i.e., 10, through a logarithmic scaling normalization strategy. 6 -10 8 The following is the mapping formula; (10) In the formula: For the i-th particle or particle contact pair, the actual contact stiffness is given. This represents the original value of the i-th particle; These represent the minimum and maximum values ​​of the field in the current simulation, respectively; ε is a small constant to prevent the denominator from being zero. Then, the particle ID, coordinates, radius, group, and corresponding random parameter values ​​are output as a structured file assigned-values.txt.

9. The discrete element method for simulating the spatial variability of riprap based on a three-dimensional random field according to claim 8, characterized in that, Step 3.2 specifically involves: After obtaining the structured file assigned-values.txt of the particles, the random normal and tangential stiffness contact parameters are read in the PFC software and mapped to each particle and particle contact in the discrete element model to form a loop. To ensure the continuity of parameter distribution and numerical stability, the contact stiffness is taken from the arithmetic mean of the stiffness of adjacent particles. After the assignment is completed, model samples containing spatial variability are exported.

Images