Rock slope discrete element modeling method based on random field theory

Abstract

The invention discloses a rock slope discrete element modeling method based on a random field theory, and the method comprises the steps: generating an initial homogeneous particle sample according to physical parameters and geometric boundaries, and applying target confining pressure through servo control, so that the sample reaches a static equilibrium state; performing state resetting on the sample, endowing particle contact with a bonding characteristic, and constructing a homogeneous rock model; extracting a unique identifier and a space coordinate of each contact, constructing a random field representing spatial variability of cohesion and tensile strength based on the unique identifier and the space coordinate, and distributing a corresponding strength parameter for each contact; and applying a gravitational field amplified in proportion to the model to simulate prototype stress, executing graded excavation simulation, and evaluating the slope stability by monitoring system displacement data. According to the method, high-precision mapping of random field parameters to the discrete element model is realized, and strict consistency of contact parameters in the whole simulation process under complex working conditions is ensured.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation technology for slope engineering, and in particular relates to a discrete element modeling method for rock slopes based on random field theory. Background Technology

[0002] Rock slope engineering is a key component of geotechnical engineering, widely used in major projects such as transportation, water conservancy, mining, and infrastructure construction. As engineering projects become larger and more complex, rock slope stability has become a critical factor affecting project safety and economic benefits. To deepen our understanding of rock slope instability mechanisms and improve the accuracy of predictions, numerical simulation technology has been widely applied in this field.

[0003] The Discrete Element Method (DEM) has been widely used in the stability analysis of rock slope engineering due to its ability to effectively simulate discontinuous deformation processes such as rock mass fracture propagation and block separation. However, current DEM modeling of rock slopes still commonly employs homogeneous or simplified heterogeneous parameter assignment methods, such as zonal assignment or layered parameter allocation based on geological survey results. These modeling methods fail to accurately reflect the continuous spatial variability and correlation of rock mass parameters (such as cohesion c and internal friction angle φ), leading to significant discrepancies between simulation results and actual engineering performance in terms of structural stress transmission and failure path evolution.

[0004] In recent years, random field theory has provided a rigorous theoretical foundation and effective tools for describing the spatial variability of geological medium parameters. In particular, spatial random field modeling methods based on the Karhunen-Loève (KL) expansion method can efficiently construct physical parameter distribution fields with predetermined spatial correlation structures, given the known statistical characteristics of the parameters. However, current research on the application of this method in DEM simulation is relatively limited. Existing attempts typically use the spatial coordinates of the cells as the index for parameter matching. This approach carries the risk of inaccurate parameter mapping due to errors in particle (or cell) coordinate calculations. Furthermore, it lacks a mechanism for tracking and matching parameters after particle displacement under complex conditions, making it difficult to ensure strict consistency between parameters and the model throughout the simulation process.

[0005] Therefore, there is an urgent need to develop a modeling method that integrates random field theory and discrete element method (DEM), which can achieve high-precision parameter assignment through a contact parameter mapping mechanism based on unique identifiers, and ensure strict consistency of contact parameters throughout the simulation process under complex working conditions. Summary of the Invention

[0006] The purpose of this invention is to provide a discrete element modeling method for rock slopes based on random field theory, which solves the problems in the existing technology that cannot truly reflect the spatial variability of rock mass parameters, inaccurate mapping due to coordinate matching, and lack of parameter tracking mechanism after displacement. It achieves high-precision mapping of random field parameters to discrete element models and ensures strict consistency of contact parameters throughout the simulation process under complex working conditions.

[0007] The technical solution adopted in this invention is a discrete element modeling method for rock slopes based on random field theory, comprising the following steps:

[0008] Step S1: Generate an initial homogeneous particle sample based on the set physical parameters and geometric boundaries;

[0009] Step S2: Apply the target confining pressure to the initial homogeneous particle sample, and use servo control to make the sample reach static equilibrium under the target confining pressure;

[0010] Step S3: Reset the state of the specimen after stress initialization and assign bonding properties to the contact between particles to construct a homogeneous rock model;

[0011] Step S4: Generate and export identification information containing a unique identifier and spatial location coordinates for each contact in the homogeneous rock model;

[0012] Step S5: Based on the identification information, construct a random field characterizing the spatial variability of cohesion and tensile strength parameters, and assign corresponding cohesion and tensile strength parameters to each contact;

[0013] Step S6: Apply a scaled-up gravity field to the model with the corresponding parameters to simulate the prototype stress, and perform a staged excavation simulation until the model reaches a new equilibrium state after excavation, and obtain displacement monitoring data for stability analysis.

[0014] Further, step S1 includes:

[0015] Based on the set physical parameters and geometric boundaries, a randomly distributed set of particles is generated within the boundary; a contact model and parameters are specified for the contact between particles in the set, and static equilibrium is achieved through iterative cycles and kinetic energy reduction; particles located outside the boundary are deleted to form an initial homogeneous particle sample.

[0016] Further, step S2 includes:

[0017] The average normal stress on the model boundary is used as the monitoring variable, and the moving speed of the boundary wall is used as the control variable. Based on the difference between the target confining pressure and the monitoring variable, combined with the servo gain parameter, the control variable is adjusted to make the monitoring variable approach the target confining pressure.

[0018] The servo gain parameter is determined based on the model's total normal contact stiffness in the pressure direction, the geometric dimensions of the boundary wall, and the calculation time step.

[0019] Furthermore, the criteria for determining the static equilibrium state include: the relative error between the monitored variable and the target confining pressure is less than a preset threshold, and the unbalanced force ratio of the particle system within the model is less than a preset threshold or the control variable approaches zero.

[0020] Further, step S3 includes: zeroing the displacement and rotation angle of all particles in the sample; and defining a combined contact model and its parameters for all contacts, including linear and parallel bonded portions.

[0021] Further, step S4 includes: traversing all contacts in the model, assigning a unique identifier to each contact and recording its spatial coordinates; and writing the identifier and spatial coordinates into a data file.

[0022] Further, step S5 includes:

[0023] Based on unique identifiers and spatial coordinates, autocovariance kernels and cross-covariance kernels characterizing the spatial correlation between rock mass cohesion and tensile strength are constructed to build a covariance matrix.

[0024] Based on monitoring data, the parameter mean vector and covariance matrix are updated by constructing a conditional random field. If there is no monitoring data, the unconditional mean and covariance matrix are used directly.

[0025] Further, in step S5, the covariance matrix is ​​decomposed into eigenvalues, the number of truncated modes is determined based on the cumulative contribution rate of eigenvalues, and a random field sample vector containing the cohesion and tensile strength values ​​of each contact node is generated using KL expansion.

[0026] Further, in step S5, the random field sample vector is split into a cohesion parameter field and a tensile strength parameter field, and organized with the unique identifier of the contact and exported as a parameter file;

[0027] The parameter file is read, the corresponding parameter value is found according to the unique identifier of the contact, and the bonding strength parameter attribute of the corresponding contact in the discrete element model is assigned respectively.

[0028] Further, step S6 includes:

[0029] The increased gravitational acceleration is set according to the geometric scaling factor of the model and prototype sizes, so that the model reaches initial static equilibrium under the corresponding gravitational field.

[0030] Virtual monitoring points are set up at key locations on the slope to record the displacement and depth data of particles at each monitoring point;

[0031] By defining the boundary of the excavation area and deleting the particles within it, iterative calculations are performed. The equilibrium convergence criterion is that the unbalanced force ratio of the model system is less than a preset threshold or the displacement rate of the monitoring point approaches zero. This results in the acquisition of the new equilibrium state and displacement monitoring data.

[0032] The beneficial effects of this invention are:

[0033] 1. This invention proposes a discrete element method (DEM) modeling method based on two-dimensional binary random field theory. It constructs spatially correlated random fields by integrating key mechanical parameters such as rock mass cohesion and tensile strength at the contact layer. Compared to existing DEMs that often employ homogeneous parameters, partitioned assignments, or simple perturbations, which struggle to reflect the spatial covariance structure and correlation of parameters, this method constructs autocovariance and cross-covariance matrices using a covariance kernel. It then combines Karhunen–Loève expansion and Latin hypercube sampling to generate random field samples. If borehole or monitoring data is available, a conditional random field can be constructed based on the Gaussian condition formula, ensuring that parameters are consistent with observed values ​​at measurement points while maintaining a reasonable spatial variation structure in unmeasured areas. This method achieves, for the first time, integrated random field modeling of "multi-parameter correlation + conditional update" in discrete element methods, thereby more meticulously characterizing weak interlayers, locally weakened zones, and multi-scale heterogeneous structures in rock slopes, significantly improving the accuracy of numerical simulation in reproducing real geological conditions.

[0034] 2. This invention innovatively establishes an independent and stable unique identifier (ID) for each contact, addressing the shortcomings of existing random field-DEM coupling methods that rely on real-time spatial coordinates or grid position assignments. Specifically, during large deformations and contact reconstruction, parameters experience mapping errors and "drift" due to coordinate changes or contact ID updates. The ID and its corresponding spatial coordinates are exported to an external program to generate a random field, and then written back to the discrete element model in the form of ID-parameter values. This ID-based lookup and assignment mechanism eliminates the dependence on real-time coordinates, establishing a stable index framework that remains unchanged with displacement and contact reconstruction, thus achieving precise and persistent binding between random field parameters and contact elements. Even under complex conditions such as self-weight loading, excavation unloading, and earthquakes, where particles undergo significant displacement or stress redistribution, each contact maintains a one-to-one correspondence with its initial random field parameters, fundamentally overcoming the problems of parameter inconsistency and drift in existing technologies.

[0035] 3. Application results show that by introducing a random field at the contact level, this invention enables the microscopic parameters such as slope cohesion and tensile strength to exhibit a continuously correlated random distribution in space. Under different random field implementations, the failure initiation location, crack propagation path, and final slip surface morphology all exhibit randomness and diversity. It can naturally simulate complex instability modes such as preferential failure of local weak zones and bending and bifurcation of slip zones. Compared with the homogeneous discrete element model, which only forms a single regular slip surface, it is closer to the actual progressive failure and multi-path evolution process of rock slopes. Attached Figure Description

[0036] 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 these drawings without creative effort.

[0037] Figure 1 This is a flowchart of the discrete element modeling of rock slopes based on random field theory in an embodiment of the present invention.

[0038] Figure 2 This is a comparison diagram of the discrete element modeling process and effect of rock slope in an embodiment of the present invention. (a) is the initial homogeneous particle sample, (b) is the sample and contact information extraction after preloading equilibrium, (c) is the contact bond strength distribution with random field assignment, (d) is the slope geometry after random field assignment, and (e) is the distribution of homogeneous model parameters.

[0039] Figure 3(a) is a diagram of slope failure morphology in the random field model of an embodiment of the present invention.

[0040] Figure 3(b) is a diagram of the slope failure morphology of the homogeneous model in the embodiment of the present invention.

[0041] Figure 3(c) is a distribution diagram of the displacement modulus field of the random field model in an embodiment of the present invention.

[0042] Figure 3(d) is a displacement modulus field distribution diagram of the homogeneous model in the embodiment of the present invention. Detailed Implementation

[0043] 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 some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0044] like Figure 1 As shown in the example, a discrete element method for modeling rock slopes based on random field theory includes the following steps:

[0045] Step S1: Rock Sample Establishment: This step is used to generate an initial homogeneous particle sample that meets the size and physical property requirements. Specifically, it includes parameter setting, particle distribution, model specification, and sample optimization operations.

[0046] The basic physical parameters of the sample, such as geometric dimensions, particle density, porosity, damping coefficient, and particle radius range, are set. The computational domain boundary is defined and a rectangular wall boundary is established, while the random seed number is set. Within the established rectangular boundary, particles are randomly distributed according to the preset porosity and particle radius range by calling the ball distribute command to generate an initial homogeneous particle sample. A linear contact model is uniformly assigned to the sample, and the contact modulus, stiffness ratio, and damping coefficient are assigned. The initial unbalanced force of the system is released through iterative cycles and kinetic energy reduction to achieve static equilibrium. Excess particles located outside the boundary are deleted to obtain a homogeneous rock sample with standardized dimensions and stable structure.

[0047] Step S2: Rock Stress Relief and Pre-compression Treatment: This step applies the target confining pressure to the established homogeneous sample, eliminating initial unbalanced stress through a servo control mechanism, and stabilizing the sample under a predetermined stress state. The specific implementation process is as follows:

[0048] S21: Servo Control Mechanism and Algorithm

[0049] Define the average normal stress on the model boundary. To monitor the variable, this stress is obtained by dividing the total normal contact force acting on the model boundary by the boundary area; where the total normal contact force originates from the interaction between the boundary wall (or the dedicated pressure-applying particle group) and the particles inside the specimen; the moving velocity of the model boundary wall is defined. For control variables; in a proportional control law, each calculation step According to the target confining pressure Compared with current monitored stress Stress error between Adjust the movement speed of the boundary wall; the adjustment formula is:

[0050]

[0051] in, express The movement speed of the boundary wall at any given time. Indicates the previous moment The movement speed of the boundary wall express The average normal stress on the model boundary monitored at all times, where G is the servo gain parameter.

[0052] To ensure the stability of the servo control system and achieve fast, oscillation-free convergence, the servo gain parameter... It is not set arbitrarily, but is determined by calculation based on the overall stiffness characteristics of the model; the specific calculation method is as follows: First, the total normal stiffness of the boundary wall (e.g., left / right side wall) in contact with the sample particles in the direction of the applied confining pressure is calculated. This total stiffness reflects the overall deformation resistance of the system; combined with the length of the wall action... and calculation time step The servo gain is determined using the following formula:

[0053]

[0054]

[0055] in, , These represent the servo gains in the X and Y directions, respectively. , These represent the total normal contact stiffness in the X and Y directions, respectively; , These represent the geometric dimensions of the corresponding wall in the X and Y directions, respectively; The adjustment coefficient is used to fine-tune the damping characteristics of the control system to ensure that the servo process reaches a near-critical damping state, thereby eliminating stress oscillations and achieving high efficiency and stability. Its value range is generally 0.1~0.3.

[0056] S22: Stress Balance and Stability Criterion: The servo control process continues until the system reaches the specified static equilibrium state; the equilibrium criterion is set to simultaneously satisfy two conditions: stress convergence and kinetic energy stability.

[0057] Stress convergence: Monitoring stress in boundary walls Containment pressure of the target The relative error is less than the preset small threshold. ,Right now:

[0058]

[0059] Kinetic stability: The average unbalanced force ratio of the particle system within the model is less than a preset threshold of 1. 10 -5 or the movement speed of the boundary wall Approaching zero.

[0060] When the above conditions are met, it is considered that the specimen has successfully completed stress initialization under the target confining pressure, entered a stable state, and the wall position is kept fixed or the stress is locked by a servo mechanism.

[0061] Step S3: Homogeneous Model Construction: This step is used to reset the pre-compression treated sample to a stable initial state and construct a contact constitutive model with cohesive properties to simulate the mechanical behavior of the intact rock mass. The specific implementation process is as follows:

[0062] The displacement and rotation angles of all particles (balls) in the sample are zeroed, and the entire system is transitioned from the initial state after pre-compression to a stable state that can be used to assign material model values ​​through kinetic energy reduction operation. The parameters required for the parallel bond model are defined, including elastic modulus, stiffness ratio, tensile strength, bond force, bond angle and bond gap. At the same time, linear contact frames are specified for the contact between particles (ball-ball) and between particles and the wall (ball-facet), and their deformation calculation methods and resistance parameters are set. The mean contact model is initialized, and a "linear and parallel bond" combination model is uniformly specified for all contacts. The corresponding parameters of the linear part and the parallel bond part are assigned in sequence, and the bond gap is set.

[0063] Step S4: Contact Information Identification and Export: This step is used to create and export unique identification information and spatial location data for all contacts in the model, establishing an index foundation for subsequent high-precision mapping of random field parameters. The specific implementation process is as follows:

[0064] Create a target text file contact_info.txt in write mode to establish a stable data output channel and traverse all contacts within the current discrete element model. Assign a unique identifier ID to each contact, which is independent of the contact.id that may be refreshed during calculations within the discrete element software. Store this ID stably in a user-defined additional attribute field reserved for that contact, such as the contact.extra(,1) field. Simultaneously, obtain the spatial coordinates (x, y) of each contact. Concatenate the unique identifier ID, x-coordinate, and y-coordinate into a single line of string data according to a predetermined order and format, and write it line by line into the opened text file to ensure orderly data recording. After traversal, close the file to ensure that all data has been physically written. Use a Python script to read the generated contact_info.txt file and convert it into a structured CSV file for efficient and accurate reading and retrieval by subsequent random field analysis programs.

[0065] Step S5: Random Field Parameter Assignment Based on Identification Information: This step is used to achieve accurate mapping of external random field parameters to contact parameters of the discrete element model. By establishing a spatially correlated random field and assigning parameters based on unique identifiers (IDs), high-fidelity representation of parameter spatial variability is ensured. The specific implementation process is as follows:

[0066] S51: Constructing a bivariate covariance matrix using a covariance kernel:

[0067] Read the unique identifiers (IDs) of all contacts and their corresponding spatial coordinates (x, y) from the contact_info.csv file, and treat these coordinate points as nodes of a two-dimensional random field; use a Gaussian correlation function (i.e., a squared exponential kernel) to characterize the spatial correlation structure of the parameters.

[0068] Define rock mass cohesion autocovariance kernel Its expression is:

[0069]

[0070]

[0071]

[0072] in, , Let i and j represent the spatial coordinates of the i-th and j-th contact nodes, respectively. This represents the spatial distance between two points in the grid. Represents an exponential function; for The variance; yes Standard deviation; yes The mean; yes coefficient of variation; It is the relevant length of c.

[0073] Similarly, the tensile strength of rock mass is defined. autocovariance kernel Cross-covariance kernel between cohesion c and tensile strength t :

[0074]

[0075]

[0076] in, It is the standard deviation of the tensile strength t. Its variance; express and The correlation coefficient is generally taken as 0.2 to 0.6; It is the cross-correlation length.

[0077] Construct the covariance matrix :

[0078]

[0079] in, for The transpose of the matrix indicates that the cross covariance is symmetric; express It is a real matrix with 2n rows and 2n columns.

[0080] S52: Constructing a Conditional Random Field Based on Monitoring Data: This step aims to address how to use existing engineering monitoring data to constrain a random field, ensuring that the generated results are consistent with measured values ​​at known points, thereby improving the model's representation accuracy. If no monitoring data is available, an unconditional random field is used. The specific implementation process is as follows:

[0081] If monitoring data obtained through drilling, field testing, or other methods exists, and the cohesion (c) and tensile strength (t) values ​​at a specific spatial location (measuring point) can be derived, then proceed to this step to generate a conditional random field; otherwise, directly use the unconditional mean. and unconditional covariance Proceed to the next steps.

[0082] The monitoring data from m measuring points are used to construct a 2m-dimensional random vector of measuring points. Its definition is:

[0083]

[0084] in, Let represent the spatial coordinates of the i-th measuring point, and denote its specific actual observed value (i.e., known data) as . ; Indicates the location at coordinates The cohesive force at that location is a random variable; Indicates the location at coordinates The tensile strength at a given location is a random variable.

[0085] According to Gaussian process theory, with known measurement point data... Under the condition of [condition], the random field vector formed by all unknown nodes (i.e., all contact points) It follows a conditional Gaussian distribution:

[0086]

[0087] Among them, conditional mean Conditional covariance It is obtained from the following formula:

[0088]

[0089]

[0090] in, Let be a 2n-dimensional joint random field vector consisting of the cohesion and tensile strength of all n contact nodes (unknown points); for The corresponding actual observed value; for The unconditional mean vector; for The unconditional mean vector; for The unconditional covariance matrix; for The unconditional covariance matrix; for and The cross-covariance block matrix between them for The transpose of .

[0091] To simplify the subsequent expressions, we define the final parameters of the conditional random field as the mean vector. Covariance matrix for:

[0092]

[0093] If there is no monitoring data, then: .

[0094] S53: Random field sample generation based on KL expansion:

[0095] The final covariance matrix (i.e., conditional covariance) obtained in step S52 or unconditional covariance Perform eigenvalue decomposition to solve for its eigenvalues ​​and eigenvectors:

[0096]

[0097] in, For the i-th eigenvalue, This is the corresponding feature vector.

[0098] To balance computational efficiency and accuracy, a truncation criterion based on the cumulative contribution rate of eigenvalues ​​is used to determine the number of dominant modes to be retained. That is, before calculation The ratio of the sum of the largest eigenvalues ​​to the sum of all eigenvalues; when this ratio reaches a preset energy threshold... (For example, when taking 90%~95%), determine The value:

[0099]

[0100] This criterion ensures that the truncated feature space can encompass most of the fluctuation energy and variability information of the original random field.

[0101] Construct the truncated eigenvector matrix and the square root of eigenvalues :

[0102]

[0103]

[0104] in, To truncate the mode matrix, the eigenvectors corresponding to the first r largest eigenvalues ​​are... A matrix formed by arranging columns; To truncate the diagonal matrix of eigenvalue square roots, the diagonal elements are the square roots of the first r eigenvalues. ; This indicates that a diagonal matrix is ​​constructed using the elements within the parentheses.

[0105] Random field samples Generated using the following KL (Karhunen-Loève) expansion:

[0106]

[0107] in, This represents the k-th random field sample generated; This is the final mean vector obtained in step S52; Let KL be the coefficient vector of the k-th group, which follows an r-dimensional independent standard normal distribution, i.e. ;in, It is an r-order identity matrix.

[0108] In generating coefficient vectors To improve sampling efficiency and ensure the representativeness of the samples in the probability space, Latin hypercube sampling (LHS) is employed. LHS is applied to independent random variables in the KL expansion, and its specific implementation involves: Each dimension of the 3D random space is divided into equally spaced layers according to probability, ensuring that exactly one sample point is drawn in each layer; first, LHS samples following a uniform distribution of [0,1] are generated, and then they are mapped to a distribution following a uniform distribution of [0,1] using the inverse cumulative distribution function of the standard normal distribution. samples Compared to simple random sampling, this method can effectively avoid sample clustering and ensure the uniformity of sampling space coverage, thereby accurately reproducing the statistical characteristics of the random field with fewer simulations.

[0109] The generated 2n-dimensional sample vector It contains the cohesion and tensile strength values ​​for all n contact nodes. It needs to be separated into two independent parameter fields:

[0110]

[0111]

[0112] in, , These represent the coordinates in the k-th random field implementation. The cohesion and tensile strength values ​​at the contact nodes; This represents the sample vector of the k-th random field realization. The i-th element in This represents the sample vector of the k-th random field realization. The (n+i)th element in the array.

[0113] S54: Parameter Sequence Organization and Derivation:

[0114] Read the contact_info.csv file and generate a linear index sequence based on its original row order (which corresponds one-to-one with the contact IDs). The index sequence is compared with the cohesion random field sample values ​​obtained in step S53. Combined, they form an n x 2 two-dimensional array with the structure: [index i, corresponding cohesion value] Export this array data to a text file c.txt. To enhance readability, write the variable name (e.g., "c") on the first line of the file, and then write the array data on each subsequent line. This file completely stores all the pb_coh (parallel bond cohesion) random field information encountered in this implementation. Using the exact same process as above, combine the index sequence with the tensile strength random field sample values. The data are combined into an array and exported to a file named t.txt, which stores the random field information of pb_ten (parallel bond tensile strength).

[0115] S55: Contact parameter assignment:

[0116] The total number of current contacts is counted through the contact set and stored as an internal variable to initialize the size of the array and table. Sequences of length num are read from c.txt and t.txt and stored in the array. Two lookup tables are created: Table 1 stores the ID-c(i) correspondence and Table 2 stores the ID-t(i) correspondence. All contacts are traversed, and each contact's stable unique identifier ID is used as the key and its corresponding value in the parameter array is used as the value to fill the lookup table. All contacts are traversed again, and their stable IDs are read. The corresponding random field values ​​are queried in Table 1 and Table 2 with these IDs, and the query results are assigned to the pb_coh and pb_ten attributes of the contact, respectively. Finally, the accurate mapping and consistent assignment of the external two-dimensional random field to the discrete element contact parameters are achieved.

[0117] Step S6: Self-weight loading and slope geometry shaping: This step is based on the stress similarity principle of geotechnical centrifuge simulation. By applying a proportionally amplified gravity acceleration field, the self-weight stress state of the prototype high slope is simulated, and staged excavation and equilibrium calculations are performed. The specific implementation process is as follows:

[0118] S61: Gravitational Similarity Criterion and Stress Initialization: Based on the Geometric Dimensions of the Discrete Element Model Compared with the actual slope prototype size Determine the geometric scaling factor based on the proportional relationship. To ensure the stress level inside the model Consistent with the prototype (i.e.) The loading was performed using the "gravity increase method"; particle density was maintained. Consistent with the macroscopic strength parameters and the prototype material, the model environment's gravitational acceleration is used. Set as:

[0119]

[0120] in, Take 9.8 m / s 2 The system uses a servo mechanism or damping iteration to bring the model to an initial static equilibrium state under the gravity field, thereby reproducing the self-weight stress distribution of the real deep slope in a small-scale model.

[0121] S62: Definition of Virtual Inclinometer Deployment and Monitoring:

[0122] To monitor internal slope deformation and identify potential sliding surfaces, a virtual inclinometer monitoring system is deployed at key locations in the model (slope top, slope surface, and slope toe). Specifically, monitoring points are set at 0.5m intervals along the depth direction (y-direction) at the top (x=0), slope surface (x=2.5), and slope toe (x=5). Twenty monitoring points are continuously deployed along the y-direction at the slope top and slope surface, and four monitoring points are continuously deployed at the slope toe. During calculations, the horizontal displacement and vertical depth of particles at each monitoring point are recorded at a fixed time step, generating a deep displacement curve similar to that of an inclinometer in engineering practice. This curve is used to analyze deformation development and the final sliding surface location. The coordinates and spacing mentioned are examples; the actual deployment should be adjusted according to the specific geometric dimensions of the model and the area of ​​interest.

[0123] S63: Slope Cutting and Equilibrium Criteria: The "regional instantaneous deletion method" is adopted, that is, the geometric boundary of the excavation area is defined, and all particles within the boundary are deleted at once; after the slope cutting is completed, the model needs to perform iterative calculations to reach equilibrium again; the convergence criterion for equilibrium is set as: the average unbalanced force ratio of the current system is less than a preset threshold of 1. 10 -5 The slope is considered stable only when the system meets the above convergence conditions, or the displacement rate of the monitoring point approaches zero. If the slope cannot converge within the specified time step and the particle displacement continues to diverge, the slope is considered unstable.

[0124] In a specific application scenario of this embodiment (corresponding to Figures 3(a)-3(d)), the total height of the slope prototype is set to 125 m. Based on the geometric similarity criterion, the model construction height is 12.5 m, and the geometric scaling factor N is 10. Accordingly, the model's gravitational acceleration... The slope was set to 98 m / s² (i.e., 10 × 9.8 m / s²). The slope cutting simulation used a single-stage excavation method with a designed slope angle of 70°. Based on the geometric scaling relationship, the prototype slope cutting depth was 125 m, corresponding to a model slope cutting depth of 12.5 m. The model adopted a linear parallel bond constitutive model, with the following key mesoscopic parameters: particle density of 2500 kg / m³, local damping coefficient of 0.7, and effective contact modulus of... The ratio of normal to tangential stiffness is 1.0, the coefficient of interparticle friction is 0.5, and the average cohesion and tensile strength of parallel bonding are both set to [value missing]. Under this parameter combination, the established undisturbed rock slope model converged within the specified calculation time step and was identified as a stable state. To compare and analyze the differences in failure modes between the homogeneous model and the model considering parameter spatial variability, and to verify the influence of parameter spatial variability on slope stability, the gravitational acceleration was further increased to 490 m / s² until both types of models exhibited the significant failure modes shown in Figures 3(a)-3(d).

[0125] In this embodiment of the invention, by innovatively introducing random field theory into the discrete element modeling process of rock slopes, high-precision simulation of the spatial variability of rock mass parameters is effectively achieved. Figure 2 The illustration shows the main operations and effect comparisons of each stage of the modeling process in this invention. Figure 2 (a) shows the generation of the initial homogeneous particle sample and the setting of physical parameters. Figure 2 (b) Demonstrates the preloading and stress adjustment of the specimen, as well as the extraction and export of contact information (including unique identifier ID and spatial coordinates). Figure 2 (c) illustrates the construction of the external random field and the process of mapping its parameters to the contact points. Figure 2 (d) Figure 2 (e) illustrates the final slope geometry and self-weight loading diagram after parameter assignment. Figure 2 (a)- Figure 2 (d) clearly presents the complete logical structure and implementation path of this invention, from the construction of the original sample and the random assignment of parameters in the space to the numerical simulation. Meanwhile, the comparison of homogeneous and spatially variable slope effects shown in the figure clearly demonstrates: Figure 2 (c)- Figure 2 (e) Together, they demonstrate the distribution characteristics of contact bond strength (contact pb_coh). The homogeneous model (corresponding to the uniform blue area) exhibits a single fixed value, failing to reflect the inherent randomness and weak interlayer distribution characteristics of real rock masses. In contrast, the heterogeneous model, after random field assignment (corresponding to areas with multi-colored banded fluctuations in red, yellow, and blue), reveals multi-scale spatial differences in parameters and a significant correlation structure. This comparison greatly enhances the ability to characterize the heterogeneous structure, weakened zones, and locally vulnerable areas of rock masses, further validating the significant advantages of this method. The above comparison fully demonstrates the beneficial effects of this invention in improving the realism and reliability of slope numerical models.

[0126] Figures 3(a)-3(d) show the comparative results of numerical simulations of slope excavation based on random field theory and mean value model under 50 times gravity conditions, systematically verifying the advantages of the present invention from three aspects: intensity distribution, failure mode, and displacement response. In Figure 3(a), the contact intensity presented by the random field model has obvious strip-like and block-like weakened structures, and the failure path is strictly controlled by the low-intensity area, developing along the weak zone. In contrast, in the mean value model of Figure 3(b), due to the complete homogenization of the intensity parameters and the lack of weak zones in the actual rock mass for guidance, the resulting failure line is singular, localized, and significantly smaller in scale. Furthermore, the contrast between Figures 3(c) and 3(d) is even more profound in terms of displacement field response: in the random field model, the presence of local weak zones triggers large-scale chain slippage, the failure zone exhibits multi-stage expansion characteristics, and ultimately presents an overall unstable state; while the displacement of the mean value model is concentrated in a limited range near the toe of the slope, and the overall degree of failure is significantly weaker. This difference indicates that, under the same geometric and load conditions, the mean model systematically overestimates slope stability by eliminating the actual heterogeneous weakened regions in the rock mass; while the random field model can more realistically reflect the failure sensitivity and instability scale in heterogeneous rock masses.

[0127] This invention innovatively introduces two-dimensional binary random field theory into discrete element modeling of rock slopes, achieving refined simulation of the spatial variability of key mechanical parameters such as rock mass cohesion and tensile strength at the discrete element contact level. Its core lies in constructing autocovariance and cross-covariance matrices using covariance kernel functions, and efficiently generating random field samples by combining Karhunen-Loève (KL) expansion and Latin hypercube sampling (LHS). This allows for simultaneous consideration of parameter mean, variance, coefficient of variation, correlation length, and inter-parameter correlation at the contact level, establishing an integrated modeling method for multi-parameter spatial randomness. To fundamentally overcome the parameter mapping relationship failure and "parameter drift" problems caused by the reliance on real-time spatial coordinate assignment in traditional random field-DEM coupling methods during large deformations or contact reconstruction, this invention proposes a parameter stabilization assignment mechanism based on unique identifiers (IDs). This mechanism assigns a stable and unique identifier (ID) to each contact element, independent of the DEM's internal numbering. An ID-spatial coordinate mapping table is established and exported to an external program to generate a random field, which is then written back to the discrete element model according to the ID-parameter value mapping table. This ID-based lookup and assignment method ensures that even under complex conditions such as self-weight loading, excavation unloading, or earthquakes, and despite significant large deformations or contact reconstruction, each contact element maintains a precise correspondence with its initial random field parameters, achieving accurate parameter mapping and full-process tracking. Regarding the random field generation method, this invention flexibly constructs a two-dimensional binary covariance matrix and a cross-covariance matrix using a covariance kernel. If measured data from boreholes, monitoring, etc., are available, a conditional random field can be further constructed using the Gaussian conditional formula, ensuring that the random field strictly matches the observed values ​​at the measurement points and maintains a reasonable spatial variation structure in unmeasured areas. This is the first time that DEM contact parameters possess the capability of random field modeling with "multi-parameter correlation + conditional update."

[0128] Application results show that the method of this invention enables the microscopic parameters of the slope to exhibit a continuously correlated random distribution in space. As shown in Figure 3(a), under different random field implementations, the initiation location of failure, the crack propagation path, and the final slip surface morphology all exhibit significant randomness and diversity, and can naturally simulate complex instability modes such as preferential failure of local weak zones and bending and bifurcation of slip zones. Compared with homogeneous models that can only produce a single regular slip surface, this invention is closer to the gradual failure and multi-path evolution process of actual rock slopes, and promotes the failure process dominated by the actual weak surface from a physical mechanism perspective. It effectively overcomes the shortcomings of mean models and provides a more realistic prediction of instability morphology and scale for slope stability assessment.

[0129] The various embodiments in this specification are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.

[0130] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of protection of the present invention.

Claims

1. A discrete element modeling method for rock slopes based on random field theory, characterized in that, Includes the following steps: Step S1: Generate an initial homogeneous particle sample based on the set physical parameters and geometric boundaries; Step S2: Apply the target confining pressure to the initial homogeneous particle sample, and use servo control to make the sample reach static equilibrium under the target confining pressure; Step S3: Reset the state of the specimen after stress initialization and assign bonding properties to the contact between particles to construct a homogeneous rock model; Step S4: Generate and export identification information containing a unique identifier and spatial location coordinates for each contact in the homogeneous rock model; Step S5: Based on the identification information, construct a random field characterizing the spatial variability of cohesion and tensile strength parameters, and assign corresponding cohesion and tensile strength parameters to each contact; Step S6: Apply a scaled-up gravity field to the model with the corresponding parameters to simulate the prototype stress, and perform a staged excavation simulation until the model reaches a new equilibrium state after excavation, and obtain displacement monitoring data for stability analysis.

2. The discrete element modeling method for rock slopes based on random field theory according to claim 1, characterized in that, Step S1 includes: Based on the set physical parameters and geometric boundaries, a randomly distributed set of particles is generated within the boundary; a contact model and parameters are specified for the contact between particles in the set, and static equilibrium is achieved through iterative cycles and kinetic energy reduction; particles located outside the boundary are deleted to form an initial homogeneous particle sample.

3. The discrete element method for modeling rock slopes based on random field theory according to claim 1, characterized in that, Step S2 includes: The average normal stress on the model boundary is used as the monitoring variable, and the moving speed of the boundary wall is used as the control variable. Based on the difference between the target confining pressure and the monitoring variable, combined with the servo gain parameter, the control variable is adjusted to make the monitoring variable approach the target confining pressure. The servo gain parameter is determined based on the model's total normal contact stiffness in the pressure direction, the geometric dimensions of the boundary wall, and the calculation time step.

4. The discrete element modeling method for rock slopes based on random field theory according to claim 3, characterized in that, The criteria for determining the static equilibrium state include: the relative error between the monitored variable and the target confining pressure is less than a preset threshold, and the unbalanced force ratio of the particle system within the model is less than a preset threshold or the control variable approaches zero.

5. The discrete element method for modeling rock slopes based on random field theory according to claim 1, characterized in that, Step S3 includes: zeroing the displacement and rotation angle of all particles in the sample; defining a combined contact model and its parameters for all contacts, including linear and parallel bonded parts.

6. The discrete element method for modeling rock slopes based on random field theory according to claim 1, characterized in that, Step S4 includes: traversing all contacts in the model, assigning a unique identifier to each contact and recording its spatial coordinates; and writing the identifier and spatial coordinates into a data file.

7. The discrete element modeling method for rock slopes based on random field theory according to claim 1, characterized in that, Step S5 includes: Based on unique identifiers and spatial coordinates, autocovariance kernels and cross-covariance kernels characterizing the spatial correlation between rock mass cohesion and tensile strength are constructed to build a covariance matrix. Based on monitoring data, the parameter mean vector and covariance matrix are updated by constructing a conditional random field. If there is no monitoring data, the unconditional mean and covariance matrix are used directly.

8. The discrete element method for modeling rock slopes based on random field theory according to claim 7, characterized in that, In step S5, the covariance matrix is ​​decomposed into eigenvalues, the number of truncated modes is determined based on the cumulative contribution rate of eigenvalues, and a random field sample vector containing the cohesion and tensile strength values ​​of each contact node is generated using KL expansion.

9. The discrete element modeling method for rock slopes based on random field theory according to claim 8, characterized in that, In step S5, the random field sample vector is split into a cohesion parameter field and a tensile strength parameter field, and then organized with the unique identifier of the contact and exported as a parameter file. The parameter file is read, the corresponding parameter value is found according to the unique identifier of the contact, and the bonding strength parameter attribute of the corresponding contact in the discrete element model is assigned respectively.

10. The discrete element method for modeling rock slopes based on random field theory according to claim 1, characterized in that, Step S6 includes: The increased gravitational acceleration is set according to the geometric scaling factor of the model and prototype sizes, so that the model reaches initial static equilibrium under the corresponding gravitational field. Virtual monitoring points are set up at key locations on the slope to record the displacement and depth data of particles at each monitoring point; By defining the boundary of the excavation area and deleting the particles within it, iterative calculations are performed. The equilibrium convergence criterion is that the unbalanced force ratio of the model system is less than a preset threshold or the displacement rate of the monitoring point approaches zero. This results in the acquisition of the new equilibrium state and displacement monitoring data.

Images