A method and system for predicting rock failure process based on discontinuous deformation analysis

Abstract

The present invention belongs to the field of rock and soil mechanics technology, and specifically discloses a method and system for predicting the rock failure process based on discontinuous deformation analysis. It includes: discretizing and modeling the rock in different regions, applying parallel spring units between adjacent discrete blocks based on elastic theory, and establishing a nonlinear failure criterion based on fracture mechanics theory; inputting a time step, calculating the normal phase and tangential deformation of the bonded spring, and updating the stress-strain information of the block in real time; calculating the damage coefficient, judging the degree of damage to the bonded spring, removing the damaged bonded spring unit, and updating the stiffness of the undamaged spring; and visualizing the rock failure process. The present invention takes into account the stress softening behavior in the rock fracture process, overcomes the problem of insufficient accuracy of the existing discontinuous deformation analysis method (DDA) in simulating the nonlinear stress softening behavior of rock, and can more accurately predict the failure process and stress evolution process of rock under tension and shear, and predict the strength and fracture characteristics of rock.

Description

Technical Field

[0001] The present invention belongs to the technical field of rock and soil mechanics, and more specifically, relates to a method and system for predicting rock failure process based on discontinuous deformation analysis. Background Art

[0002] As a common geological material, rock has very complex deformation and failure behaviors and mechanical responses, which directly affect the stability and safety of mines, tunnels and other projects. Therefore, predicting and analyzing the strength and fracture characteristics of rock is of great significance for the stability evaluation of rock engineering. The macroscopic mechanical properties of rock are closely related to its mesostructure. Studying the evolution law of rock mesostructure can reveal the nature of its macroscopic strength and deformation. In recent years, numerical simulation has become an important means to solve this complex problem.

[0003] Discrete element method (DEM) and discontinuous deformation analysis (DDA) are commonly used methods for analyzing the mechanical properties of rock materials. Discontinuous deformation analysis (DDA) has obvious advantages in calculating long time steps and quasi-static problems because it can directly process irregular blocks and consider block deformation, use larger time steps, and ensure unconditional stability of calculations. Therefore, it has become an important method for analyzing rock failure laws and mechanical properties. When using DDA to simulate rock materials, rocks are usually discretized into deformable blocks, and spring units are applied at the contact positions of the blocks. The springs not only transfer stress between blocks, but also simulate the initiation and expansion of cracks in rocks through their own destruction.

[0004] The original DDA adopts the brittle failure criterion and does not consider the stress softening behavior of the material. On this basis, some scholars have made improvements. For example, the paper: Xia M, Chen G, Yu P, et al. Improvement of DDA with a New Unified Tensile Fracture Model for Rock Fragmentation and its Application on Dynamic Seismic Landslides [J]. Rock Mechanics and Rock Engineering. 2021, 54 (3): 1055-75. A unified tensile constitutive model is proposed to improve the accuracy of fracture simulation of brittle materials under tension; the paper: Hu C, Gong S, Zhuang D, et al. Improved discontinuous deformation analysis method for modeling tensile cracking in quasi-brittle materials [J]. Theoretical and Applied Fracture Mechanics. 2023, 124 (2023): 103748. A multi-segment linear tensile failure model is introduced, and two linear models are used to approximate the stress-strain relationship after failure. However, the existing improved DDA model generally adopts a linear failure criterion, which has the problem of insufficient accuracy in describing the true stress-strain characteristics of rock. Experiments show that the softening process of rock is nonlinear. On the other hand, the existing DDA model generally uses two sets of contact spring units to transfer stress and strain between blocks, which cannot withstand tensile stress and is difficult to accurately calculate tensile strain.

[0005] In summary, establishing a numerical method that can more accurately simulate the rock failure process and analyze the stress softening behavior of rocks is of great significance for predicting the strength and fracture characteristics of rocks and evaluating the stability of rock engineering. Summary of the invention

[0006] In view of the above defects or improvement needs of the prior art, the present invention provides a method and system for predicting the rock failure process based on discontinuous deformation analysis. On the one hand, it utilizes the discontinuous deformation analysis method, which can simulate irregular blocks and has higher calculation efficiency. On the other hand, it introduces the exponential tensile failure constitutive relationship, which can more accurately represent the stress softening behavior of rocks; on the other hand, three parallel bonded spring units are established between blocks to transfer stress and strain between blocks, which can more accurately capture the properties of joints. This method can capture the crack initiation position, crack initiation process, stress-strain curve and stress evolution characteristics of rocks, with higher accuracy, and the key strength parameters can be determined by experiments, providing a more accurate method for predicting and evaluating the strength and fracture characteristics of rocks.

[0007] To achieve the above object, according to one aspect of the present invention, a method for predicting rock failure process based on discontinuous deformation analysis is proposed, comprising the following steps:

[0008] S100: Obtain the geometric and physical parameters of the rock, conduct regional discretization modeling based on the geometric characteristics and microstructural development characteristics of the rock, and impose boundary conditions on the model;

[0009] S200: Based on elastic theory, a bonding spring unit is applied between adjacent discrete blocks to transmit stress and strain information between the blocks;

[0010] S300: Based on the fracture mechanics theory, a nonlinear failure criterion is established and embedded in the DDA algorithm to determine the degree of damage of the bonded spring;

[0011] S400: Enter the time step to start the calculation, calculate the normal phase and tangential deformation of the bonded spring unit, update the stress and strain information of the block in real time, then calculate the damage coefficient according to the nonlinear failure criterion described in step S300, and judge the damage degree of the bonded spring unit, remove the bonded spring unit whose damage degree reaches the failure condition, and update the undamaged spring stiffness, and update the block deformation according to the spring stiffness;

[0012] S500: According to the calculation result of step S400, the rock crack initiation position is captured, and the crack initiation process, stress-strain curve and stress evolution process are generated.

[0013] As a further preferred embodiment, in step S100, the regional discretization modeling includes:

[0014] The degree of discretization is determined according to the degree of microstructure development and geometric characteristics. For areas with relatively developed microstructures or complex geometric shapes, small-sized block discretization is used, and large-sized block discretization is used in other areas to establish a discretized numerical model.

[0015] As further preferred, in step S100, the physical parameters of the rock include: stiffness, density, Poisson's ratio, tensile strength, cohesive strength, and friction angle.

[0016] As further preferred, in step S200, three groups of parallel bonded spring units are included, and each group of bonded spring units includes two bonded springs, a normal bonded spring and a tangential bonded spring.

[0017] As a further preference, in step S200, the material parameters of the bonded spring unit need to be input, and the material parameters include bonded tensile strength, bonded cohesive strength, bonded tensile fracture energy, bonded shear fracture energy, and bonded spring normal stiffness, and the above material parameters are calculated and determined based on the physical parameters of the rock.

[0018] As a further preference, in step S300, a nonlinear failure criterion is used to judge the damage and fracture process of the reaction rock. The nonlinear failure criterion is manifested as an exponential tensile failure constitutive relationship, which includes two stages, pre-peak and post-peak. In the pre-peak stage, the cementation is not damaged, and it mainly manifests as elastic deformation, and the stress-strain relationship is linear; in the post-peak stage, the cementation begins to be damaged, the strength gradually decreases, and the stress-strain relationship is exponential.

[0019] As further preferred, step S400 specifically includes the following steps:

[0020] S410: Input the time step and start calculation;

[0021] S420: Calculate bonded spring deformation: Calculate block deformation using the DDA control equation, and then calculate the bonded spring normal based on the block deformation. , Tangential deformation The bonded spring phase , Tangential deformation The computational models include:

[0022] ,

[0023] in, and are the normal and tangential lengths at the beginning of the current time step, and are the normal and tangent unit vectors, and Spring action points i and j displacement;

[0024] S430: Calculate the damage coefficient of the bonded spring and determine the degree of damage to the bonded spring: Calculate the damage coefficient based on the deformation of the bonded spring combined with the nonlinear failure criterion , using the damage factor To reflect the degree of bond damage, the damage state of the bonded spring is determined according to the damage coefficient: =0, in the elastic deformation stage, >0, in the stress softening stage, >0.999, the bonding is considered to be damaged and the spring unit is removed. The calculation model of the damage coefficient includes:

[0025] ,

[0026] In the formula, is the deformation of the bonding unit, r = n , represents the normal deformation, r = s , represents the tangential deformation, The calculated value of is equal to the average value of the three groups of spring deformations, and is the elastic critical displacement and damage parameter;

[0027] S440: Calculate and update the bond spring stiffness based on the damage coefficient Calculate and update spring normal and tangential stiffness ; Updated spring normal, tangential stiffness include:

[0028] ,

[0029] in, is the input bond stiffness;

[0030] S450: Reconstruct the DDA control equation and calculate and update the block displacement.

[0031] As further preferred, in step S450, the reconstructing the DDA control equation includes:

[0032] S451: Calculate the spring coefficient matrix and force matrix based on the minimum potential energy theorem;

[0033] S452: Reconstruct the DDA control equation according to the bond spring stiffness matrix and the force matrix, and calculate and update the block displacement. The calculation model of the DDA control equation is:

[0034]

[0035] in, is the mass matrix; is the damping matrix; is the elasticity matrix; is the contact spring coefficient matrix; is the bonding spring coefficient matrix; is the particle contact force and load matrix; is the bonding spring force matrix; and is the velocity and acceleration matrix, and D is the displacement matrix.

[0036] In the above steps, the spring coefficient matrix and moment matrix are calculated as:

[0037]

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] Where T is the displacement transformation matrix, , , , is the spring coefficient matrix; is the total energy generated by the cemented spring; , Spring action points , The displacement transformation matrix at ; , and , They are the degrees of freedom of the two blocks acted upon by springs; and are the normal and tangential stiffness of the cemented spring, respectively; and are the normal and tangential unit vectors of the bonded spring; Fi, Fj are the spring force matrices; and are the normal and tangential lengths at the beginning of the current time step; the subscripts i and j are the block numbers.

[0044] As a further preferred embodiment, step S500 specifically includes the following steps: converting and storing the calculated geometric structure and corresponding physical quantity data according to the standard data architecture of the VTK file, importing and parsing the generated VTK file through the ParaView platform, and performing visual analysis on the crack initiation position, crack initiation process, stress-strain curve and stress evolution process.

[0045] According to another aspect of the present invention, a system for predicting rock failure process based on discontinuous deformation analysis is provided, comprising:

[0046] The first main control module is used to obtain the geometric and physical parameters of the rock, conduct regional discretization modeling based on the geometric characteristics and microstructural development characteristics of the rock, and impose boundary conditions on the model;

[0047] The second main control module is used to apply a bonding spring unit between adjacent discrete blocks for transmitting stress and strain information between the blocks based on elasticity theory;

[0048] The third main control module is used to establish a nonlinear failure criterion based on fracture mechanics theory and embed it into the DDA algorithm to determine the degree of damage to the bonded spring;

[0049] The fourth main control module inputs a time step to start calculation, calculates the normal phase and tangential deformation of the bonded spring unit, updates the stress and strain information of the block in real time, and then calculates the damage coefficient according to the nonlinear failure criterion described in S300, and determines the damage degree of the bonded spring unit, removes the bonded spring unit whose damage degree reaches the failure condition, and updates the undamaged spring stiffness, and updates the block deformation according to the spring stiffness;

[0050] The fifth main control module is used to capture the rock crack initiation position according to the calculation results of the fourth main control module, and generate the crack initiation process, stress-strain curve and stress evolution process.

[0051] In general, the above technical solution conceived by the present invention has the following technical advantages compared with the prior art:

[0052] 1. The present invention obtains the geometric parameters, physical parameters, geometric characteristics and microstructural development characteristics of rocks to perform regional discretization modeling, and applies boundary conditions, applies parallel spring units between blocks to transmit stress and strain information, establishes nonlinear failure criteria, and is used to judge the degree of damage to the cemented spring. During calculation, input the time step to start the program, calculate and update the deformation of the cemented spring and the stress and strain information of the block; calculate the damage coefficient, judge the degree of damage to the cemented spring, remove the damaged cemented spring unit, and update the stiffness of the undamaged spring; extract the simulation results for visual analysis to obtain the crack initiation position, crack initiation process, stress-strain curve and stress evolution process, which overcomes the problem of insufficient accuracy of the existing discontinuous deformation analysis method in simulating the nonlinear stress softening behavior of rocks, and can more accurately predict the failure process and stress evolution process of rocks under tension and shear, and predict the strength and fracture characteristics of rocks.

[0053] 2. The present invention adopts a discontinuous deformation analysis method, which can simulate irregular blocks such as rocks, and can adopt a long time step with high calculation efficiency.

[0054] 3. The present invention introduces an exponential tensile failure model into the discontinuous deformation analysis method, which can effectively overcome the problem that the existing discontinuous deformation method cannot accurately simulate the nonlinear softening behavior of rock, and the simulation results are more in line with the actual situation.

[0055] 4. The present invention uses three parallel spring units between blocks to transfer stress and strain, which can overcome the problem that the contact spring in the existing discontinuous deformation analysis method cannot withstand tensile stress and accurately calculate tensile strain. The method proposed by the present invention can accurately calculate tensile stress and tensile strain, thereby accurately capturing joints.

[0056] 5. The present invention can more accurately capture the crack initiation position, crack initiation process, stress-strain curve and stress evolution of rocks. The simulated rock failure process is more consistent with the experimental results and can effectively predict and evaluate the strength and fracture characteristics of rocks.

[0057] 6. The present invention is an organic combination of the discontinuous deformation method and the post-processing software ParaView, which can visualize the crack development process and the stress and strain evolution process, and has a good visualization effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 A flow chart of a method for predicting rock failure process based on discontinuous deformation analysis in an embodiment of the present invention;

[0059] Figure 2 A schematic diagram of a method for predicting rock failure process based on discontinuous deformation analysis in an embodiment of the present invention;

[0060] Figure 3 is the nonlinear failure criterion involved in the embodiment of the present invention;

[0061] Figure 4 (a) is a schematic diagram of the initial conditions of the sample in the embodiment of the present invention. Figure 4 (b) is a schematic diagram of the initial discretization model of the sample in the embodiment of the present invention;

[0062] Figure 5 (a) is the final failure form of the test result. Figure 5 (b) is the final destruction form simulated by the embodiment of the present invention;

[0063] Figure 6 is a force-displacement curve in one embodiment of the present invention;

[0064] Figure 7 (a) is the evolution diagram of the principal stress field in the initial state in the embodiment of the present invention. Figure 7 (b) is the evolution diagram of the principal stress field in the development state. Figure 7(c) is the evolution diagram of the principal stress field in the final state;

[0065] Figure 8 (a) is the displacement field evolution diagram of the initial state in the embodiment of the present invention, Figure 8 (b) is the displacement field evolution diagram of the development state. Figure 8 (c) in the figure is the displacement field evolution diagram of the final state. DETAILED DESCRIPTION

[0066] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0067] like Figure 1 As shown, a method for predicting rock failure process based on discontinuous deformation analysis provided by an embodiment of the present invention includes the following steps:

[0068] S100: Based on field investigation and indoor tests, the geometric and physical parameters of the rock are obtained and the parameters are input; the rock is discretized and modeled in different regions according to its geometric characteristics and microstructural development characteristics; according to the actual working conditions, boundary conditions are imposed on the discretized model, including load conditions and detection point locations.

[0069] In this step, the numerical model discretization process is: determine the degree of discretization according to the degree of microstructural development and geometric characteristics, use smaller blocks for discretization in areas with more developed microstructures or complex geometric shapes, and use larger blocks for discretization in other areas, and then establish a discretized numerical model. The physical parameters of the rock include stiffness, density, Poisson's ratio, tensile strength, cohesive strength, and friction angle, and the above physical parameters are obtained based on indoor tests.

[0070] S200: Based on elasticity theory, three sets of parallel spring units are applied between adjacent discrete blocks to transmit stress and strain information between blocks.

[0071] The initial spring unit includes three groups of parallel cemented spring units, each group of spring units includes two cemented springs in the normal direction and the tangential direction. The initial spring unit applied needs to input its material parameters, including cemented tensile strength, cemented cohesive strength, cemented tensile fracture energy, cemented shear fracture energy, and cemented spring normal stiffness. The above material parameters are calculated and determined according to the physical parameters of the rock.

[0072] S300: Based on the fracture mechanics theory, a nonlinear failure criterion is established to judge the degree of damage of the bonding spring, and is embedded in the DDA algorithm to judge the degree of damage of the bonding spring, thereby simulating the damage and nonlinear stress softening of rocks. The nonlinear failure criterion used is used to judge the damage and fracture process of the reaction rock. The application of this model overcomes the shortcomings of the existing discontinuous deformation analysis method that cannot describe the nonlinear softening behavior of the material. The nonlinear failure criterion used is characterized by: the model is an exponential tensile failure constitutive relationship, which includes two stages, pre-peak and post-peak. In the pre-peak stage, the bonding does not fail, and it mainly manifests as elastic deformation, and the stress-strain relationship is linear; in the post-peak stage, the bonding begins to be damaged, the strength gradually decreases, and the stress-strain relationship is exponential.

[0073] S400: Input the time step, start the program, calculate the normal phase and tangential deformation of the cemented spring, and update the stress and strain information of the block in real time; calculate the damage coefficient according to the nonlinear failure criterion described in step S300, and judge the damage degree of the cemented spring, remove the cemented spring unit whose damage degree reaches the failure condition, and update the undamaged spring stiffness. This step is used to simulate the rock failure process.

[0074] More specifically, in step S400, the steps of simulating the rock destruction process are:

[0075] S410: Input the time step and start calculation.

[0076] S420: Calculate bonded spring deformation: Calculate block deformation using the DDA control equation, and then calculate the bonded spring normal based on the block deformation. , Tangential deformation The bonded spring phase , Tangential deformation The computational models include:

[0077]

[0078] in, and are the normal and tangential lengths at the beginning of the current time step, and are the normal and tangent unit vectors, and Spring action points i and j displacement.

[0079] S430: Calculate the damage coefficient of the bonded spring and determine the degree of damage to the bonded spring: Calculate the damage coefficient based on the deformation of the bonded spring , using the damage factor To reflect the degree of bond damage. The damage coefficient is calculated as:

[0080]

[0081] In the formula, is the deformation of the bonding unit, r = n , represents the normal deformation, r = s , which represents the tangential deformation, is equal to the average value of the deformation of the three groups of springs. The bonded spring normal phase and tangential deformation are referred to in the following text. r The meaning of is deduced from this; and is the elastic critical displacement and damage parameter.

[0082] Determine the damage state of the bonded spring based on the damage coefficient: =0, in the elastic deformation stage, >0, in the stress softening stage, >0.999, the bond is considered to be broken and the spring element is removed.

[0083] S440: Calculate and update bond spring stiffness: based on damage factor Calculate and update spring normal and tangential stiffness , calculated as:

[0084]

[0085] in, is the input bond stiffness.

[0086] S450: Reconstruct the DDA control equation.

[0087] In step S450, the method for reconstructing the DDA control equation is:

[0088] S451: Calculate the spring coefficient matrix and force matrix based on the minimum potential energy theorem, and the calculation is:

[0089]

[0090]

[0091]

[0092]

[0093]

[0094]

[0095] in, T is the displacement transformation matrix, , , , is the spring coefficient matrix; is the total energy generated by the cemented spring; , Spring action points , The displacement transformation matrix at ; , and , They are the degrees of freedom of the two blocks acted upon by springs; and are the normal and tangential stiffness of the cemented spring, respectively; and are the normal and tangential unit vectors of the bonded spring; Fi, Fj are the spring force matrices; and are the normal and tangential lengths at the beginning of the current time step; the subscripts i and j are the block numbers.

[0096] S452: Substitute the bond spring stiffness matrix and force matrix, reconstruct the DDA control equation, and calculate and update the block displacement. The control equation is:

[0097]

[0098] in, is the mass matrix; is the damping matrix; is the elasticity matrix; is the contact spring coefficient matrix; is the bonding spring coefficient matrix; is the particle contact force and load matrix; is the bonding spring force matrix; and is the velocity and acceleration matrix; D is the displacement matrix.

[0099] S500: Analyze the simulation results, call the dedicated module, convert and store the calculation results according to the standard data structure of the Visualization Toolkit (VTK) file, and use the ParaView software to parse the exported VTK file to visualize the crack initiation position, crack initiation process, stress-strain curve and stress evolution process.

[0100] In this step, the calculated geometric structure and corresponding physical quantity data are converted and stored according to the standard data architecture of the VTK file by calling a dedicated data conversion module. The generated VTK file is imported and parsed through the Paraview platform to visualize the crack initiation position, crack initiation process, stress-strain curve and stress evolution process.

[0101] In a specific embodiment of the present invention, the flow chart of steps S1-S5 is as follows: Figure 1 As shown in the figure, the principle of its numerical algorithm Figure 2 shown.

[0102] The nonlinear failure criterion in step S300 is that the deformation before the peak point is elastic, damage begins after the peak point, and the stress-strain relationship is exponential.

[0103] In this embodiment, a tensile test is performed on a rock double-notched rectangular plate. The test chamber reflects a standard case of material crack initiation and evolution. The specific experimental steps are as follows:

[0104] S11: In this embodiment, the sample is 60 mm long, 120 mm high, and 10 mm thick. Two cuts of 10 mm long and 2 mm wide are preset in the middle of the two edges of the rectangular plate, and the distance between the two cuts is 5 mm. The geometric figure is drawn using CAD software and imported into GID software.

[0105] S12: Discretize the rock in different regions in the GID software according to its geometric characteristics and microstructural development characteristics, and then export the discretized model as a TXT file. In this embodiment, 1 mm blocks are used to refine the middle part of the rectangular plate, and 20 mm blocks are used for the rest.

[0106] S13: According to the experiment, the physical parameters of the rock double-cut rectangular plate are obtained and input into a TXT file. The physical parameters of this embodiment are shown in Table 1.

[0107] Table 1 Physical parameters of rock double-cut rectangular plate

[0108]

[0109] S14: Write the boundary conditions into the TXT file. In this case, the bottom of the rectangular plate is fixed, and a uniform load with a speed of 0.3 mm / s is applied to the top. Two monitoring points, m1 and m2, are set 15 mm above the two cuts to record the deformation of the specimen during the test.

[0110] The rock double-cut rectangular plate tensile model established in steps S11-S14 is as follows: Figure 4 shown.

[0111] The steps S2-S5 in the present invention are based on the discontinuous deformation theory and are written in a self-compiled M-DDA (improved discontinuous deformation method) program using C++ language.

[0112] Use the self-written M-DDA program to read the TXT file created by S11-S14 and generate the M-DDA model.

[0113] S2: Based on elasticity theory, three groups of parallel spring units are applied between adjacent blocks of the M-DDA model. Each group of spring units includes two normal and tangential bonding springs for transferring stress and strain between discrete blocks. The bonding spring parameters used in this embodiment are shown in Table 2.

[0114] Table 2 Parameters of bonded spring

[0115]

[0116] S3: Based on the fracture mechanics theory, a nonlinear failure criterion is established and introduced into the M-DDA model to simulate the deformation and fracture process of the bonded spring. The failure criterion includes two stages: pre-peak and post-peak. In the pre-peak stage, the bond does not fail, and mainly exhibits elastic deformation, with a linear relationship between stress and strain. In the post-peak stage, the bond begins to be damaged, the strength gradually decreases, and the stress-strain relationship is exponential. The application of the nonlinear failure criterion overcomes the shortcomings of existing discontinuous deformation analysis methods that cannot describe the nonlinear softening behavior of materials. The nonlinear failure criterion in the present invention is as follows: Figure 3 shown.

[0117] S41: Enter the time step and start calculation.

[0118] S42: Calculation of block deformations using the M-DDA governing equations :

[0119]

[0120] in, and is the displacement of the center of the block in the x and y directions; is the rotational displacement of the block center; and is the strain of the block in the x and y directions; is the shear strain of the block.

[0121] Calculate the displacement of the spring action point and , i and j represents the two points of action of a set of springs, and The calculation method is the same as For example:

[0122]

[0123] Calculate the displacement transformation matrix based on the coordinates of the centroid and the spring action point .

[0124]

[0125] in, are the coordinates of the block centroid, are the coordinates of the spring action point.

[0126] According to the displacement of the spring action point and , calculate the bonded spring normal phase deformation , Tangential deformation :

[0127]

[0128] in, and are the normal and tangential lengths at the beginning of the current time step, and are the normal and tangent unit vectors, and Spring action points i and j displacement.

[0129] S43: According to the deformation of each group of springs and , calculate the average deformation of the three sets of parallel spring units , r = n , represents the normal deformation, r = s , represents the tangential deformation, and in the following text r The meaning of is deduced from this.

[0130] According to the average deformation Calculation of elastic critical displacement :

[0131]

[0132] in, is the input bond stiffness; (r=n,s) is the bond spring strength, and A is the bond action surface area.

[0133] Further calculation of tensile fracture energy G:

[0134]

[0135] in, is the bonding spring force:

[0136]

[0137] in, is the deformation of the cementing unit; and is the elastic critical displacement and damage parameter; is the input bond stiffness

[0138] Calculating the damage parameters of bonded springs :

[0139]

[0140] Calculating the damage factor of a bonded spring :

[0141]

[0142] In the formula, and is the elastic critical displacement and damage parameter, is the deformation of the bonding unit.

[0143] Determine the damage state of the bonded spring based on the damage coefficient: =0, in the elastic deformation stage, >0, in the stress softening stage, >0.999, the bond is considered to be broken and the spring element is removed.

[0144] S44: According to the damage coefficient Calculating Spring Rate , and the bond spring stiffness is updated in M-DDA, calculated as:

[0145]

[0146] S45: Calculate the spring coefficient matrix and force matrix based on the minimum potential energy theorem, as follows:

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153] in, T is the displacement transformation matrix, , , , is the spring coefficient matrix; is the total energy generated by the cemented spring; , Spring action points , The displacement transformation matrix at ; , and , They are the degrees of freedom of the two blocks acted upon by springs; and are the normal and tangential stiffness of the cemented spring, respectively; and are the normal and tangential unit vectors of the bonded spring; Fi, Fj are the spring force matrices; and are the normal and tangential lengths at the beginning of the current time step; the subscripts i and j are the block numbers.

[0154] S46: Substitute the bond spring stiffness matrix and force matrix into the control equation, and the control equation is:

[0155]

[0156] in, is the mass matrix; is the damping matrix; is the elasticity matrix; is the contact spring coefficient matrix; is the bonding spring coefficient matrix; is the particle contact force and load matrix; is the bonding spring force matrix; and is the velocity and acceleration matrix; D is the displacement matrix.

[0157] The block deformation is calculated and updated according to the updated control equation, and then the calculation loop is returned to step S42.

[0158] S5: Analyze the simulation results, call the dedicated module, convert and store the calculation results according to the standard data architecture of the Visualization Toolkit (VTK) file, and use the ParaView software to parse the exported VTK file to visualize the crack initiation position, crack initiation process, stress-strain curve and stress evolution process.

[0159] The final destruction form of this embodiment is as follows Figure 5 As shown, it is basically consistent with the experimental results. The force-displacement curve obtained in this embodiment is shown in Figure 6 As shown, it can be seen that the method proposed by the present invention shows a nonlinear stress softening phenomenon after the force reaches the peak value, and the stress-strain curve is basically consistent with the test; while the model using the brittle failure criterion is difficult to show stress softening phenomenon after the peak value, and the model using the linear softening criterion is quite different from the test results. The principal stress field evolution diagram obtained in this embodiment is shown in Figure 7 As shown in FIG. 1 , the stress is concentrated at the front end of the crack, driving the crack to expand forward. The displacement field evolution diagram obtained in this embodiment is shown in FIG. Figure 8 shown.

[0160] It will be easily understood by those skilled in the art that the above description is only 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 in the protection scope of the present invention.

Claims

1. A method for predicting rock failure process based on discontinuous deformation analysis, characterized in that: The following steps are involved: S100: Obtain the geometric and physical parameters of the rock, conduct regional discretization modeling based on the geometric characteristics and microstructural development characteristics of the rock, and impose boundary conditions on the model; S200: Based on elastic theory, a bonding spring unit is applied between adjacent discrete blocks to transmit stress and strain information between the blocks; S300: Based on the fracture mechanics theory, a nonlinear failure criterion is established and embedded in the DDA algorithm to determine the degree of damage of the bonded spring; S400: Enter the time step to start the calculation, calculate the normal phase and tangential deformation of the bonded spring unit, update the stress and strain information of the block in real time, then calculate the damage coefficient according to the nonlinear failure criterion described in step S300, and judge the damage degree of the bonded spring unit, remove the bonded spring unit whose damage degree reaches the failure condition, and update the undamaged spring stiffness, and update the block deformation according to the spring stiffness; S5 00: According to the calculation result of step S400, the rock crack initiation position is captured, and the crack initiation process, stress-strain curve and stress evolution process are generated; In step S100, the regional discretization modeling includes: The discretization degree is determined according to the microstructure development degree and geometric characteristics. For areas with relatively developed microstructures or complex geometric shapes, small-sized blocks are used for discretization, and large-sized blocks are used for other areas to establish a discretized numerical model. In step S100, the physical parameters of the rock include: stiffness, density, Poisson's ratio, tensile strength, cohesive strength, and friction angle.

2. A method for predicting rock failure process based on discontinuous deformation analysis according to claim 1, characterized in that: In step S200, three groups of parallel bonded spring units are included, and each group of bonded spring units includes two bonded springs, one in the normal direction and the other in the tangential direction.

3. A method for predicting rock failure process based on discontinuous deformation analysis according to claim 2, characterized in that: In step S200, the material parameters of the bonded spring unit include bonded tensile strength, bonded cohesive strength, bonded tensile fracture energy, bonded shear fracture energy, and bonded spring normal stiffness.

4. The method for predicting rock failure process based on discontinuous deformation analysis according to claim 1 is characterized in that: In step S300, a nonlinear failure criterion is used to judge the damage and fracture process of the reaction rock. The nonlinear failure criterion is expressed as an exponential tensile failure constitutive relationship, which includes two stages: pre-peak and post-peak. In the pre-peak stage, the cementation is not damaged, mainly manifested as elastic deformation, and the stress-strain relationship is linear; in the post-peak stage, the cementation begins to be damaged, the strength gradually decreases, and the stress-strain relationship is exponential.

5. The method for predicting rock failure process based on discontinuous deformation analysis according to claim 4 is characterized in that: Step S400 specifically includes the following steps: S410: Input the time step and start calculation; S420: Calculate bonded spring deformation: Calculate block deformation using the DDA control equation, and then calculate the bonded spring normal based on the block deformation. , tangential deformation The bonded spring phase , tangential deformation The computational models include: in, and are the normal and tangential lengths at the beginning of the current time step, and are the normal and tangent unit vectors, and Spring action points i and j displacement; S430: Calculate the damage coefficient of the bonded spring and determine the degree of damage to the bonded spring: Calculate the damage coefficient based on the deformation of the bonded spring combined with the nonlinear failure criterion , using the damage factor To reflect the degree of bond damage, the damage state of the bonded spring is determined according to the damage coefficient: =0, in the elastic deformation stage, >0, in the stress softening stage, >0.999, the bonding is considered to be damaged and the spring unit is removed. The calculation model of the damage coefficient includes: In the formula, is the deformation of the bonding unit, r = n , represents the normal deformation, r = s , represents the tangential deformation, The calculated value of is equal to the average value of the three groups of spring deformations, and is the elastic critical displacement and damage parameter; S440: Calculate and update the bond spring stiffness based on the damage coefficient Calculate and update spring normal and tangential stiffness ; Updated spring normal, tangential stiffness include: in, is the input bond stiffness; S450: Reconstruct the DDA control equation and calculate and update the block displacement.

6. The method for predicting rock failure process based on discontinuous deformation analysis according to claim 1 is characterized in that: In step S450, the reconstructing the DDA control equation includes: S451: Calculate the spring coefficient matrix and force matrix based on the minimum potential energy theorem; S452: Reconstruct the DDA control equation according to the bond spring stiffness matrix and the force matrix, and calculate and update the block displacement. The calculation model of the DDA control equation is: in, is the mass matrix; is the damping matrix; is the elasticity matrix; is the contact spring coefficient matrix; is the bonding spring coefficient matrix; is the particle contact force and load matrix; is the bonding spring force matrix; and is the velocity and acceleration matrix; D is the displacement matrix.

7. The method for predicting rock failure process based on discontinuous deformation analysis according to claim 1 is characterized in that: Step S500 specifically includes the following steps: converting and storing the calculated geometric structure and corresponding physical quantity data according to the standard data architecture of the VTK file, importing and parsing the generated VTK file through the ParaView platform, and visually analyzing the crack initiation position, crack initiation process, stress-strain curve and stress evolution process.

8. A system for predicting rock failure process based on discontinuous deformation analysis, characterized in that: include: The first main control module is used to obtain the geometric and physical parameters of the rock, conduct regional discretization modeling based on the geometric characteristics and microstructural development characteristics of the rock, and impose boundary conditions on the model; The second main control module is used to apply a bonding spring unit between adjacent discrete blocks for transmitting stress and strain information between the blocks based on elasticity theory; The third main control module is used to establish a nonlinear failure criterion based on fracture mechanics theory and embed it into the DDA algorithm to determine the degree of damage to the bonded spring; The fourth main control module is used to input the time step to start the calculation, calculate the normal phase and tangential deformation of the bonded spring unit, update the stress and strain information of the block in real time, and then calculate the damage coefficient according to the nonlinear failure criterion, determine the damage degree of the bonded spring unit, remove the bonded spring unit whose damage degree reaches the failure condition, and update the undamaged spring stiffness, and update the block deformation according to the spring stiffness; The fifth main control module is used to capture the rock crack initiation position according to the calculation results of the fourth main control module, and generate the crack initiation process, stress-strain curve and stress evolution process; The regional discretization modeling includes: The discretization degree is determined according to the microstructure development degree and geometric characteristics. For areas with relatively developed microstructures or complex geometric shapes, small-sized blocks are used for discretization, and large-sized blocks are used for other areas to establish a discretized numerical model. The physical parameters of the rock include: stiffness, density, Poisson's ratio, tensile strength, cohesive strength, and friction angle.

Images