Rock mass local frictional contact modeling method based on fracture-contact self-adaption

Abstract

The invention discloses a fracture-contact self-adaption-based rock mass local frictional contact modeling method and system. The method comprises the following steps: completing mathematical pre-modeling of a discontinuous surface by adopting a mass point volume representation-based non-penetrating constraint and frictional contact method; constructing a double-domain KD-tree data structure of a local contact domain-non-local action domain, dividing double domains of each calculation mass point through a KD-tree index, updating mass point coordinates in each step to trigger KD-tree local reconstruction, and completing double-domain parallel search; for mass point pairs at the two ends of the fracture surface, the contact force density is decomposed into short-distance repulsive force and friction force, and fracture and contact behaviors are coupled; the method comprises the following steps: simplifying the contact of splashing mass points with disconnected connecting keys into one-dimensional spring vibration damping control, limiting the speed during momentum updating, calculating a damping coefficient by combining rock mass characteristics, supplementing reverse damping contact force, ensuring mass conservation and optimizing collision clastic simulation; and a more reliable simulation method is provided.

Description

Technical Field

[0001] This invention belongs to the field of rock mechanics numerical simulation and engineering disaster prevention and control technology, and more specifically, relates to a rock mass local friction contact modeling method based on fracture-contact adaptive modeling. Background Technology

[0002] In rock engineering construction and operation, rock masses, as natural geological bodies, generally possess joints, fissures, and various structural surfaces. These weak surfaces directly control the stability of the rock mass. In engineering practice, rock masses often face complex physical and mechanical processes such as fracturing, frictional slippage of structural surfaces, and multi-body contact and collision of fractured rock fragments. The coupling effect of these processes often triggers major engineering disasters such as slope instability, tunnel collapse, and rockfall, which not only threaten the lives of construction workers but also cause huge economic losses and ecological damage. Therefore, accurately characterizing these complex processes and achieving effective simulation and prediction of rock engineering disasters is a core requirement for ensuring engineering safety.

[0003] Traditionally, continuum mechanics (CCM) has been the mainstream method for analyzing the mechanical behavior of rock masses, describing the mechanical response of rock mass micro-elements by solving differential equations. However, this method relies on the core assumption of "continuous medium." When the rock mass fractures and forms discontinuities, its governing equations exhibit singularity problems with no solution. To address this issue, existing techniques often require pre-setting crack propagation paths and introducing additional criteria to simulate crack propagation. This results in an inability to realistically reflect the spontaneous initiation and random propagation of cracks in engineering projects, making it difficult to meet the simulation requirements of complex rock mass failure processes.

[0004] The emergence of peri-field dynamics (PD) theory has provided a new approach to solving discontinuity problems. Based on the idea of ​​nonlocal interactions, it uses spatial integral equations to construct models without relying on continuity assumptions. It can naturally and accurately describe the discontinuous failure process of rock masses, while combining the advantages of molecular dynamics and meshless methods. It breaks through the limitations of classical molecular dynamics in terms of computational scale and can be applied to rock mechanics analysis at multiple scales, from macroscopic to microscopic, demonstrating significant advantages in the field of rock mass fracture modeling.

[0005] However, existing near-field dynamic models still have significant shortcomings when applied to macroscopic rock engineering. On the one hand, they are insufficient in simulating the frictional contact of weak surfaces such as rock joints and fissures, making it difficult to accurately reproduce the slip characteristics of structural surfaces. On the other hand, for the large amount of debris formed after rock fracturing, the model cannot effectively handle multi-body contact collisions between debris, and cannot track the contact effects of newly formed interfaces in real time. This results in a discontinuity in the simulation of the entire "fracture-contact" process, failing to fully present the entire process of engineering disasters from occurrence to development.

[0006] While the Discrete Element Method (DEM) has advantages in simulating particle contact and can effectively describe friction and collision behavior, it struggles to account for the nonlocal mechanical properties of rock fracture processes. Therefore, integrating the advantages of different methods to achieve accurate simulation of the entire "fracture-friction-contact" chain process in rock masses, overcoming the challenge of unified modeling of discontinuous rock failure and complex contact behavior, and providing reliable predictive data for rock engineering disasters, has become a critical issue urgently needing to be addressed in the field of rock engineering to ensure engineering safety and reduce disaster risks. This has significant practical engineering value and scientific significance. Summary of the Invention

[0007] This invention aims to solve the core challenge of unifying "local contact-non-local fracture" in existing rock engineering modeling. Addressing the coupled processes of rock mass fracture, structural surface friction, and multi-body contact of aggregates, it overcomes the limitation of traditional models that cannot fully simulate the entire "fracture-contact" process. By integrating relevant theories and algorithms, it achieves accurate modeling of the mechanical characterization of weak rock surfaces, dynamic updating of contact surfaces, and control of splash particles. This provides a more reliable simulation method for engineering disasters such as slope instability, helping to improve the safety assurance capabilities of rock engineering projects.

[0008] To address the aforementioned deficiencies or improvement needs of existing technologies, as a first aspect of this invention, the present invention provides a method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling, comprising: S1. Using a non-penetrating constraint and frictional contact method based on particle volume characterization, the weak surfaces of rock mass, including joints, fissures and structural planes, are mechanically characterized, and mathematical pre-modeling of discontinuous surfaces is achieved by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. S2. Construct a dual-domain KD-tree data structure of local contact domain and non-local scope domain; for each computational mass point, divide the local contact domain and non-local scope domain by KD-tree index; at each time step, only update the spatial coordinates of the mass point to trigger local reconstruction of the KD-tree and complete the parallel search of the dual domains. S3. For the particle pairs at both ends of the fracture surface, the contact force density is decomposed into short-range repulsive force and frictional force density, thus completing the coupling of fracture and contact behavior; S4. For splash particles whose connecting keys are completely broken, simplify their contact with other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, first limit the particle velocity, and then supplement the damping contact force opposite to the direction of motion by combining the damping coefficient dynamically calculated by rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

[0009] Furthermore, the frictional contact method in S1 is specifically as follows: In the discrete element contact model of a sphere, the contact forces between particle pairs within the local contact domain are defined, including normal contact forces and tangential frictional forces, to simulate the mechanical interactions at weak rock surfaces such as joints, fissures, and structural planes; among them, the normal contact force... The calculation is performed using the following formula: , in, For point mass The normal contact force acting on it; The point mass under study is the object subjected to the normal contact force. For key connection and The relative position vector between particles; For contact stiffness; The size of the nonlocal scope; For point mass Outer normal direction; It is a nonlocal scope; This is a local contact domain.

[0010] Furthermore, the method for calculating the tangential friction force is as follows: , in, For point mass The tangential frictional force experienced; Normal contact force Size; It is the tensor product; The coefficient of friction; The velocity of the particle; It is a unit matrix.

[0011] Furthermore, the method for calculating the contact force density in S3 is as follows: , in, This represents the total contact force density. This represents the short-range repulsive force density. Friction density; This represents the relative position vector between particles connected by a bond.

[0012] Furthermore, the method for calculating the short-range repulsive force in S3 is as follows: , in, This represents the short-range repulsive force density. For contact stiffness; The embedding distance of the particle to the spherical characterization volume range; This is the lower bond position vector for the current configuration; It is a unit direction vector; , in, For embedding displacement vector The length of the module.

[0013] Furthermore, the method for calculating the frictional force in S3 is as follows: , in, Friction density; The coefficient of friction; This is a sign function used to determine the trend and direction of relative motion between particles.

[0014] Furthermore, the damping coefficient in S4 is calculated as follows: , in, The damping coefficient; The damping ratio; The density of the rock mass; For point mass The contact area centered on the target; For point mass With particles within the scope Spatial distance between them; Parameters related to the equivalent contact area of ​​a particle; The characteristic dimension representing the volume of a particle; For point mass The characterization of the volumetric infinitesimal element.

[0015] Furthermore, the damping contact force in S4 that is opposite to the direction of motion is: , in, For damping contact force; This is the fundamental force acting when splashing particles come into contact; The direction of motion of the splashing particles; is the velocity vector of the splashing particles.

[0016] As a second aspect of the present invention, the present invention provides a rock mass local friction contact modeling system based on fracture-contact adaptive modeling, comprising: The rock mass weak surface mathematical pre-modeling unit is used to mechanically characterize the rock mass weak surface, including joints, fissures and structural surfaces, using non-penetrating constraints and frictional contact methods based on particle volume representation. It also realizes mathematical pre-modeling of discontinuous surfaces by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. The dual-domain KD-tree construction and search unit is used to construct a dual-domain KD-tree data structure of local contact domain and non-local scope. For each computational mass point, the local contact domain and non-local scope are divided by KD-tree index. At each time step, the local reconstruction of KD-tree is triggered by updating the spatial coordinates of the mass point, thus completing the parallel search of the dual domain. The fracture-contact behavior coupling modeling unit is used to decompose the contact force density into short-range repulsive force and friction force density for the mass pairs at both ends of the fracture surface, thus completing the coupling between fracture and contact behavior. The splash particle optimization unit is used to simplify the contact between splash particles with completely broken connecting keys and other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, the particle velocity is first limited, and then the damping contact force opposite to the direction of motion is supplemented by the damping coefficient dynamically calculated by combining the rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

[0017] As a third aspect of the invention, the invention provides a computer-readable storage medium having a computer program stored thereon, the computer program being executed by a processor of any step of the described method for modeling local frictional contact in rock mass based on fracture-contact adaptation.

[0018] In summary, compared with the prior art, the above-described technical solutions conceived by this invention can achieve the following beneficial effects: 1. The rock mass local friction contact modeling method based on fracture-contact adaptation of the present invention constructs a dual-domain KD-tree data structure of local contact domain and non-local action domain. For each computational particle, the dual-domain partitioning is completed using the KD-tree index, and each time step triggers local KD-tree reconstruction only by updating the spatial coordinates of the particle, achieving parallel search of the dual domains. This technique avoids the redundant computation of traditional full-domain search, significantly improving the efficiency of particle action domain partitioning and search. Simultaneously, dynamic local reconstruction ensures that the dual-domain range always matches the actual action state when the spatial position of the particle changes during rock mass deformation and fracture, providing accurate and efficient spatial index support for subsequent fracture-contact behavior modeling, and guaranteeing the real-time performance and accuracy of the entire calculation process.

[0019] 2. The rock mass local friction contact modeling method based on fracture-contact adaptation of the present invention decomposes the contact force density into short-range repulsive force and friction force for the mass pairs at both ends of the fracture surface, thereby achieving coupled modeling of fracture and contact behavior. This technology overcomes the limitations of traditional models that simulate fracture and contact behavior separately, accurately depicting the true mechanical state of the mass pairs at both ends of the fracture surface during rock mass fracture, where there is both short-range repulsive force due to relative embedding and frictional force due to relative sliding. This makes the simulation results more consistent with the actual fracture-contact response law of rock mass, effectively solving the problem of mechanical behavior distortion caused by separate modeling, and improving the reliability of the simulation of fracture and contact coupled behavior during rock mass disasters.

[0020] 3. The rock mass local friction contact modeling method based on fracture-contact adaptation of the present invention simplifies the contact between splash particles with completely broken connecting bonds and other particles or the main model as a one-dimensional spring vibration damping control problem. During momentum update, the particle velocity is first constrained, and then the damping coefficient is dynamically calculated based on rock mass characteristics to supplement the damping contact force opposite to the direction of motion. This technology specifically solves the problem of non-conservation of model mass caused by the unconstrained motion of splash particles. Simultaneously, by simulating the collision behavior of splash particles through a spring-damping system, it optimizes the simulation effect of rock mass collision fragmentation, avoiding the phenomenon of false bouncing or uncontrolled motion of splash particles in traditional simulations. This makes the simulation of debris movement and accumulation after rock mass fracture more consistent with engineering reality, further improving the simulation accuracy of the entire rock mass disaster process. Attached Figure Description

[0021] Figure 1 This is a flowchart of the rock mass local friction contact modeling method based on fracture-contact adaptive method according to an embodiment of the present invention; Figure 2 This is a schematic diagram of a partial contact model according to an embodiment of the present invention; Figure 3 This is a schematic diagram of the outward normal direction of a particle according to an embodiment of the present invention; Figure 4 This is a schematic diagram of the discrete element contact model of a sphere according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the local search of the particle contact domain using the KDTree algorithm according to an embodiment of the present invention; Figure 6 This is a logical relationship diagram between the fracture contact adaptive method, the nonlocal PD model, and the local contact model in this embodiment of the invention. Figure 7 This is a schematic diagram of a splash particle control method according to an embodiment of the present invention; Figure 8 This is a schematic diagram of a uniaxial compression model of a fractured material according to an embodiment of the present invention; Figure 9This is a schematic diagram comparing the samples before and after damage according to an embodiment of the present invention; Figure 10 This is a uniaxial compressive force-displacement curve of the fractured material according to an embodiment of the present invention; Figure 11 This is a system unit diagram of an embodiment of the present invention. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0023] Example 1 Please refer to Figure 1 This embodiment 1 provides a method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling, including: S1. Using a non-penetrating constraint and frictional contact method based on particle volume characterization, the weak surfaces of rock mass, including joints, fissures and structural planes, are mechanically characterized, and mathematical pre-modeling of discontinuous surfaces is achieved by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. S2. Construct a dual-domain KD-tree data structure of local contact domain and non-local scope domain; for each computational mass point, divide the local contact domain and non-local scope domain by KD-tree index; at each time step, only update the spatial coordinates of the mass point to trigger local reconstruction of the KD-tree and complete the parallel search of the dual domains. S3. For the particle pairs at both ends of the fracture surface, the contact force density is decomposed into short-range repulsive force and frictional force density, thus completing the coupling of fracture and contact behavior; S4. For splash particles whose connecting keys are completely broken, simplify their contact with other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, first limit the particle velocity, and then supplement the damping contact force opposite to the direction of motion by combining the damping coefficient dynamically calculated by rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

[0024] This embodiment 1 further elaborates on the above steps.

[0025] Peridynamics (PD) is an emerging method that uses nonlocal interactions to model and solve spatial integral equations to describe the mechanical behavior of matter. It combines the advantages of molecular dynamics and meshless methods, avoiding the singularity of traditional macroscopic methods based on continuity assumptions and solving spatial differential equations when faced with discontinuous problems. It also overcomes the computational limitations of classical molecular dynamics, demonstrating high accuracy and efficiency in analyzing both macroscopic and microscopic discontinuous mechanical problems. The shift in PD from local contact stress (surface force density) to nonlocal, non-contact bond forces represents a significant development in the modeling paradigm of continuum mechanics. Specifically: Continuum mechanics: Continuum mechanics obtains the distribution of infinitesimal elements in space and their changes over time through differential forms, as well as the mechanical responses resulting from these changes. The governing equations are undefined at discontinuities, and solving them at these discontinuities can lead to singularities. These singularities usually require the introduction of additional algorithms to handle, such as pre-defining the crack propagation path and introducing additional criteria to simulate the crack propagation behavior.

[0026] Perifield dynamics: Perifield dynamics reconstructs the theoretical framework of traditional continuum mechanics using spatial integrals. Unlike continuum mechanics, the integral governing equations of perifield dynamics remain defined at discontinuities, thus handling discontinuity issues more naturally. This allows perifield dynamics to better describe phenomena such as the spontaneous initiation and propagation of cracks in materials.

[0027] Strain energy density of particle X based on PD model Compared with the strain energy density of CCM micro-element in traditional continuum mechanics The equivalence principle connects perifield dynamics theory with traditional continuum mechanics methods. When the perifield size approaches zero, perifield dynamics can converge to classical continuum mechanics; therefore, PD is a useful extension and supplement to CCM.

[0028] In a preferred embodiment, the specific process of traditional continuum mechanics is associated with the principle of strain energy density equivalence as follows: In conventional continuum mechanics (CCM), the strain energy density of a rock micro-element under stress is given by... In bonded near-field dynamics BB-PD, the strain energy density of the particle at the corresponding location is: ; Based on the principle that the strain energy density of the same rock under the same stress is unique, let the two be equal: , For linearly elastic rocks, in CCM: , in, These are components of the stress tensor; These are the components of the strain tensor; In BB-PD: , in, For is a point mass The characterization volume; For is key stiffness; For is key Elongation; By combining the equivalent formulas, we can derive: , in, It is a fourth-order elastic constant tensor; strain tensor components For tensor indices, the value is... Corresponding spatial rectangular coordinate system direction; By using the above equivalence relations, the "bond stiffness" required by the PD model can be calculated in reverse by directly using the existing CCM parameters in the project, without having to carry out complex parameter calibration tests separately for the PD model.

[0029] (1) Mathematical pre-modeling of weak surfaces in rock mass Please refer to Figure 2 For a brittle material system involving complex failure processes such as fracture, friction, and multi-body contact, this paper attempts to achieve a unified modeling of local contact and nonlocal fracture at a more macroscopic engineering scale, based on nonlocal PD theory and incorporating a contact model using the spherical discrete element method and the KDTree contact search algorithm. First, for the fracture process of brittle materials, nonlocal PD theory (BB elastic-brittle or SB elastic-plastic) is used for modeling. Second, for the frictional contact problem of discontinuous surfaces, non-penetrating constraints and a discrete element frictional contact model are introduced. Finally, for the multi-body contact problem caused by material fracture, an efficient contact search algorithm and multi-body contact model based on updating the particle positions are developed using the KDtree algorithm.

[0030] Please refer to Figure 3 The contact interface is characterized within the PD model framework by calculating the outward normal direction of the mass point. And mark the object The volume is represented by the outer contour particles, and non-penetrating constraints are set to perform mathematical pre-modeling of discontinuous surfaces.

[0031] In a preferred embodiment, the method for calculating the direction of the outward normal of the particle is as follows: , in, It is a point mass Scope internal particles Deformation vector exist Component of direction; Represents a point mass Deformation vector in Components in direction, and Correspondingly.

[0032] The relative spacing between particle pairs is calculated by representing the volume of the particles, and a penalty function is introduced to set a non-penetrating constraint on the particle pairs; in a preferred embodiment, the non-penetrating constraint based on the particle volume representation is specifically as follows: in, Indicators representing the criteria for non-penetration constraints. It is a point mass Local contact domain internal particles Deformation vector, The distance between particles; Represents a point mass The position vector describes the point mass. Its position in space; It indicates the direction of the outward normal of the particle.

[0033] Simultaneously, a frictional contact method is introduced; please refer to [reference needed]. Figure 4 In the discrete element contact model of a sphere, the contact force is defined for the particle pairs within the local contact domain: the normal contact force comprehensively considers the relevant parameters between particles, material mechanical parameters, characteristic values ​​of the action range, and the direction of the outward normal, and is obtained by summing the effects of the particles within the local contact domain, and is used to simulate the squeezing effect at the weak surface; the tangential friction force is related to the magnitude of the normal contact force, the friction coefficient, and the direction of the relative velocity, and the tangential component is extracted through tensor operations to ensure that it is along the tangential direction of the relative motion, so as to restore the friction effect at the weak surface.

[0034] In a preferred embodiment, the friction contact method specifically includes: In the discrete element contact model of a sphere, the contact forces between particle pairs within the local contact domain are defined, including normal contact forces and tangential frictional forces, to simulate the mechanical interactions at weak rock surfaces such as joints, fissures, and structural planes; among them, the normal contact force... The calculation is performed using the following formula: , in, For point mass The normal contact force acting on it; The point mass under study is the object subjected to the normal contact force. For key connection and The relative position vector between particles; For contact stiffness; The size of the nonlocal scope; For point mass Outer normal direction; It is a nonlocal scope; This is a local contact domain.

[0035] In a preferred embodiment, the tangential friction force is calculated as follows: , in, For point mass The tangential frictional force experienced; Normal contact force Size; It is the tensor product; The coefficient of friction; The velocity of the particle; It is a unit matrix.

[0036] By calculating the outward normal direction of the mass and the volume of the mass on the outer contour of the object, a mathematical pre-modeling of the discontinuous surface is achieved, and the mechanical characterization of the weak surface is completed by combining the two methods mentioned above.

[0037] This process solves the problem of insufficient characterization of weak surfaces in traditional models, enabling the model to realistically reflect the normal constraint and tangential friction characteristics of weak surfaces under stress. It provides a reliable weak surface characterization basis for accurately simulating the fracture path, strength degradation and overall mechanical response of rocks with weak surfaces, and also makes subsequent rock fracture simulation more consistent with the mechanical behavior of rock masses in actual engineering.

[0038] (2) Construction and Search of Two-Domain KD-tree Please refer to Figure 5 In fact, nonlocal PD models and local friction contact models can only solve contact problems with pre-defined boundaries. However, in actual macroscopic engineering problems, it is more necessary to consider the catastrophic mechanics problems after failure. The aforementioned local friction contact model requires the outer boundary of the object to be determined in advance. For catastrophic mechanics problems, the newly exposed fracture interface under the crack path also needs to participate in the contact calculation in the next time step. Here, we developed an efficient algorithm for local search of the particle contact domain and nonlocal search of the action domain, based on the KDTree algorithm, for the new contact surface formed during the fracture process. This algorithm achieves real-time tracking and judgment of multi-body contact by updating the particle position alone.

[0039] By constructing a KD-tree data structure from the coordinates of all mass points, a spatial index is established using a recursive partitioning method, so that each node corresponds to a specific subset of mass points. Based on this, two types of scopes are defined for each mass point: a local contact domain is formed by searching for mass points within a short distance based on a set local contact domain radius; and a non-local scope is formed by searching for mass points at a greater distance based on a larger non-local scope radius. The distance is determined using the Euclidean norm.

[0040] In a preferred embodiment, the specific method for dividing local contact regions and non-local scope regions using a KD-tree index is as follows: Set of coordinates of all particles Constructing a KD-tree, where each node corresponds to a subset of data points, and implementing spatial indexing through recursive partitioning, denoted as […]. .

[0041] Given the radius of the local contact domain For a point mass KD-tree search satisfies All particles , forming a local contact domain ,Right now: , in, To indicate the current research topic The spatial coordinates of a point mass; Represents a set any point in Spatial coordinates; Given the radius of the nonlocal scope , KD-tree search satisfies All particles , constitute a nonlocal scope ,Right now: , in, It is the Euclidean norm, used to calculate point masses. and point mass The spatial distance between them.

[0042] At each time step, the total force acting on the particle is first calculated—consisting of the contact forces within the local contact domain and the bond forces within the non-local action domain. Then, the acceleration is derived according to Newton's second law, and time integration is performed using the central difference method to update the particle's velocity and coordinates. Boundary particles are adjusted according to the constraint conditions; fixed-constraint particles maintain their coordinates, while displacement-constrained particles are updated according to preset increments.

[0043] In a preferred embodiment, the method for updating the spatial coordinates of the mass point is as follows: Let the mass be a point At any moment The coordinates are The quality is ; point mass Total force Local contact force Nonlocal bond forces vector sum: , , , in, It is a local contact force; These are nonlocal bond forces; For local contact domain; It is a nonlocal scope; This represents the total contact force density. It is the nonlocal bond force density; Let be an integral infinitesimal element, representing a point mass within the domain of action. Characterizing volume elements; According to Newton's second law, a point mass... At any moment acceleration for: , The central difference method is used for time integration, with a time step of . : , , in, show mass point At any moment speed; The current moment is a time variable; For point mass At any moment speed; For point mass At any moment Spatial coordinates; For a boundary mass point, if it is subject to fixed constraints, then If subject to displacement constraints, the adjustment will be made according to the preset displacement amount, i.e. ,in This is the preset displacement increment.

[0044] After the coordinates are updated, the KD-tree automatically triggers local reconstruction, enabling parallel search across two domains and real-time tracking of the multi-body contact state of the rock mass. This provides a precise spatial relationship basis for subsequent fracture-contact behavior coupling and special particle control.

[0045] (3) Coupled modeling of fracture-contact behavior Please refer to Figure 6 To address the complex multi-body contact problem arising after fracture, this paper leverages the advantage of the clear physical entity representation of the particle-based method. The boundaries of multiple objects are characterized using particle normals and volume occupation. Based on this, the aforementioned frictional contact algorithm is introduced to solve the multi-body contact problem. Within the local contact domain, the contact forces and frictional forces between particle pairs at both ends of the fracture surface are directly calculated as force density through the broken bonds. By cascading the fracture-contact adaptive methods, a true unification of local contact and non-local fracture modeling can be achieved.

[0046] For the particle pairs at both ends of the fracture surface, the contact force density is decomposed into two parts: short-range repulsive force density and frictional force density. The two together constitute the total contact force density, thus fully characterizing the mechanical effect during contact.

[0047] In a preferred embodiment, the contact force density is calculated as follows: , in, This represents the total contact force density. This represents the short-range repulsive force density. Friction density; This represents the relative position vector between particles connected by a bond.

[0048] Among them, the short-range repulsive force is used to resist the mutual embedding between particles. Its calculation needs to be based on the embedding distance of the particles to the spherical characterization volume range: when the modulus of the embedding displacement vector (i.e. the actual embedding distance) is less than the set embedding distance threshold, the short-range repulsive force density is related to the difference between the "embedding distance threshold minus the actual embedding distance", and the direction is determined by the unit direction vector, which is consistent with the direction of the embedding displacement vector; if the actual embedding distance reaches or exceeds the threshold, the short-range repulsive force density is zero and no longer produces a repulsive effect.

[0049] In a preferred embodiment, the method for calculating short-range repulsive force is as follows: , in, This represents the short-range repulsive force density. For contact stiffness; The embedding distance of the particle to the spherical characterization volume range; This is the lower bond position vector for the current configuration; It is a unit direction vector; , in, For embedding displacement vector The length of the module.

[0050] Friction is used to impede the relative sliding between particles. Its calculation is based on the short-range repulsive force density, and the magnitude of the friction is controlled by the friction coefficient, which is related to characteristics such as the roughness of the contact surface. Simultaneously, a sign function is used to determine the trend direction of the relative motion between particles—the result of the sign function is determined by the time rate of change of the embedded displacement vector magnitude. This ensures that the direction of the friction is always opposite to the trend direction of the relative motion, guaranteeing that the friction effectively impedes the relative sliding of the particles. In this way, the two core mechanical forces of particle contact after fracture are integrated, achieving coupled modeling of rock mass fracture behavior and contact behavior.

[0051] In a preferred embodiment, the friction force is calculated as follows: , in, Friction density; The coefficient of friction; This is a sign function used to determine the trend and direction of relative motion between particles.

[0052] (4) Splash particle optimization During rock mass catastrophe, the contact behavior of ejected particles whose connecting bonds are completely broken (such as rock fragments generated by blasting or rock pieces detached from the main body when the slope is unstable) after they are freed from their original mechanical constraints directly affects the realism of the catastrophe evolution. If such particles lack effective constraints, they are prone to uncontrolled movement (such as penetrating other structures or scattering without restriction), resulting in an imbalance in the conservation of model mass and failing to reflect the energy dissipation and debris accumulation characteristics of rock collisions in actual engineering.

[0053] Please refer to Figure 7 From a mechanical perspective, the contact between splashing particles is essentially a brief and discontinuous impact. A one-dimensional spring vibration damping control model perfectly captures this characteristic: the spring effect simulates the elastic recovery during contact (such as the deformation and rebound of rock blocks during collision), while the damping effect reflects energy dissipation (such as the conversion of kinetic energy into thermal or fragmentation energy during the collision). This simplification retains the core characteristics of contact mechanics while avoiding the computational redundancy caused by complex multi-directional constraints.

[0054] In actual calculations, the momentum update phase first applies a threshold limit to the velocity of the splash particles, i.e. ,in The maximum speed threshold is set to prevent calculation instability caused by abnormal speed (such as numerical oscillation at high speeds), which is consistent with the reality in engineering where the collision speed of rock blocks is constrained by factors such as material strength and gravity.

[0055] The dynamic calculation of the damping coefficient is deeply related to the physical properties of the rock mass: the density of the rock mass determines the inertial characteristics of the particles, the spacing between particles in the contact domain affects the interaction intensity, and the equivalent contact area and the characterizing volume parameters transform the geometric characteristics of discrete particles into mechanical parameters, so that the damping effect can accurately match the energy dissipation law of different lithologies (such as the brittle collision of granite and the plastic collision of sandstone).

[0056] In a preferred embodiment, the damping coefficient is calculated as follows: , in, The damping coefficient; The damping ratio; The density of the rock mass; For point mass The contact area centered on the target; For point mass With particles within the scope Spatial distance between them; Parameters related to the equivalent contact area of ​​a particle; The characteristic dimension representing the volume of a particle; For point mass The characterization of the volumetric infinitesimal element.

[0057] The resulting damping contact force is opposite to the direction of the particle's motion. In a preferred embodiment, the damping contact force opposite to the direction of motion is: , in, For damping contact force; This is the fundamental force acting when splashing particles come into contact; The direction of motion of the splashing particles; is the velocity vector of the splashing particles.

[0058] It not only constrains the disordered motion of splashed particles through resistance (ensuring mass conservation), but also simulates deceleration and stagnation in real collisions through energy dissipation (such as the rolling deceleration of rock blocks on a slope after blasting and the kinetic energy decay during debris accumulation), enabling the simulation results to more reliably predict the distribution range of debris and key engineering parameters such as impact load after a disaster.

[0059] Meanwhile, the following experimental procedures were carried out based on the method of Embodiment 1: Traditional fracture modeling methods destroy all bonding at the fracture surface or delete material particles near the fracture, which leads to two problems: first, the fracture surface cannot transmit stress or waves, which is barely feasible under tensile-shear stress but not applicable under compressive-shear stress; second, the fracture tip will be in a non-closed state. Under tensile-shear stress, cracks tend to open, and there is no stress on the fracture surface. However, under compressive-shear stress, cracks tend to close, and the contact surface needs to transmit compressive stress or stress waves. Here, we attempt to apply the contact model mentioned earlier to the fracture surface to establish a frictional contact model of the fracture surface.

[0060] Please refer to Figure 8 As shown, there is a height millimeters, width A rectangular specimen, millimeters in diameter, contains randomly distributed cracks. Closed cracks are located inside the specimen and intersect its boundaries. The specimen is discretized into 11,760 material points (MPs) with a spacing between the material points. Millimeters, influence field size millimeter, that is The explicit calculation time step is set to 5e-9 seconds.

[0061] The specific values ​​of each input calculation parameter are listed in Table 1 below.

[0062] , Calculation results of uniaxial compressive failure of fractured materials are as follows Figure 9 As shown, the complete failure process of fractured materials, from initial fragmentation and joint expansion to macroscopic penetration, was successfully reproduced. Figure 8 As shown, the equivalent elastic modulus of the fractured material specimen, calculated from the elastic stage of the force-deformation curve, is 5.65 GPa, slightly lower than the input elastic modulus of the intact brittle material (7.2 GPa). With increasing stress, the specimen first undergoes slip deformation along the two main fracture surfaces (upper right-lower left orientation). The block on the right side of the specimen, cut by two sets of joints, is squeezed and extruded towards the free surface. Simultaneously with the fracture of the right-side block, the cumulative slip deformation on the main fracture surfaces leads to the formation of airfoil tensile cracks at the tips of the two main fracture surfaces. As the load increases, a continuous fracture surface forms along the main fracture direction at the peak stress stage. The continuity and slippage of the fracture network cause a rapid release of stored energy within the specimen, inducing the ejection and splashing of fracture fragments.

[0063] In the post-peak stage, crack propagation gradually becomes interconnected, and the failure mode shifts from tensile-dominated to a shear-tensile hybrid mode, ultimately forming the main shear fracture zone. This is highly consistent with the evolution process from wing cracks to main shear fractures summarized in laboratory tests. In the post-peak stage, the axial load rapidly transfers from the discontinuities caused by crack propagation within the specimen to the contact transfer between the testing machine's loading plate and the specimen; the contact effect becomes the dominant mechanism for stress transfer. On the stress-strain curve, the stress drops rapidly, but the rate of drop gradually slows down, and the stress-strain curve enters a softening plateau that gradually approaches the residual strength. Figure 10 As shown, this corresponds to a hyperbolic attenuation, meaning that after the stress reaches its peak, it gradually and asymptotically approaches the residual strength with strain as the independent variable.

[0064] Example 2 Please refer to Figure 11 This embodiment 2 provides a rock mass local friction contact modeling system based on fracture-contact adaptive modeling, including: The rock mass weak surface mathematical pre-modeling unit is used to mechanically characterize the rock mass weak surface, including joints, fissures and structural surfaces, using non-penetrating constraints and frictional contact methods based on particle volume representation. It also realizes mathematical pre-modeling of discontinuous surfaces by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. The dual-domain KD-tree construction and search unit is used to construct a dual-domain KD-tree data structure of local contact domain and non-local scope. For each computational mass point, the local contact domain and non-local scope are divided by KD-tree index. At each time step, the local reconstruction of KD-tree is triggered by updating the spatial coordinates of the mass point, thus completing the parallel search of the dual domain. The fracture-contact behavior coupling modeling unit is used to decompose the contact force density into short-range repulsive force and friction force density for the mass pairs at both ends of the fracture surface, thus completing the coupling between fracture and contact behavior. The splash particle optimization unit is used to simplify the contact between splash particles with completely broken connecting keys and other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, the particle velocity is first limited, and then the damping contact force opposite to the direction of motion is supplemented by the damping coefficient dynamically calculated by combining the rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

[0065] Example 3 This embodiment 3 also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, can implement any step of a rock mass local friction contact modeling method based on fracture-contact adaptive method.

[0066] The computer-readable storage medium may include various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0067] For a description of the computer-readable storage medium provided in this application, please refer to the above method embodiments; further details will not be repeated here.

[0068] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling, characterized in that, include: S1. Using a non-penetrating constraint and frictional contact method based on particle volume characterization, the weak surfaces of rock mass, including joints, fissures and structural planes, are mechanically characterized, and mathematical pre-modeling of discontinuous surfaces is achieved by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. S2. Construct a dual-domain KD-tree data structure of local contact domain and non-local scope domain; for each computational mass point, divide the local contact domain and non-local scope domain by KD-tree index; at each time step, only update the spatial coordinates of the mass point to trigger local reconstruction of the KD-tree and complete the parallel search of the dual domains. S3. For the particle pairs at both ends of the fracture surface, the contact force density is decomposed into short-range repulsive force and frictional force density, thus completing the coupling of fracture and contact behavior; S4. For splash particles whose connecting keys are completely broken, simplify their contact with other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, first limit the particle velocity, and then supplement the damping contact force opposite to the direction of motion by combining the damping coefficient dynamically calculated by rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

2. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling as described in claim 1, characterized in that, The frictional contact method in S1 is specifically as follows: In the discrete element contact model of a sphere, the contact forces between particle pairs within the local contact domain are defined, including normal contact forces and tangential frictional forces, to simulate the mechanical interactions at weak rock surfaces such as joints, fissures, and structural planes; among them, the normal contact force... The calculation is performed using the following formula: ， in, For point mass The normal contact force acting on it; The point mass under study is the object subjected to the normal contact force. For key connection and The relative position vector between particles; For contact stiffness; The size of the nonlocal scope; For point mass Outer normal direction; It is a nonlocal scope; This is a local contact domain.

3. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling according to claim 2, characterized in that, The method for calculating the tangential friction force is as follows: ， in, For point mass The tangential frictional force experienced; Normal contact force Size; It is the tensor product; The coefficient of friction; The velocity of the particle; It is a unit matrix.

4. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling as described in claim 1, characterized in that, The method for calculating the contact force density in S3 is as follows: ， in, This represents the total contact force density. This represents the short-range repulsive force density. Friction density; This represents the relative position vector between particles connected by a bond.

5. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling according to claim 4, characterized in that, The method for calculating the short-range repulsive force in S3 is as follows: ， in, This represents the short-range repulsive force density. For contact stiffness; The embedding distance of the particle to the spherical characterization volume range; This is the lower bond position vector for the current configuration; It is a unit direction vector; ， in, For embedding displacement vector The length of the module.

6. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling according to claim 5, characterized in that, The method for calculating the friction density in S3 is as follows: ， in, Friction density; The coefficient of friction; This is a sign function used to determine the trend and direction of relative motion between particles.

7. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling according to claim 1, characterized in that, The method for calculating the damping coefficient in S4 is as follows: ， in, The damping coefficient; The damping ratio; The density of the rock mass; For point mass The contact area centered on the target; For point mass With particles within the scope Spatial distance between them; Parameters related to the equivalent contact area of ​​a particle; The characteristic dimension representing the volume of a particle; For point mass The characterization of the volumetric infinitesimal element.

8. The method for modeling local frictional contact in rock mass based on fracture-contact adaptive modeling according to claim 7, characterized in that, The damping contact force in S4 that is opposite to the direction of motion is: ， in, For damping contact force; This is the fundamental force acting when splashing particles come into contact; The direction of motion of the splashing particles; is the velocity vector of the splashing particles.

9. A rock mass local friction contact modeling system based on fracture-contact adaptive modeling, characterized in that, include: The rock mass weak surface mathematical pre-modeling unit is used to mechanically characterize the rock mass weak surface, including joints, fissures and structural surfaces, using non-penetrating constraints and frictional contact methods based on particle volume representation. It also realizes mathematical pre-modeling of discontinuous surfaces by calculating the outward normal direction of particles and the volume of particles on the outer contour of the object. The dual-domain KD-tree construction and search unit is used to construct a dual-domain KD-tree data structure of local contact domain and non-local scope. For each computational mass point, the local contact domain and non-local scope are divided by KD-tree index. At each time step, the local reconstruction of KD-tree is triggered by updating the spatial coordinates of the mass point, thus completing the parallel search of the dual domain. The fracture-contact behavior coupling modeling unit is used to decompose the contact force density into short-range repulsive force and friction force density for the mass pairs at both ends of the fracture surface, thus completing the coupling between fracture and contact behavior. The splash particle optimization unit is used to simplify the contact between splash particles with completely broken connecting keys and other particles or the main model into a one-dimensional spring vibration damping control problem. When updating momentum, the particle velocity is first limited, and then the damping contact force opposite to the direction of motion is supplemented by the damping coefficient dynamically calculated by combining the rock mass characteristics. This ensures the overall mass conservation of the calculation model and optimizes the simulation of rock mass collision fragmentation.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, The computer program is executed by a processor using the rock mass local friction contact modeling method based on fracture-contact adaptation as described in any one of claims 1-8.

Images