A method, device, medium and product for simulating a three-dimensional dynamic fracturing process of a rock mass

Abstract

The application discloses a rock mass three-dimensional dynamic cracking process simulation method, equipment, medium and product, and relates to the technical field of rock mass strength and surrounding rock stability evaluation. Firstly, a plurality of three-dimensional elastic-plastic finite strain element finite element models are constructed, then finite element calculation simulation is performed on the finite element models, and whether cracks are generated is determined; then a cohesive crack element is generated between adjacent three-dimensional elastic-plastic finite strain elements of the cracks and is added to a cohesive crack element set; in a multi-thread parallel calculation mode, state parameters and separation states of each cohesive crack element in the cohesive crack element set are obtained, the finite element model is updated, the above process is repeated until convergence. The application simultaneously considers the finite strain element elastic-plasticity and the strain rate effect of the cohesive crack element under the continuous-discontinuous numerical method framework, which has higher precision for simulating the strain rate effect and crack propagation of the rock mass under the action of the power.

Description

Technical Field

[0001] This application relates to the field of rock mass strength and surrounding rock stability evaluation technology, and in particular to a method, equipment, medium and product for simulating the three-dimensional dynamic fracturing process of rock mass. Background Technology

[0002] Currently, blasting excavation remains the most common method for achieving efficient rock mass fracturing, and it has been widely applied to various engineering projects for decades, including hydropower station cavern excavation, nuclear power plant excavation, open-pit / underground mining, and tunnel excavation. However, the testing methods available for observing / quantifying the transient process of fracture within the rock mass are limited, and the understanding of the dynamic fracture process in rock blasting is not deep enough, making the optimization of rock blasting technology extremely challenging. With the rapid development of computational methods and hardware performance in recent decades, continuous-discontinuous numerical methods have been used to quantify the effects of rock blasting. Under dynamic loads such as blasting, fracture analysis based solely on bonded crack elements may not effectively release rock mass stress and reflect the stress state of the rock mass, severely overestimating crack length and crack propagation rate. Simultaneously, under dynamic loads such as blasting, the rock mass is subjected to high strain rate loads, and the mechanical properties of joints exhibit significant strain rate effects; that is, their strength, deformation, and failure modes change with the loading rate. Summary of the Invention

[0003] The purpose of this application is to provide a method, equipment, medium, and product for simulating the three-dimensional dynamic fracturing process of rock mass, so as to improve the accuracy of rock mass fracturing process simulation.

[0004] To achieve the above objectives, this application provides the following solution.

[0005] Firstly, this application provides a method for simulating the rock mass fracture process, including: Finite element modeling is performed on the rock mass to obtain a finite element model; the finite element model includes multiple three-dimensional elastic-plastic finite strain elements; Finite element simulation was performed on the finite element model to determine whether cracks were generated between adjacent three-dimensional elastoplastic finite strain elements. Bonded crack elements are generated between adjacent three-dimensional elastoplastic finite strain elements of the crack and added to the bonded crack element set; The system employs a multi-threaded parallel computing approach to obtain the state parameters and separation states of each bonded crack element in the bonded crack element set, updates the finite element model, and returns to the step of "performing finite element calculation simulation on the finite element model and determining whether cracks are generated between each adjacent three-dimensional elastoplastic finite strain element" until convergence.

[0006] Secondly, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above-described method for simulating the rock mass fracturing process.

[0007] Thirdly, this application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for simulating the rock mass fracturing process.

[0008] Fourthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the above-described method for simulating the rock mass fracturing process.

[0009] According to the specific embodiments provided in this application, this application has the following technical effects.

[0010] This application provides a method, equipment, medium, and product for simulating the three-dimensional dynamic fracturing process of rock mass. Under the framework of continuous-discontinuous numerical methods, this application simultaneously considers the strain rate effects of finite strain elements and bonded crack elements, which provides higher accuracy for simulating the strain rate effect and crack propagation of rock mass under dynamic loading. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0012] Figure 1 This is a schematic flowchart of a rock mass fracture process simulation method provided in an embodiment of this application.

[0013] Figure 2 A schematic diagram of the shape of the Mohr-Coulomb criterion in the π plane and principal stress space provided in an embodiment of this application.

[0014] Figure 3 This is a schematic diagram of a three-dimensional elastoplastic finite strain element provided in an embodiment of this application.

[0015] Figure 4 This is a schematic diagram of a three-dimensional rate-dependent bond crack element provided in an embodiment of this application.

[0016] Figure 5 This is a schematic diagram illustrating the coupling process between a three-dimensional elastoplastic finite strain element and a bonded crack element, provided in an embodiment of this application.

[0017] Figure 6The image shows the uniaxial compression simulation results of a three-dimensional rock mass sample provided in an embodiment of this application.

[0018] Figure 7 This is a schematic diagram illustrating the principle of contact judgment and contact force calculation provided in an embodiment of this application.

[0019] Figure 8 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0020] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0021] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0022] In one exemplary embodiment, a method for simulating rock mass fracturing processes is provided, such as... Figure 1 As shown, it includes the following steps 101-104.

[0023] Step 101: Perform finite element modeling on the rock mass to obtain a finite element model; the finite element model includes multiple three-dimensional elastoplastic finite strain elements; Step 102: Perform finite element calculation simulation on the finite element model and determine whether cracks are generated between adjacent three-dimensional elastoplastic finite strain elements. Step 103: Generate bonded crack elements between adjacent three-dimensional elastoplastic finite strain elements of the crack and add them to the bonded crack element set; Step 104: Using multi-threaded parallel computing, obtain the state parameters and separation states of each bonded crack element in the bonded crack element set, update the finite element model, and return to the step of "performing finite element calculation simulation on the finite element model and determining whether cracks are generated between each adjacent three-dimensional elastoplastic finite strain element" until convergence.

[0024] Implementing steps 101-104 above can achieve the following technical effects.

[0025] Continuous-discontinuous numerical methods must consider both elasto-plastic and strain rate effects, which are crucial for assessing the stability of surrounding rock in deep caverns under static-dynamic coupling. This application, within the framework of continuous-discontinuous numerical methods, simultaneously considers the elasto-plastic effect of finite strain elements and the strain rate effect of bonded crack elements. This provides higher accuracy for simulating the strain rate effect and crack propagation of rock mass under dynamic loading, and is significant for accurately characterizing the dynamic strength properties and crack propagation mechanism of rock mass, revealing the influence of dynamic disturbance on the rock mass fracturing process. It also has important engineering application value for the stability evaluation and support design of surrounding rock in deep caverns.

[0026] In another exemplary embodiment, the specific implementation process of step 101 above is as follows: By embedding the tensioned truncated Mohr-Coulomb model into the finite strain element (tetrahedral element) within the framework of continuous-discontinuous numerical methods, a finite element model comprising multiple three-dimensional elastoplastic finite strain elements is obtained. However, in the principal stress space, the intersection of the tensile truncated surface and the Mohr-Coulomb yield surface forms an irregular hexagon, such as... Figure 2 As shown, Figure 2 (a) is the main view. Figure 2 (b) in the figure is a 3D image. To improve numerical convergence, a piecewise linear-spline hybrid smoothing algorithm is used to construct a local smoothing function in the neighborhood of the corner point.

[0027] To improve computational efficiency, a GPU parallel processing method is constructed. First, the corner coordinates are pre-stored. The coordinates of the six corner points of the irregular hexagon are encapsulated into a 6×3 coordinate array, along with the internal friction angle, cohesion, tensile strength, neighborhood radius, and attenuation coefficient. All of these are then input into the GPU and cached globally in read-only mode. Next, Gaussian integration points (i.e., ...) of all three-dimensional elastoplastic finite strain elements are extracted in batches from the numerical simulation software. Figure 3 The three principal stress values ​​of N material points (a total of N points) are encapsulated into an N×3 array, along with key parameters (corner neighborhood radius, weight decay parameter), and input into the GPU at once. Then, the GPU automatically distributes the N integration points in the stress array to N threads for parallel computation. The computation includes calculating the Euclidean distance from the Gaussian integration point to the six corner points. The calculation is shown in Formula (1); the weight calculation is shown in Formula (2). Within the corner neighborhood, an exponentially decaying weight is used; outside the corner neighborhood, no smoothing is performed, and the original stress is used directly. The stress value of the current integration point is calculated by weighted average based on the six principal stresses at the corner points, as shown in Formula (3). It is determined whether the current integration point is within the corner neighborhood and the corresponding smoothing function is selected, as shown in Formula (4), where piecewise linear interpolation is used outside the corner neighborhood, and a cubic B-spline local smoothing function is used within the corner neighborhood. The stress value is calculated after all threads have completed their calculations. The stress array is batch-cached into the GPU's global memory, and then synchronized to the integration points of the numerical simulation software via the CPU for subsequent stress updates and iterative calculations. The constructed three-dimensional elastoplastic finite strain element is as follows: Figure 3 As shown.

[0028] (1); (2); (3); (4); in, The principal stress vector at the Gaussian integration point. Let K be the principal stress vector at the k-th corner point. , , These are the three principal stress values ​​at the Gaussian integration point. , , These are the three principal stress values ​​at the k-th corner point. The weight of the k-th corner point. for, for, The smoothed reference stress is obtained by weighted interpolation using the principal stresses at the six corner points of the hexagon as the reference field. To smooth the corner stress after transformation, It is a cubic B-spline smoothing function.

[0029] In another exemplary embodiment, in step 102 above, material points are set within the finite strain element to track the material's historical information (such as plastic strain). The stress and strain at the material points of the first-order tetrahedral element are calculated using the finite element method. Inverse extrapolation is performed using shape functions. Because the stress is uniform across the entire element field, the stress at the integration point equals the stress across the entire element field. All four nodes are assigned the same stress and strain values. Potential damage nodes are determined (first threshold). (Generally taken as 0.7) and the actual crack initiation judgment (second threshold) ) controls the crack initiation location, namely two threshold parameters ( and The stress and strain state equations are controlled. First, the state of each node is checked to determine whether it meets the first threshold condition (e.g., whether the stress or strain exceeds a certain limit). Nodes that meet this condition are marked as "potentially damaged nodes", as shown in Equations (5) and (6). For each potentially damaged node, the stress (normal stress and shear stress) at each connection edge / face of the node is calculated to determine whether it meets the second threshold condition, so as to determine the location and direction of the actual crack, as shown in Equations (7) and (8). If it is determined to be an initiating crack, the original node is split into multiple new nodes, and at the same time, bonded crack elements are generated at the damage site (edge ​​of finite strain element). The generated three-dimensional rate-dependent bonded crack elements are as follows: Figure 4 As shown.

[0030] (5) (6) in, This is the function for determining potential damaged nodes. Let be the strain determination function. This is the tensile strength criterion function. This is the shear strength determination function. The equivalent plastic strain of the node is denoted as . The characteristic plastic strain of the node, This represents the maximum principal stress at the node (negative under compression, positive under tension). The first threshold, For the tensile strength of the material, The maximum shear stress at the node. For the shear strength of the material, The normal stress is perpendicular to the direction of maximum shear stress. is the coefficient of friction.

[0031] (7); (8); in, This is the function for determining the potential damage surface. This is the tensile strength criterion function. This is the shear strength determination function. The normal stress of the target connection surface (negative under compression, positive under tension). The second threshold, For the tensile strength of the material, The shear stress at the target connection surface. For the shear strength of the material, is the coefficient of friction.

[0032] In another exemplary embodiment, step 102 described above may be implemented using steps 201-205.

[0033] Step 201: Perform finite element calculation simulation on the finite element model to obtain the stress and strain of the material points of each three-dimensional elastoplastic finite strain element; Step 202: Determine the stress and strain of each node of each three-dimensional elastoplastic finite strain element based on the stress and strain of the material points of each element. Step 203: Determine whether the stress and strain of each node of each three-dimensional elastoplastic finite strain element satisfy the first threshold condition (as shown in formulas (5) and (6) above). Step 204: Determine the nodes that meet the first threshold condition as potential damaged nodes; Step 205: Determine whether the normal stress and shear stress of the target connection surface where each potential damage node is located meet the second threshold condition (as shown in formulas (7) and (8) above); the target connection surface is the surface that connects two adjacent three-dimensional elastoplastic finite strain elements; Step 206: Determine that cracks are generated between adjacent three-dimensional elastoplastic finite strain elements corresponding to the target connection surface that meets the second threshold condition.

[0034] In another exemplary embodiment, step 207 above is the process of calculating strain rate and updating dynamic mechanical parameters. To improve calculation speed, strain rate is calculated only after the bonded crack element is generated, and the pre-calculated strain rate array is stored in global memory. A kernel function for strain rate is designed, which calculates the rate state of all joint surfaces. One thread processes one bonded crack element, and each CUDA thread block can process multiple bonded crack elements. Adjacent threads within a thread bundle access data of consecutive integration points within the same element. If some bonded crack elements are completely destroyed, the mechanical update still needs to be executed but is skipped, ensuring a unified execution path for the thread bundle. The state of the bonded crack element is updated based on the mechanical state of the finite strain element, dynamic mechanical parameters, and fracture criteria. The coupling process between the three-dimensional elastoplastic finite strain element and the bonded crack element is as follows: Figure 5 As shown, the uniaxial compression simulation results of the three-dimensional rock mass sample are as follows: Figure 6 As shown.

[0035] In another exemplary embodiment, the above-mentioned state parameters include joint separation rate, strain rate and strength parameters, and the separation state is before cracking, during cracking and after cracking; the above-mentioned multi-threaded parallel computing method is to set up multiple threads, each thread calculates one bonded crack element in the set of bonded crack elements, wherein the execution process of one thread includes the following steps 301-303.

[0036] Step 301: Calculate the joint separation rate, strain rate, and strength parameters of the bonded crack element based on the strain rate of each node on the joint surface connecting the bonded crack element and the two three-dimensional elastoplastic finite strain elements. Step 302: Based on the joint separation rate, strain rate, and strength parameters of the bonded crack element, the separation state of the bonded crack element is determined using fracture criteria; the separation state is before cracking, during cracking, and after cracking. Step 303: Update the state of the two three-dimensional elastoplastic finite strain elements connected to the bonded crack element based on the joint separation rate, strain rate, strength parameters and separation state of the bonded crack element.

[0037] In another exemplary embodiment, in step 301 above, the strain rate is calculated by dividing the separation rate of the tetrahedral finite strain elements on both sides of the bonded crack element (i.e., the three-dimensional elastoplastic finite strain element mentioned above; in this embodiment, the three-dimensional elastoplastic finite strain element is a tetrahedral finite strain element) by the characteristic size of the bonded crack element. A bonded crack element includes six nodes: N0, N1, and N2 on the upper surface, and N3, N4, and N5 on the lower surface. The velocity of each joint element face (triangle) is calculated using the velocities of the three nodes and the shape function. Interpolation yields that the surface nodal velocities of the joint element inherit the velocities of the corresponding nodes of the original tetrahedron, and the shape function satisfies... and The joint separation rate is calculated as shown in equation (9), and the strain rate is calculated as shown in equation (10).

[0038] (9); (10); in, The joint separation rate of the bonded crack element. , , The strain velocities are represented by three nodes on the joint surface where a bonded crack element is connected to a three-dimensional elastoplastic finite strain element. , , The strain velocities are represented by the strain velocities at three nodes on the joint surface connecting the bonded crack element to another three-dimensional elastoplastic finite strain element. , , To generate the shape function that yields the velocity of the joint surface; For the strain rate of the bonded crack element, The characteristic dimension of the bonded crack element is the average of the sum of the side lengths of the tetrahedral element and the triangular element connected to the bonded crack element.

[0039] Then, the strength parameters of the bonded crack element, including cohesion, are updated based on the strain rate. ),tensile strength( Type I fracture energy, Type II fracture energy. Dynamic cohesion ( The relationship between the cohesion and the strain rate is shown in equation (11), and the dynamic enhancement coefficient of cohesion at each time step is: As shown in equation (12), the dynamic enhancement factor of tensile strength is also set to Dynamic tensile strength The calculation is shown in equation (13), while ensuring that the updated tensile strength does not exceed the theoretical upper limit, i.e. The dynamic enhancement coefficients for type I and type II fracture energies are set according to empirical formulas. .

[0040] (11) (12) (13) In the formula, For dynamic enhancement coefficient, An empirical constant related to the size of the bonded crack element. For the strain rate of the bonded crack element, The static strain rate, For quasi-static cohesion, For dynamic tensile strength, This represents the quasi-static tensile strength of the bonded crack element.

[0041] In another exemplary embodiment, the newly generated bonded crack elements are given the constitutive relation of the cohesive crack model, and the behavior of the crack surface is controlled by the traction-separation law: ① Before cracking: It can transmit stress and maintain the interaction between the materials on both sides of the crack.

[0042] ② During the cracking process: Simulate the softening and damage evolution of materials.

[0043] The fracture criterion is as follows: when the joint opening displacement of the bonded crack element is less than the critical displacement, the separation state is determined to be before cracking; when the joint opening displacement of the bonded crack element is greater than or equal to the critical displacement and less than the limit displacement, the separation state is determined to be during cracking. Figure 5 As shown, the specific implementation method of this fracture criterion is as follows: If the joint surface of the bonded crack element opens and shifts ( ) exceeding the tensile strength ( The critical displacement defined by () ),Right now Tensile yielding failure occurs, at which point cracking begins. When the joint surface opening displacement of the bonded crack element exceeds the ultimate displacement ( ),Right now The joint element is completely destroyed, the cracking is completed, and the two adjacent three-dimensional elastoplastic finite strain elements are completely separated into two discrete elements.

[0044] If the sliding displacement of the joint element ( ) exceeding the shear strength ( The critical displacement defined by () ),Right now Then shear failure occurs. When the sliding displacement of the joint element exceeds the limit displacement ( When the bond is separated tangentially, the tangential adhesive stress becomes the frictional force defined by the Coulomb model.

[0045] set up The material softening coefficient: (14) In the formula, a , b and q These are the empirical curve fitting parameters, which are 0.63, 1.8, and 6.0, respectively. D It is the damage coefficient, and its calculation method is as follows: (15) (16) (17) (18) (19) in, The normal stress at the target connection surface. This represents the shear stress on the joint surface of the bonded crack element.

[0046] After complete cracking: the transmitted stress drops to zero, and the two sides of the crack become independent free surfaces, allowing contact mechanical behaviors such as collision, slippage, and separation to occur.

[0047] Finally, the entire model domain will consist of a continuous medium domain, a discrete medium domain, and a medium domain where damage is about to occur. The continuous medium domain still follows step 102; the discrete medium domain requires contact determination and contact force calculation. The contact determination uses the NBS method, such as... Figure 7As shown, the normal contact force is calculated using a potential-based method, as shown in formula (14). The frictional force between the contact pairs is calculated using Coulomb's law of friction, as shown in formulas (15) and (14). The damage to the medium domain is calculated using the two threshold parameter methods in steps 201-206.

[0048] (20) In the formula, For normal contact force, For normal contact stiffness, Let be the contact potential gradient inside the Contactor block. The contact potential gradient inside the Target block is given. Let be any infinitesimal element in the contact region.

[0049] (twenty one) In the formula, For the frictional force between the contact pairs, For tangential contact stiffness, For friction, For the Rodriguez rotation tensor, The relative displacement increment at the current time step.

[0050] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 8 As shown, this computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides the environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The I / O interfaces are used for exchanging information between the processor and external devices. The communication interface is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements a method for simulating a rock mass fracturing process.

[0051] Those skilled in the art will understand that Figure 8 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0052] In one exemplary embodiment, a computer device is provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0053] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0054] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0055] According to the specific embodiments provided in this application, this application has the following technical effects.

[0056] This application presents an efficient continuous-discontinuous numerical method that considers the strain rate effects of finite strain elements and bonded crack elements. Within the framework of continuous-discontinuous numerical methods, it is the first to simultaneously consider the strain rate effects of finite strain elements and bonded crack elements, which provides higher accuracy for simulating the strain rate effect of rock mass and crack propagation under dynamic loading.

[0057] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties.

[0058] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0059] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, distributed databases based on blockchain. The processors involved in the embodiments provided in this application may be, but are not limited to, general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc.

[0060] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0061] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A method for simulating the three-dimensional dynamic fracturing process of rock mass, characterized in that, include: Finite element modeling is performed on the rock mass to obtain a finite element model; The finite element model includes multiple three-dimensional elastic-plastic finite strain elements; Finite element simulation was performed on the finite element model to determine whether cracks were generated between adjacent three-dimensional elastoplastic finite strain elements. Bonded crack elements are generated between adjacent three-dimensional elastoplastic finite strain elements of the crack and added to the bonded crack element set; The system employs a multi-threaded parallel computing approach to obtain the state parameters and separation states of each bonded crack element in the bonded crack element set, updates the finite element model, and returns to the step of "performing finite element calculation simulation on the finite element model and determining whether cracks are generated between each adjacent three-dimensional elastoplastic finite strain element" until convergence.

2. The method for simulating rock mass fracturing process according to claim 1, characterized in that, Finite element analysis was performed on the finite element model to determine whether cracks were generated between adjacent three-dimensional elastoplastic finite strain elements. Specifically, this included: Finite element simulation is performed on the finite element model to obtain the stress and strain of the material points of each three-dimensional elastoplastic finite strain element; Based on the stress and strain of the material points of each three-dimensional elastoplastic finite strain element, determine the stress and strain of each node of each three-dimensional elastoplastic finite strain element. Determine whether the stress and strain at each node of each three-dimensional elastoplastic finite strain element satisfy the first threshold condition; Nodes that meet the first threshold condition are identified as potentially damaged nodes; Determine whether the normal stress and shear stress of the target connection surface where each potential damage node is located meet the second threshold condition; the target connection surface is the surface connecting two adjacent three-dimensional elastoplastic finite strain elements; It is determined that cracks are generated between adjacent three-dimensional elastoplastic finite strain elements corresponding to the target connection surface that meets the second threshold condition.

3. The method for simulating rock mass fracturing process according to claim 2, characterized in that, The first threshold condition is: ； ； in, This is the function for determining potential damaged nodes. Let be the strain determination function. This is the tensile strength criterion function. This is the shear strength determination function. The equivalent plastic strain of the node is denoted as . The characteristic plastic strain of the node, The maximum principal stress at the node is... The first threshold, For the tensile strength of the material, The maximum shear stress at the node. For the shear strength of the material, The normal stress is perpendicular to the direction of maximum shear stress. is the coefficient of friction.

4. The method for simulating rock mass fracturing process according to claim 2, characterized in that, The second threshold condition is: ； ； in, This is the function for determining the potential damage surface. This is the tensile strength criterion function. This is the shear strength determination function. The normal stress at the target connection surface, The second threshold, For the tensile strength of the material, The shear stress at the target connection surface. For the shear strength of the material, is the coefficient of friction.

5. The method for simulating rock mass fracturing process according to claim 1, characterized in that, The state parameters include joint separation rate, strain rate and strength parameters, and the separation state is before cracking, during cracking and after cracking; The multi-threaded parallel computing method involves setting up multiple threads, each of which calculates one bonded crack element from the set of bonded crack elements. The execution process of one thread is as follows: Based on the strain rate of each node on the joint surface where the bonded crack element is connected to two three-dimensional elastoplastic finite strain elements, the joint separation rate, strain rate and strength parameters of the bonded crack element are calculated. Based on the joint separation rate, strain rate, and strength parameters of the bonded crack element, a fracture criterion is used to determine the separation state of the bonded crack element; the separation state is defined as before cracking and during cracking; the fracture criterion is as follows: when the joint opening displacement of the bonded crack element is less than the critical displacement, the separation state is determined to be before cracking; when the joint opening displacement of the bonded crack element is greater than or equal to the critical displacement and less than the ultimate displacement, the separation state is determined to be during cracking. The states of the two three-dimensional elastoplastic finite strain elements connected to the bonded crack element are updated based on the joint separation rate, strain rate, strength parameters, and separation state of the bonded crack element.

6. The method for simulating rock mass fracturing process according to claim 5, characterized in that, The formula for calculating the joint separation rate is: ； in, The joint separation rate of the bonded crack element. , , The strain velocities are represented by three nodes on the joint surface where a bonded crack element is connected to a three-dimensional elastoplastic finite strain element. , , The velocities of three nodes on the joint surface connecting the bonded crack element to another three-dimensional elastoplastic finite strain element are given. , , Let be the shape function of the three-node joint surface; The formula for calculating strain rate is: ； in, For the strain rate of the bonded crack element, This refers to the dimensional characteristics of the bonded crack element.

7. The method for simulating rock mass fracturing process according to claim 5, characterized in that, The strength parameters are the cohesion dynamic enhancement coefficient, tensile strength dynamic enhancement coefficient, Type I fracture energy dynamic enhancement coefficient, Type II fracture energy dynamic enhancement coefficient, and dynamic tensile strength. The cohesion dynamic enhancement coefficient and tensile strength dynamic enhancement coefficient are the same, both being dynamic enhancement coefficients. The Type I fracture energy dynamic enhancement coefficient and Type II fracture energy dynamic enhancement coefficient are the same, both being the square of the first dynamic enhancement coefficient. The dynamic strengthening factor of cohesion and the dynamic strengthening factor of tensile strength are calculated using the following formulas: ； in, For dynamic enhancement coefficient, An empirical constant related to the size of the bonded crack element. For the strain rate of the bonded crack element, This is the static strain rate.

8. A computer device, comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the rock mass fracturing process simulation method according to any one of claims 1-7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the rock mass fracture process simulation method as described in any one of claims 1-7.

10. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the rock mass fracture process simulation method as described in any one of claims 1-7.

Images