A rock failure simulation method based on three-dimensional particle DDA

Abstract

This invention discloses a rock failure simulation method based on three-dimensional particle dynamics analysis (DDA), relating to the field of rock mechanics. The study domain is discretized into a set of multiple rigid spherical particles, and their initial positions, radii, material properties, and control parameters are recorded. A contact determination scheme is used to determine the spatial positions of contact points between particles, and a distance threshold is used to control the contact activation conditions. Based on the principle of minimum potential energy, stiffness sub-matrices corresponding to four types of springs (normal, shear, torsional, and rolling) and friction force sub-matrices corresponding to frictional forces are generated. These matrices are assembled into a system of simultaneous equations, and the correct contact state is determined using open-closed iteration. During this process, a set of contact states is first predicted, the simultaneous equations are solved, and the solved contact states are updated. If the solved contact states are consistent with the predicted contact states, the iteration converges; otherwise, the updated contact states are used to update the predicted contact states, and this process is repeated until convergence. This solves the problem of maintaining consistent macroscopic mechanical responses at different particle scales.

Description

Technical Field

[0001] This invention relates to the field of rock mechanics, and more specifically to a method for simulating rock failure based on three-dimensional particle DDA. Background Technology

[0002] In the field of rock mechanics, numerical simulation methods are mainly divided into two categories: continuous methods and discontinuous methods. Traditional continuous medium methods, such as the Finite Element Method (FEM), Finite Difference Method (FDM), and Finite Volume Method (FVM), have formed a mature system for simulating intact rock masses or small deformation problems. However, when faced with complex contact problems such as highly discontinuous materials, large displacements, multiple crack initiation and propagation, and particle flow, the assumption of continuity limits their applicability. Numerical methods in discontinuous medium mechanics are suitable for simulating large displacements and failures of granular materials, jointed rock masses, and soil-rock masses. Representative methods include the Discrete Element Method (DEM), Discontinuous Deformation Analysis (DDA), and Combined Finite-Discrete Element Method (FDEM). When the discontinuous medium method simulates the contact mechanics behavior between discrete blocks or particles in rocks, its constitutive relation directly controls the macroscopic strength, failure mode and energy dissipation mechanism. Constructing a reasonable discontinuous constitutive model is the core of accurately reproducing the entire process of discontinuous deformation, large displacement and fracture of rocks.

[0003] Although current research methods provide multiple pathways for simulating contact or collision processes in particle groups, they still have many shortcomings. To overcome these shortcomings, further research by those skilled in the art is still needed. Summary of the Invention

[0004] In view of this, the present invention provides a rock failure simulation method based on three-dimensional particle DDA to address one of the deficiencies in the prior art.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: A rock failure simulation method based on three-dimensional particle DDA includes the following steps: The study domain is discretized into a set of multiple rigid spherical particles, and their initial positions, radii, material properties, and control parameters are recorded. A contact determination scheme is adopted to determine the spatial location of the contact point between particles, and a reference distance threshold is introduced to control the contact activation condition; Based on the principle of minimum potential energy, the stiffness sub-matrices corresponding to four types of springs (normal, shear, torsion, and rolling) and the friction sub-matrices corresponding to friction force are generated. Assemble the matrices to form a system of simultaneous equations, and use open-closed iteration to determine the correct contact state; During the process, a set of contact states is first predicted, a system of simultaneous equations is solved, and the contact states after the solution are updated. If the solution is consistent with the predicted contact states, the iteration converges; otherwise, the predicted contact states are updated using the updated contact states, and the iteration is repeated until convergence.

[0006] Optionally, the particle is treated as a rigid body, and its configuration is determined by the position of its center of mass and its radius. The displacement of the particle is approximated by the degree of freedom, which includes translational and rotational degrees of freedom.

[0007] Optionally, the contact point is defined on the line connecting the centers of the two particles, and located at the center of the overlapping area of ​​the two particles. The calculation formula is as follows: ; From particles i and granules j The contact pairs formed have their centroids located at the following positions: and , d 0 for and The distance between them for and The distance between them for and The distance radius between them and and displacement vector and .

[0008] Optionally, the potential energy corresponding to the normal spring sub-matrix is ​​expressed as: ; Among them, particles i and granules j The displacement vectors of the contact pairs are respectively and , It is the stiffness of the normal spring. It is the extended normal vector. It is the normal relative displacement, i.e., the embedding depth. It is the embedding depth at the end of the previous time step.

[0009] Optionally, the potential energy corresponding to the shear spring sub-matrix is ​​expressed as: ; in, It is the stiffness coefficient of the shear spring. It is the shear force at the end of the previous time step. It is a tangential relative displacement.

[0010] Optionally, the potential energy corresponding to the friction submatrix is ​​expressed as: ; Among them, particles i and granules j The displacement vectors of the contact pairs are respectively and , It is a tangential relative displacement. It is Coulomb friction. It is the shear force at the end of the previous time step. It is a tangential projection operator used to project a vector onto a plane perpendicular to the contact normal, thereby extracting the components used to calculate friction, shear force, and torsion.

[0011] Optionally, the potential energy corresponding to the torsional spring sub-matrix is ​​expressed as: ; in, It is torsional stiffness. It is the cumulative torque of the previous step. The torsional component is along the normal direction.

[0012] Optionally, the potential energy corresponding to the submatrix of the rolling spring can be expressed as: ; in It is torsional stiffness. It is the rolling torque of the previous step. It is the difference between the rotation vectors.

[0013] Optionally, the bonding bond is equivalent to an Euler–Bernoulli elastic beam, the stiffness of which is determined by the beam's geometry and material properties, and the failure criterion is determined based on the maximum tensile stress of the beam section and the Mohr-Coulomb shear strength condition.

[0014] Optionally, it also includes updating information such as particle position, velocity, and contact, and proceeding to the next time step of the calculation.

[0015] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a rock failure simulation method based on three-dimensional particle DDA, which has the following beneficial effects: 1. This method effectively overcomes the shortcomings of traditional three-dimensional particle discontinuous deformation analysis and discrete element rock failure simulation methods, where spring stiffness is heavily dependent on particle size, and solves the problem of maintaining consistent macroscopic mechanical response across different particle scales. In this method, macroscopic elastic modulus and Poisson's ratio can be directly used for stiffness estimation, significantly reducing parameter calibration costs.

[0016] 2. The contact point and contact activation distance are defined in the contact determination scheme. This can significantly improve simulation accuracy and numerical stability, specifically as follows: 1) Allowing bonds to form between particles with small initial spacing helps control the condition number of the system; 2) The reference gap mechanism effectively prevents excessive internal forces caused by the initial interlocking of particles, thus enhancing the robustness of the calculation process.

[0017] 3) Eliminate false torques; ensure clear contact force transmission paths and explicit physical meaning.

[0018] 3. The model in this invention can distinguish between particle rolling and torsional behavior, and achieves a high-fidelity characterization of particle contact rotation dynamics by assigning independent spring stiffness to each. This mechanism significantly improves the model's prediction accuracy in simulating rock bending deformation, shear failure, and crack propagation, making the micromechanical response and macroscopic fracture characteristics more reasonable.

[0019] 4. It fills the theoretical gap in the derivation of the contact submatrix system for simulating the failure of discontinuously deformed rocks with three-dimensional particles under complex contact conditions. All submatrices are expressed in vectorized form, laying the foundation for high-efficiency GPU parallel computing.

[0020] 5. This proposal overcomes the limitations of traditional models that oversimplify post-bonding failure processing, such as directly deleting or ignoring friction, and accurately describes the residual friction intensity at the fracture surface and the sliding behavior after shear failure. This proposal can more realistically simulate the entire process of rock from intact to fractured to frictional slip, making the simulated failure mode highly consistent with the rock's physical mechanism, significantly improving the credibility of the fracture process simulation results. Attached Figure Description

[0021] 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 embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0022] Figure 1 A schematic diagram of contact geometry; Figure 2A schematic diagram showing the force and bending moment decomposition between two bonded particles; Figure 3 The diagram shows the failure criteria in PDDA3D, where (a) represents tensile failure and (b) represents shear failure. Figure 4 The following are the deformation verification model diagrams, where (a) is the beam model and (b) is the equivalent PDDA3D model; Figure 5 This is a comparison chart of the results for working condition 1; Figure 6 This is a comparison chart of the results for working condition 2; Figure 7 This is a comparison chart of the results for working condition 3; Figure 8 The diagram shows the destructive verification model, where (a) represents tension and (b) represents shear. Figure 9 Comparison of tensile failure results; Figure 10 This is a comparison chart of shear failure results; Figure 11 The method flowchart provided by the present invention. Detailed Implementation

[0023] 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.

[0024] This invention discloses a method for simulating rock failure based on three-dimensional particle DDA, such as... Figure 11 As shown, it includes the following steps: The study domain is discretized into a set of multiple rigid spherical particles, and their initial positions, radii, material properties, and control parameters are recorded. A contact determination scheme is adopted to determine the spatial location of the contact point between particles, and a reference distance threshold is introduced to control the contact activation condition; Based on the principle of minimum potential energy, the stiffness sub-matrices corresponding to four types of springs (normal, shear, torsion, and rolling) and the friction sub-matrices corresponding to friction force are generated. Assemble the matrices to form a system of simultaneous equations, and use open-closed iteration to determine the correct contact state; During the process, a set of contact states is first predicted, a system of simultaneous equations is solved, and the contact states after the solution are updated. If the solution is consistent with the predicted contact states, the iteration converges; otherwise, the predicted contact states are updated using the updated contact states, and the iteration is repeated until convergence. Update information such as particle position, velocity, and contact, and proceed to the next time step calculation.

[0025] Specifically, this embodiment constructs a bonded particle model within the framework of Three-dimensional Particle-Based Discontinuous Deformation Analysis (PDDA3D). The computational domain is discretized as an assembly of rigid particles, with normal force, tangential force, rolling torque, and torsional torque transmitted between particles through bonding bonds. In this invention, the motion of particles and contact points is calculated using displacement approximation, and the contribution of the stiffness matrix to the overall equations is rigorously derived using the principle of minimum potential energy. This invention equates the bonds between particles to cylindrical beams and derives the initial normal stiffness based on Euler–Bernoulli beam theory. Tangential stiffness Rotational stiffness With torsional stiffness The function relating to macroscopic material parameters, and the determination of bond failure based on tensile strength and the Mohr-Coulomb criterion, are as follows: Step 1: Discretize the study domain into a collection of multiple rigid spherical particles and record their initial positions, radii, material properties, and control parameters.

[0026] Step 2: The particle is treated as a rigid body, and its configuration is completely determined by the position and radius of its center of mass. The displacement of the particle is then approximated by displacement. Step 3: Use a contact determination scheme to determine the spatial location of the contact point between particles, and introduce a reference distance threshold to control the contact activation condition; Step 4: Generate the stiffness sub-matrix and friction force sub-matrix corresponding to the four types of springs (normal, shear, torsion, and rolling) based on the principle of minimum potential energy; Step 5: Establish a model of the bonded particles, treat the bond between particles as an elastic beam, calculate the spring stiffness based on beam theory, and use the tension and shear criteria to determine the bond failure and the contact mode after failure. Step 6: Establish and solve the system equations. The overall stiffness equation is formed by minimizing the total potential energy. After determining the contact state through open-closed iteration, the particle displacement is solved. Step 7: Cyclic simulation and termination determination. Repeat steps S2 to S5 until the preset simulation termination condition is met, and complete the simulation of the entire rock fracture process.

[0027] Furthermore, in step 1, the particle is considered a rigid sphere, and its motion is described by three translational degrees of freedom and three rotational degrees of freedom. The particle displacement is expressed as the product of the particle's displacement transformation matrix and its displacement vector.

[0028] The basic unit of PDDA3D is the rigid particle. In three-dimensional space, each particle is a sphere, and its motion is characterized by six degrees of freedom. For a center coordinate of... = and radius particles Its displacement vector Defined as: (1) in It is the translational displacement vector and Represents the rotation vector (in axis angle form).

[0029] particle i On point p The displacement can be calculated as: (2) (3) in, This is called the displacement transformation matrix, and its function is to approximate the particle displacement using a linear method. generalized displacement vector Mapping to spatial points Local displacement at that location The last three columns encode the displacement components caused by particle rotation. In subsequent derivations, all displacements are calculated at the contact point. Therefore, for simplicity, Abbreviation .

[0030] It should be noted that updating the particle boundary point coordinates using formula (2) may introduce so-called "free expansion error" under finite rotation conditions. However, in this embodiment, the particle is treated as a rigid body, and its configuration is entirely determined by the position of its center of mass and its radius, so there is no need to explicitly update the position of the boundary points. Furthermore, the specific coordinates of the contact points are not stored in the current model framework (as described later), so no positional error will accumulate at the contact level.

[0031] Furthermore, in step 2, the contact point is defined on the line connecting the centers of the two particles and is located at the center of the overlapping area of ​​the two particles.

[0032] Considering particles i and granules j The contact pairs formed have their centroids located at the following positions: and ,radius and and displacement vector and For ease of explanation, Figure 1A simplified two-dimensional diagram is given in the figure. From the particles i Pointing to particles j The unit normal vector is defined as: (4) Assuming contact point Located in the connecting particles i and j On the line of the center of mass, and located at the center of the overlapping region of the two particles, such as Figure 1 As shown. Contact point The formula for calculating the position is: (5) In particles and For reference, contact point The relative displacement is defined as: (6) Furthermore, in step 3, as Figure 2 As shown, this relative displacement can be decomposed into two independent components: normal relative displacement and tangential relative displacement. In the normal direction, it is described in full form. Here, the normal relative displacement (embedding depth) It can be calculated as follows: (7) in (8) Here, d gap This is a reference gap used to define the contact activation distance. It is the extended normal vector, used only for formula derivation. This is the embedding depth at the end of the previous time step. Except for the normal direction, the tangential and rotational directions use incremental formulas. In subsequent analysis, all relative displacements mentioned are in incremental form. The tangential relative displacement is calculated as follows:

[0033] Among them, the projection operator Project any vector onto a surface perpendicular to the unit normal vector. n The operator is symmetric in the plane. ) and idempotency ( ), and They are particles and granules The displacement transformation matrix, It is a 3×3 identity matrix.

[0034] Particles and granules The relative rotation between them is defined by the difference between their rotation vectors:

[0035] in They are particles and granules The rotation vector.

[0036] Relative rotation can be decomposed into torsional and rolling components, such as Figure 2 As shown. The torsional component is along the normal direction, and its expression is: (11) The scroll component is: (12) Since the direction of the torsional component is consistent with the normal to the contact surface, only its scalar magnitude needs to be stored. This scalar is defined as follows: .

[0037] Furthermore, in step 4, the contact-related sub-matrices are derived based on the principle of minimum potential energy, including the normal spring sub-matrices, shear spring sub-matrices, friction sub-matrices, torsional spring sub-matrices, and rolling spring sub-matrices.

[0038] First, we introduce the mechanical model of normal contact, as it determines whether the contact is activated. When a pair of contacting particles are in an bonded state, or even if not bonded... When contact is activated, a normal spring is introduced, the potential energy of which is expressed as follows: (13) in It is the stiffness of the normal spring.

[0039] Expanding equation (13) yields: (14) By minimizing about and Given: (15) When using the incremental approach, the potential energy corresponding to the shear spring can be expressed as: (16) in It is the stiffness coefficient of the shear spring. This is the shear force at the end of the previous time step. At the end of each time step, if the shear force does not exceed the frictional force, the shear force should be updated as follows: (17) here It is the shear force at the end of the current time step, which will be used as the shear force in the next time step. The initial value. If the shear force exceeds the frictional force, It will be updated to Coulomb friction.

[0040] Equation Substitute into equation (16), and use Expanding on its properties, we get: (18) By minimizing At once and In other words: (19) When the shear force exceeds the Coulomb friction force and the contact remains active, the contact state will change to sliding. At this point, the shear spring will be replaced by a pair of frictional forces, the magnitude of which is calculated as follows: (20) in Let be the friction angle. The embedding depth at the end of the previous step, as defined by equation (7). This represents the normal contact force at the end of the previous time step. The embedment depth from the previous time step is used here because the embedment depth in the current time step is still unknown. The direction of the frictional force depends on the contact point. The relative motion direction. Since the relative motion direction at the contact point is also unknown, the direction of the frictional force is determined by the shear force of the previous step. The direction is determined by.

[0041] In this case, the shear force is expressed as: (twenty one) in It is the unit shear vector, and the calculation formula is: (twenty two) The potential energy of friction is: (twenty three) By respectively and minimize We can obtain: (twenty four) The potential energy associated with a torsion spring is expressed as: (25) in It is torsional stiffness. This is the cumulative torque magnitude from the previous time step. At the end of the time step, the torque... Update as follows: (26) Using equation (25) Expanding equation (10), we get: (27) By respectively and minimize This allows us to obtain the contribution of torque to the overall stiffness matrix and force vector. Specifically, the terms related to the rotational degrees of freedom are as follows: (28) The superscript here "" indicates the stiffness matrix and force vectors that are only related to the rotational degree of freedom. For example, It is a 3×3 submatrix, corresponding to the global stiffness matrix. The 6i-2 arrive 6i row and number 6j-2t arrive 6j List; It is a 3×1 subvector, corresponding to the local force vector. of 6j-2t arrive 6j OK.

[0042] According to formula (12), the potential energy of the rolling spring is expressed as: (29) in It is torsional stiffness. This is the rolling torque of the previous time step. At the end of the time step, the rolling torque... Update as follows: (30) It should be noted that in equation (26) It is a scalar, and here it is... It is a vector.

[0043] Expanding on (29) further, we get: (31) Through the Regarding nodal displacement variables and By minimizing this, we can obtain the submatrix of the rolling spring: (32) Furthermore, in step 5, the bonding bond is equivalent to an Euler–Bernoulli elastic beam, the stiffness of which is determined by the beam's geometry and material properties, and the failure criterion is determined based on the maximum tensile stress of the beam section and the Mohr-Coulomb shear strength condition.

[0044] Connection radius is and The bond between the two particles is idealized to have a radius of... , length is A cylindrical beam. Typically, the bond radius is determined by the smaller or average of the two particle radii: (33) in The shape parameter (typically set to 1.0 or less) controls the size of the cemented bond. This parameter directly affects the bond's cross-sectional area and moment of inertia, thus influencing the overall stiffness and strength. For the rock material considered in this embodiment, It is usually set to 1.0.

[0045] Given the material parameters of the beam (elastic modulus) Compared to Poisson ) and geometric parameters ( and Under the premise that the spring stiffness of the four types can be calculated as follows: Normal spring stiffness: (34) in Let be the cross-sectional area of ​​the bond.

[0046] Shear spring stiffness: (35) in It is the shear modulus.

[0047] Torsional spring stiffness: (36) in It is the polar moment of inertia.

[0048] Rolling spring stiffness: (37) in It is the moment of inertia of the cross section.

[0049] In BPM, microcrack initiation is captured through progressive bond failure, and these microcracks eventually converge to form macroscopic cracks. In this embodiment, the bonds may fail under tensile or shear stress. For a typical bonded connection, the maximum normal stress and shear force are calculated as follows: (38) (39) in This is the torque contribution factor. Note that... This indicates that the contact is under tension. In the current work, Set it to 0 to obtain a more accurate ratio of compressive strength to tensile strength. For example... Figure 3 As shown in (a), when the following conditions are met Tensile cracking occurs at times, among which It is the tensile strength of the bond. Once tensile failure occurs, the contact force and torque drop to 0, and the contact state becomes open.

[0050] Shear failure adopts the Mohr-Coulomb type failure criterion. (40) in It is the cohesive force of the bond. It is the friction angle. After shear failure, the shear force becomes sliding friction, such as... Figure 3 As shown in (b), where This represents the cumulative shear displacement. Simultaneously, the torque also degenerates to zero.

[0051] In this embodiment, tensile failure is examined first. Once the key breaks, the tensile strength... and cohesion Reduce to zero in the following steps.

[0052] Furthermore, in step 6, the equations are established by minimizing the total potential energy of the system, which includes the potential energy contributions from the particles and the contact interface.

[0053] The derivation of the global governing equations in PDDA3D is the minimization of the system's total potential energy. This total potential energy consists of the potential energy contributed by each particle and each contact interface in the system, specifically including the contact between particles and the contact between particles and boundary surfaces. The potential energy of a single particle includes the energy generated by inertial forces, concentrated loads, volume forces (such as gravity), displacement, velocity constraints, damping, etc. For contact interfaces, the potential energy is contributed by normal springs, shear springs, rolling springs, torsional springs, and frictional forces.

[0054] Consider a n The global equilibrium equation of a system composed of particles can be expressed as: (41) in K It is a size of 6 n ×6 n The global stiffness matrix, D It is 6 n A displacement vector of ×1, F It is 6 n A force vector of ×1.

[0055] Before establishing the simultaneous equations, the contact states between each contact pair must be pre-assumed. The process of constructing and solving the equations is achieved by combining the assumed contact states with the mechanical model described in Section 4. The solution only has physical meaning when the calculated contact states are consistent with the initial assumptions. In DDA, the process of identifying the correct contact states is called "open-closed iteration".

[0056] To verify the deformation accuracy of the proposed BPM, an end-loaded cantilever beam model was used. Figure 4 (a)). The beam is 1 m long, with a circular cross-section radius of 0.025 m, and an elastic modulus of... E =68.95 GPa, shear modulus G =25.92 GPa. Three loading conditions were set (the specific application magnitude for each condition is shown in Table 1): Condition 1 along y Axial tension is applied to the shaft; operating condition 2. y Applying torque to the shaft T Condition 3 applies a perpendicular force to the free end. y Concentrated force on the shaft F The calculation uses a dynamic damping coefficient (same as traditional DDA) to dissipate energy, and the equilibrium displacement is used as the simulation result.

[0057] The simulation uses particles of equal diameter arranged in sequence. Figure 4 (b) Discretize the beam into 21 particles with a radius of 0.025 m. The PDDA3D parameters are set as follows: and E Equal tangential stiffness according to G / E =0.376. Since only elastic deformation is being verified, the strength parameter... , , All are set to sufficiently high values ​​to prevent corruption.

[0058] Table 1. Magnitude of Force (Torque) Applied Under Various Working Conditions

[0059] The analytical solutions for operating conditions 1, 2, and 3 are given by formulas (42), (43), and (44), respectively. The comparison results are as follows: Figure 5 , Figure 6 as well as Figure 7 As shown, the analytical solutions and numerical simulation results under various working conditions agree well, verifying the accuracy of PDDA3D in deformation analysis.

[0060] (42)

[0061] In the formula This refers to the axial tensile length. P It is an axial force. y The distance from the fixed end is the length. m² is the cross-sectional area of ​​the beam. E It is the elastic modulus.

[0062] (43)

[0063] In the formula To bypass y The angle of twist of the shaft, T For torque, G Shear modulus It is the polar moment of inertia.

[0064] (44)

[0065] In the formula For deflection, F To apply perpendicular to the free end y Concentrated force on the shaft, L Let be the length of the beam. / 4 represents the moment of inertia.

[0066] This embodiment verifies the correctness of the tensile and shear failure simulations through bond failure testing. For example... Figure 8 As shown, consider a particle i and granules j The contact pair consists of particles with a radius of 0.5 meters. The two particles are connected by a bond, idealized as a cylindrical beam with a cross-sectional area of... A b =0.785 m², length L b =1 m.

[0067] In tensile failure tests, particles i Fixed, while particles j A constant vertical velocity is applied. v 1 = 0.2 m / s. The macroscopic parameter is taken as... GPa and MPa. According to formulas (34) and (38), we get The particle size distribution was N / m, and the analytical tensile peak displacement was 0.001 m. During the simulation, particle size was recorded. i and j normal stress between (Normal force divided by cross-sectional area) A b ). Figure 9 The numerical solution was compared with the corresponding analytical solution, and the results show that they agree well.

[0068] In shear failure tests, such as Figure 8 As shown in (b), the particles i Fix, while also affecting the particles j Apply normal stress =1MPa. After the system reaches a steady state, the particles... j Apply constant horizontal speed v 2 = 0.2 m / s. The calculation parameters are set as follows: GPa, =0.3, =2 MPa, =30°. Calculated according to formula (35) k s =2.36×10 8 N / m. Based on From formula (40), we can obtain that the shear displacement corresponding to the peak shear stress is 0.0086 m, the peak shear stress is 2.58 MPa, and the residual stress after bond failure is 0.58 MPa. Figure 10 As shown, the numerical calculation results agree well with the analytical solution throughout the loading process.

[0069] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since it corresponds to the method disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.

[0070] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for simulating rock failure based on three-dimensional particle DDA, characterized in that, Includes the following steps: The study domain is discretized into a set of multiple rigid spherical particles, and their initial positions, radii, material properties, and control parameters are recorded. A contact determination scheme is adopted to determine the spatial location of the contact point between particles, and a reference distance threshold is introduced to control the contact activation condition; Based on the principle of minimum potential energy, the stiffness sub-matrices corresponding to four types of springs (normal, shear, torsion, and rolling) and the friction sub-matrices corresponding to friction force are generated. Assemble the matrices to form a system of simultaneous equations, and use open-closed iteration to determine the correct contact state; During the process, a set of contact states is first predicted, a system of simultaneous equations is solved, and the contact states after the solution are updated. If the solution is consistent with the predicted contact states, the iteration converges; otherwise, the predicted contact states are updated using the updated contact states, and the iteration is repeated until convergence.

2. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The particles are treated as rigid bodies, and their configuration is determined by the position of their center of mass and their radius. The displacement of the particles is approximated by their displacement, which is characterized by their degrees of freedom, including translational and rotational degrees of freedom.

3. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The contact point is defined on the line connecting the centers of the two particles, and is located at the center of the overlapping area of ​​the two particles. The calculation formula is as follows: ； From particles i and granules j The contact pairs formed have their centroids located at the following positions: and , d 0 for and The distance between them for and The distance between them for and The distance radius between them and and displacement vector and .

4. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The potential energy corresponding to the normal spring sub-matrix is ​​expressed as: ； Among them, particles i and granules j The displacement vectors of the contact pairs are respectively and , It is the stiffness of the normal spring. It is the extended normal vector. It is the normal relative displacement, i.e., the embedding depth. It is the embedding depth at the end of the previous time step.

5. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The potential energy corresponding to the shear spring sub-matrix is ​​expressed as: ； in, It is the stiffness coefficient of the shear spring. It is the shear force at the end of the previous time step. It is a tangential relative displacement.

6. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The potential energy corresponding to the friction submatrix is ​​expressed as: ； Among them, particles i and granules j The displacement vectors of the contact pairs are respectively and , It is a tangential relative displacement. It is Coulomb friction. It is the shear force at the end of the previous time step. It is a tangential projection operator.

7. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The potential energy corresponding to the torsional spring sub-matrix is ​​expressed as: ； in, It is torsional stiffness. It is the cumulative torque of the previous step. The torsional component is along the normal direction.

8. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The potential energy corresponding to the submatrix of the rolling spring is expressed as: ； in It is torsional stiffness. It is the rolling torque of the previous step. It is the difference between the rotation vectors.

9. The rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, The bond is equivalent to an Euler–Bernoulli elastic beam, whose stiffness is determined by the beam's geometry and material properties. The failure criterion is based on the maximum tensile stress in the beam section and the Mohr-Coulomb shear strength condition.

10. A rock failure simulation method based on three-dimensional particle DDA according to claim 1, characterized in that, It also includes updating information such as particle position, velocity, and contact, and proceeding to the next time step of the calculation.

Images