Three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuities

Abstract

The application provides a three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuous surface. The method comprises the following steps: generating surface domains used for representing various discontinuous surfaces in the rock mass as grid surfaces in a three-dimensional rock model region, and dividing the grid surfaces into tetrahedral solid element grids to obtain tetrahedral solid elements of the three-dimensional rock model region; updating node topological structures of the tetrahedral solid elements to obtain a plurality of joint elements in the three-dimensional rock model region; at an initial time step, determining stress information of each joint element according to a type identifier of each joint element, and adopting a stress algorithm matched with the type identifier to determine the stress information of each joint element; determining total equivalent node forces on each node in the plurality of joint elements based on the stress information of each joint element; updating the type identifier of each joint element in the plurality of joint elements based on the total equivalent node forces, and updating the total equivalent node forces of each node in the next time step to update the type identifier of each joint element.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geotechnical engineering, energy and environment, and particularly relates to a three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuities. BACKGROUND

[0002] A system of discontinuities such as fractures, bedding planes, and joints widely exists in natural rock mass. The seepage softening effect of the system directly affects the engineering stability problems caused by the attenuation of rock mass self-supporting force (such as underground cavern collapse and slope instability), and is a key factor restricting the efficiency of shale gas reservoir reconstruction and deep mining. The softening process is essentially a fluid-structure interaction process of the mutual interaction between the deterioration of mechanical properties of fracture surfaces and the deformation of the medium under the driving of fluid pressure.

[0003] However, the existing fluid-structure interaction theory model and technology of rock mass are not universal, and their completeness needs to be optimized and supplemented. SUMMARY

[0004] The present application provides a three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuities.

[0005] In a first aspect, the present application provides a three-dimensional numerical simulation method for seepage softening of rock mass discontinuities, comprising:

[0006] Generating a surface domain used to represent various discontinuities in the rock mass as a grid surface in the rock three-dimensional model region, and dividing the grid surface into tetrahedral solid element grids to obtain tetrahedral solid elements of the rock three-dimensional model region;

[0007] Updating the node topology of the tetrahedral solid elements to obtain a plurality of joint elements in the rock three-dimensional model region;

[0008] At an initial time step, according to the type identification of each joint element in the plurality of joint elements, a stress algorithm matched with the type identification is used to determine the stress information of each joint element; wherein the type of the joint element is a hydraulic joint element, and the stress information of the hydraulic joint element includes the contact force of the two walls of the hydraulic joint element calculated based on the softening process and the first fluid pressure of each node of the hydraulic joint element calculated based on the seepage;

[0009] Based on the stress information of each joint element, the total equivalent node force on each node in the plurality of joint elements is determined;

[0010] Based on the total equivalent node force on each node in the plurality of joint elements, the type identification of each joint element in the plurality of joint elements is updated, and the total equivalent node force of each node in the next time step is updated to update the type identification of each joint element.

[0011] In a second aspect, the embodiments of the present application provide a three-dimensional numerical simulation system for seepage softening of discontinuities in rock mass, comprising:

[0012] a preprocessing module configured to generate, in a three-dimensional rock model region, a surface domain used to represent various discontinuities in the rock mass as a grid surface, and subdivide the grid surface into a tetrahedral solid element grid to obtain a tetrahedral solid element of the three-dimensional rock model region;

[0013] a node topology updating module configured to update a node topology of the tetrahedral solid element to obtain a plurality of joint elements in the three-dimensional rock model region;

[0014] a calculation module configured to, at an initial time step, determine stress information of each of the joint elements according to a type identifier of each of the joint elements in the plurality of joint elements by using a stress algorithm matched with the type identifier; wherein the type of the joint element is a hydraulic joint element, and the stress information of the hydraulic joint element includes a contact force of two walls of the hydraulic joint element calculated based on a softening process and a first fluid pressure of each node of the hydraulic joint element calculated based on seepage;

[0015] the calculation module is further configured to determine a total equivalent node force on each node in the plurality of joint elements based on the stress information of each of the joint elements;

[0016] a type updating module configured to update the type identifier of each of the joint elements in the plurality of joint elements based on the total equivalent node force on each node in the plurality of joint elements, and update the type identifier of each of the joint elements by using the total equivalent node force of each node in a next time step.

[0017] In a third aspect, the embodiments of the present application provide an electronic device, comprising:

[0018] at least one processor; and

[0019] a memory connected with the at least one processor in communication; wherein

[0020] the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the method of the first aspect.

[0021] In a fourth aspect, the embodiments of the present application provide a storage medium storing instructions, when the instructions are run on an electronic device, the electronic device performs the method of the first aspect.

[0022] In a fifth aspect, an embodiment of the present application provides a program product, the program product comprising at least one of a program, an instruction, which, when executed by a processor, implements the steps of the method of the preceding first aspect.

[0023] According to the technical scheme of the present application, the time-dependent seepage-softening model considering the geometric and intensity characteristics of the discontinuous surface can intuitively reflect the whole process of the shear mechanical properties of the discontinuous surface of the rock mass in tangential sliding, and the application range of the fluid-structure coupling model of the rock mass is widened.

[0024] Additional aspects and advantages of the present application will be made apparent by the following description and the accompanying drawings. BRIEF DESCRIPTION OF DRAWINGS

[0025] The above and / or additional aspects and advantages of the present application will become apparent and be readily understood by considering the following detailed description, including the accompanying drawings, in which:

[0026] Figure 1 A flowchart of a three-dimensional numerical simulation method of seepage softening of a discontinuous surface of a rock mass provided by an embodiment of the present application is shown in the figure.

[0027] Figure 2 An effect diagram of tetrahedral mesh partitioning provided by an embodiment of the present application is shown in the figure.

[0028] Figure 3 An example diagram of a three-dimensional joint element topology provided by an embodiment of the present application is shown in the figure.

[0029] Figure 4 A definition diagram of tangential friction force provided by an embodiment of the present application is shown in the figure.

[0030] Figure 5 An example diagram of a triangular surface seepage channel in a global and local coordinate system provided by an embodiment of the present application is shown in the figure.

[0031] Figure 6 An equivalent node force diagram for fluid pressure calculation provided by an embodiment of the present application is shown in the figure.

[0032] Figure 7 A block diagram of a three-dimensional numerical simulation system of seepage softening of a discontinuous surface of a rock mass provided by an embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0033] Embodiments of the present application are described in detail below with reference to the accompanying drawings, examples of which are shown in the figures, wherein the same or similar notations represent the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to explain the present application, and cannot be understood as limiting the present application.

[0034] It should be noted that the discontinuous surface system such as fissure, bedding, joint, etc. is widely developed in natural rock mass. The seepage softening effect not only directly affects the engineering stability problem caused by the attenuation of rock mass self-supporting force (such as underground cavern collapse, slope instability), but also is the key factor restricting the efficiency of shale gas reservoir reconstruction, deep mining and other energy development. The essence of this softening process is the fluid-structure interaction process of the mutual interaction between the mechanical properties degradation of fissure surface and the deformation of medium under the driving of fluid pressure. With the progress of high-performance computing technology, discrete element method (DEM) and discontinuous deformation analysis (DDA) have become important tools for studying this problem. Existing methods usually simplify the discontinuous surface as a parallel plate model, establish the relationship between opening and flow based on the cubic law, and realize the fluid-structure coupling calculation by iteratively updating the fluid pressure and fissure opening. However, the current research still has significant gaps in the fine characterization of seepage softening mechanism, especially in the analysis of softening effect in three-dimensional complex environment. The existing numerical model mainly aims at the joint rock mass with weak water sensitivity, and its core limitations are: (1) only considering the fluid mechanics effect of fissure surface, ignoring the physical and chemical softening mechanisms such as mineral expansion and cementation weakening of high water sensitivity rock mass during seepage process; (2) using two-dimensional or quasi-three-dimensional simplified model, which cannot truly reflect the coupling effect of three-dimensional seepage field and stress field in complex fissure network; (3) lacking the dynamic tracking ability of rock mass deformation localization and fracture expansion in the softening process. When facing typical high water sensitivity rock mass such as water-rich fault zone and mudstone interbedding, the existing method cannot accurately describe the strength attenuation law of fissure surface caused by seepage, nor can it simulate the multi-scale fracture evolution process induced by softening, leading to deviation in engineering stability assessment.

[0035] In view of the above problems, it is urgent to construct a fine model that can consider the coupling of rock mass material softening characteristics and three-dimensional seepage field. The core advantage of three-dimensional seepage softening model is that it can accurately depict the spatial difference of rock mass softening effect under different seepage directions and stress states by introducing softening parameters to characterize the strength attenuation law of fissure surface under seepage action and describing three-dimensional flow field distribution based on the seepage theory of pore-fracture double medium. This model not only breaks through the limitations of traditional parallel plate assumption, but also reveals the complex mechanical behavior of fissure opening-slippage-penetration induced by seepage softening under three-dimensional stress state, providing more reliable theoretical support for deep underground engineering support design, hydraulic fracturing fracture propagation prediction, etc. Therefore, the theoretical construction and numerical implementation of three-dimensional seepage softening model is not only a key link to fill the gap in current fluid-structure coupling research, but also an important technical path to promote the transformation of rock mass engineering stability analysis from qualitative evaluation to quantitative prediction.

[0036] Fluids change the physical and mechanical properties of different rock types to different degrees. Existing fluid-structure coupling models of rock discontinuities do not focus on the modification process of fluids on different types of rocks, but uniformly apply the effect of fluids on rocks to the rock skeleton boundary as fluid pressure. Although the above method has certain applicability to some low water sensitivity or non-water sensitivity rocks (such as granite, marble, etc.), for the widely developed cemented rocks (such as mudstone, shale, argillaceous sandstone, etc.) in some mountainous areas and coal measures, the existing technology often makes the calculation results better than the actual field results, which has a direct negative impact on the prediction accuracy and control strategy selection of engineering soft rock disasters. Therefore, the existing rock fluid-structure coupling theory model and technology are not universal, and its completeness needs to be optimized and supplemented. In view of the shortcomings of the existing technology, the present application proposes a three-dimensional numerical simulation method and system for the seepage softening of rock discontinuities. This method can accurately match the corresponding rock type through a unified form of softening equation, thereby improving the robustness and applicability of the fluid-structure coupling model in the field of rock mechanics and rock engineering.

[0037] The modification of fluids on the physical and mechanical properties of rocks has significant lithology difference, and the existing fluid-structure coupling model of rock discontinuities generally ignores this core feature, which is simplified as a calculation method of uniformly acting on the rock skeleton as a fluid pressure boundary condition. Although it has certain applicability to low water sensitivity rocks such as granite and marble, it has significant defects in water-sensitive cemented rocks such as mudstone, shale, and argillaceous sandstone widely distributed in some mountainous areas. These rocks will undergo physical and chemical softening reactions such as clay mineral swelling and cement dissolution during seepage, resulting in nonlinear decay of mechanical strength, while the "homogenization" processing method of the existing model often makes the calculation results optimistic, directly affecting the prediction accuracy and control strategy effectiveness of soft rock engineering disasters (such as large deformation of tunnels and water inrush from coal seam floor).

[0038] There are three key gaps in current research: first, a multi-mechanism coupling model of fluid-rock interaction has not been established, and the softening paths (swelling type, dissolution type, ion exchange type) dominated by different mineral components (such as montmorillonite and illite) lack differentiated representation; second, two-dimensional plane or quasi-three-dimensional simplification assumptions are used, which cannot capture the anisotropic influence of three-dimensional seepage field (such as seepage direction, flow velocity gradient, pressure distribution) on the softening process in complex fracture networks; third, a "seepage-softening-deformation" fully coupled framework has not been established, making it difficult to describe the dynamic feedback process of softening-induced crack propagation, stress redistribution, and seepage channel evolution. This theoretical gap is particularly prominent in coal measure strata gas extraction, shale gas horizontal well fracturing, and other engineering, often leading to the failure of early warning of softening-induced wellbore instability and crack penetration.

[0039] To address the aforementioned challenges, constructing a three-dimensional seepage softening model has irreplaceable engineering value. This model, by introducing a lithology-sensitive softening constitutive equation, can accurately characterize the strength attenuation laws of different rock types under three-dimensional seepage conditions. Through coupled discrete fracture network (DFN) modeling technology, it achieves three-dimensional dynamic simulation of fracture initiation, propagation, and connection during the softening process. This model not only provides a precise tool for the stability analysis of surrounding rock in hydropower projects under certain complex geological conditions, but also supports the construction of a "precise prediction-targeted prevention and control" technical system for softening-type disasters in coal-bearing strata, fundamentally solving the problem of insufficient universality of existing models and promoting the leap from "simplified simulation" to "realistic representation" in rock mass fluid-structure interaction theory. Specifically, the three-dimensional numerical simulation method and system for seepage softening of rock discontinuities according to embodiments of this application are described below with reference to the accompanying drawings.

[0040] It should be noted that the execution entity of the three-dimensional numerical simulation method for seepage softening of rock discontinuities in this application embodiment can be a three-dimensional numerical simulation system for seepage softening of rock discontinuities. This system can be implemented by software and / or hardware, and can be configured in an electronic device. For example, the electronic device may include, but is not limited to, a terminal, a server, etc.

[0041] Figure 1 This is a flowchart illustrating the three-dimensional numerical simulation method for seepage softening of rock discontinuities provided in this application embodiment. Figure 1 As shown, the three-dimensional numerical simulation method for seepage softening of the discontinuous rock mass can include, but is not limited to, the following steps.

[0042] In step 101, a surface region is generated within the three-dimensional rock model region to represent various discontinuous surfaces in the rock mass as a mesh surface, and the mesh surface is divided into tetrahedral solid element meshes to obtain tetrahedral solid elements of the three-dimensional rock model region.

[0043] For example, such as Figure 2 As shown, taking a certain type of 3D rock model as an example, surface regions representing various discontinuities such as fractures, bedding, and joints in the rock mass can be generated within the 3D rock model area. These surface regions of various discontinuities are then used as mesh surfaces and divided into tetrahedral solid element meshes, thus obtaining the tetrahedral solid elements of the 3D rock model area. The rock type can be any rock type. For example, the rock type can be any of the following: mudstone, shale, argillaceous sandstone, granite, or marble.

[0044] In step 102, the node topology of the tetrahedral solid element is updated to obtain multiple joint elements in the rock 3D model region.

[0045] For example, the node topology of the tetrahedral solid elements in the region of the rock 3D model can be updated so that all the tetrahedral solid elements have independent node numbers, and then the six-node joint elements (N1, N2, N3, N4, N5, N6) shown in FIG. 1 can be embedded therebetween to simulate the continuous-discontinuous process. Figure 3 The six-node joint element (N1, N2, N3, N4, N5, N6) shown in FIG. 1 can be used to simulate the continuous-discontinuous process, and the joint element is divided into a broken joint element representing a discontinuous surface and a cohesive joint element representing a continuous surface. After calculation, the cohesive joint element can be converted into a broken joint element. Figure 3 The boundary region between the three triangular pyramids can be regarded as a joint element. For example, the boundary region between the surface composed of a1, b1, and c1 and the surface composed of a2, b2, and c2 can be regarded as a joint element. If the two surfaces represent continuous surfaces, the type of the joint element is a cohesive joint element. If the cohesive joint element is softened by fluid seepage to form cracks, expansion, and penetration, resulting in the two surfaces representing discontinuous surfaces, the type of the joint element is updated from a cohesive joint element to a broken joint element. For example, the cohesive joint element on the hydraulic boundary can be marked as a hydraulic joint element (HMJ), and only the broken joint elements connected to the HMJ are subjected to the seepage softening calculation of the discontinuous surface.

[0046] In step 103, at the initial time step, the stress information of each joint element is determined by using a stress algorithm matched with the type identification of each joint element in the plurality of joint elements.

[0047] In the embodiments of the present application, at the initial time step, the discontinuous surfaces such as cracks, bedding, and joints in the region of the rock 3D model can be marked as broken joint elements, and the continuous surfaces between the remaining tetrahedral solid elements are marked as cohesive joint elements. The cohesive joint elements on the water boundary are marked as hydraulic joint elements (HMJ), and only the broken joint elements connected to the HMJ are subjected to the seepage softening calculation of the discontinuous surface. The stress information of each joint element in the plurality of joint elements in the region of the rock 3D model can be determined by using a stress algorithm matched with the type identification of the corresponding joint element. In the embodiments of the present application, for the joint element whose type is a hydraulic joint element, the stress information of the hydraulic joint element can include the contact force between the two walls of the hydraulic joint element calculated based on the softening process and the first fluid pressure of each node of the hydraulic joint element calculated based on the seepage calculation.

[0048] In some embodiments, in the current time step, all joint elements can be traversed to determine whether the joint element belongs to the fractured joint element, i.e., whether the joint element is damaged (for example, whether the joint element is damaged can be determined by the distance between the plane composed of a1, b1, c1 and the plane composed of a2, b2, c2, such as greater than or equal to a certain threshold, it is considered that the joint element is damaged (i.e., the type of the joint element is updated to the fractured joint element), otherwise it is considered that the joint element is not damaged (i.e., it can be considered that the type of the joint element is still the cohesive joint element)). If the joint element is not damaged, the joint element is not a fractured joint element, and can be marked as a cohesive joint element, and the stress information (such as the stress state) of the joint element in the next time step is calculated using the cohesive transition zone model. As an example, the cohesive transition zone model uses the following formula to calculate the element stress state when calculating the element stress state:

[0049]

[0050] wherein, f D is the shape function, which can be expressed as D is the dimensionless damage parameter, a b and c are the fitting parameters of the empirical curve; C φ are the cohesion and internal friction angle of the rock, respectively.

[0051] In the above embodiment, if the joint element is damaged, it can be marked as a fractured joint element, and the next time step is solved using the discrete contact algorithm. However, before calculation, it is necessary to determine whether the fractured joint element is connected with the HMJ. In one possible implementation, the optional implementation of determining whether the fractured joint element is connected with the HMJ is as follows: the newly generated crack set in the current time step is marked as NJ. Select a crack i in NJ and mark it as FJ(i), i = 1, 2, 3, …, as long as one of the nodes of FJ(i) belongs to HMJ, FJ(i) is deleted from NJ and merged into HMJ. Repeat the above process until all FJ(i) in NJ that are connected with HMJ are deleted from NJ and merged into HMJ. For the above fractured joint element, the discrete contact algorithm can be used to solve the contact force (i.e., the stress information) of the two walls of the fractured joint element, i.e., it can include the normal repulsive force and the tangential friction force.

[0052] In one possible implementation, for the normal repulsive force, the distributed contact force penalty function method can be used to calculate the normal repulsive force between the contact elements (i.e., the contact pairs). This method assumes that the contact pairs can penetrate each other, by calculating the size and shape of the overlapping area, and introducing a normal contact stiffness​​​​p n i.e. the contact force can be calculated.

[0053] In one possible implementation, for the tangential friction force, the Coulomb friction law can be adopted for calculation. It should be noted that in the explicit integration framework, the friction force not only depends on the normal stress on the edge of the interaction element and the discontinuity surface characteristics (roughness and wall strength), but also depends on the relative displacement or velocity between the surfaces. Therefore, a friction model considering the characteristics of the discontinuity surface (such as Figure 4 ) can be adopted, and the tangential friction force can be calculated by the following formula:

[0054]

[0055] In the formula, h is the edge length of the tetrahedral element; is the tangential contact stiffness; is the peak shear strength, for the discontinuity surface of the three-dimensional rock mass, for example, a three-dimensional shear strength theoretical model considering the statistical roughness parameter, i.e. , wherein, is the normal stress on the interaction surface, is the tensile strength parameter of the intact rock, is the basic friction angle, is the maximum apparent inclination angle of the surface; C is the dimensionless fitting parameter; is the normalized area of the joint surface in the shear direction (i.e. the normalized area in the analysis direction is greater than zero); B is the second fitting parameter, which depends on the spatial resolution of the surface digitization; is the three-dimensional roughness in the shear direction. It should be emphasized that other optimized three-dimensional shear strength theoretical formulas are also within the protection scope of the present application; is the residual shear strength, and its expression is , wherein, is the residual friction angle. The sign of is the same as the change amount of the sliding distance on the interaction surface (also referred to as the slip distance) , wherein can be calculated by the following formula:

[0056]

[0057] In the formula, is the change amount of the sliding distance with the time step ; is the relative velocity of the interaction edge. The peak shear strength corresponds to the sliding distance , which can be expressed as:

[0058]

[0059] It is worth noting that for the fractured joint element of the non-hydraulic joint element, the normal repulsive force can be solved by using the above-mentioned distributed contact force penalty function method at each time step, and the tangential friction force can be solved by using the above-mentioned formula (3) to formula (5). In some embodiments, for the hydraulic joint element, the normal repulsive forces of the elements on both sides at each time step are calculated by using the above-mentioned distributed contact force penalty function method, and the tangential friction force needs to be introduced into the interface softening coefficient SC ( t ) and the friction angle decay function on the basis of the above-mentioned tangential friction force calculation method.

[0060] For the tangential friction force of the hydraulic joint element, at each time step, the softening process calculation of the hydraulic joint element can be performed to determine the tangential friction force acting on the two walls of the hydraulic joint element. In a possible implementation, the optional implementation of the above-mentioned softening process calculation of the hydraulic joint element to determine the tangential friction force acting on the two walls of the hydraulic joint element is as follows: obtaining the interface softening coefficient and the friction angle decay function corresponding to the rock; wherein the interface softening coefficient is a decay function of the parameter measuring the strength of the discontinuous face wall with the immersion time, and the friction angle decay function is a decay function of the friction angle with the immersion time; according to the interface softening coefficient, the basic friction angle is modified to obtain the modified basic friction angle; according to the friction angle decay function, the rock tensile strength parameter is modified to obtain the modified rock tensile strength parameter; according to the modified basic friction angle and the modified rock tensile strength parameter, the modified peak shear strength is determined; and according to the modified peak shear strength, the tangential friction force is determined.

[0061] For example, the optional implementation of the above-mentioned determination of the tangential friction force according to the modified peak shear strength is as follows: determining the sliding distance change amount of the interaction edge of the hydraulic joint element and the sliding distance corresponding to the peak shear strength; in the case that the absolute value of the sliding distance change amount is less than or equal to the absolute value of the sliding distance, the tangential friction force is determined according to the modified peak shear strength and the sliding distance change amount; in the case that the absolute value of the sliding distance change amount is greater than the sliding distance, the residual friction angle is modified according to the friction angle decay function, and the modified residual shear strength is determined according to the modified residual friction angle; and the tangential friction force is determined according to the modified peak shear strength, the sliding distance change amount and the modified residual shear strength.

[0062] It is worth noting that, in the embodiments of the present application, for the tangential friction force of the hydraulic joint unit, at each time step, the tangential friction force needs to be modified based on the tangential friction force calculation method (i.e., the above-mentioned formula (3) - formula (5)) by introducing an interface softening coefficient SC ( t ) and a friction angle decay function to the relevant parameters, i.e.

[0063]

[0064] In the formula, is the decay function of the friction angle obtained by the test with the immersion time t ; The fitting relationship between the uniaxial compression strength and the immersion time obtained in the actual test can be used, which can be a decay function of the parameter for measuring the strength of the discontinuous face wall surface with the immersion time, or wherein, D , r and t 0 are fitting parameters related to the cementation type of the rock. It should be noted that the immersion time has little effect on the change of the geometric appearance of the rock joint, so the statistical roughness parameters describing the joint appearance in formula (6) remain unchanged, and only the strength parameters and the friction mechanics parameters and change. In the specific calculation process, when a certain fractured joint unit is marked as HMJ at time step t i , then at the subsequent time step t, the rock tensile strength, friction mechanics parameters are , and .

[0065] When performing the above softening process calculation, seepage calculation also needs to be performed in the hydraulic joint unit. That is, the hydraulic joint unit can be subjected to seepage calculation to determine the first fluid pressure acting on each node of the hydraulic joint unit. In some embodiments, the optional implementation of the above seepage calculation of the hydraulic joint unit to determine the first fluid pressure acting on each node of the hydraulic joint unit is as follows: for any node cavity in each node of the hydraulic joint unit, the first fluid pressure of any node cavity at the current time step is determined according to the second fluid pressure of any node cavity at the previous time step. For example, in the case where the second fluid pressure is greater than the first threshold value, the first fluid pressure can be determined according to the second fluid pressure and the flow sum of all hydraulic joint units to which any node cavity belongs.

[0066] In some embodiments, in case the second fluid pressure is less than the first threshold value, the second fluid pressure is assigned to the first threshold value, and the second node saturation of any node cavity at the current time step is determined according to the first node saturation of any node cavity at the previous time step; wherein the second node saturation is used to determine the flux of the hydraulic-joint element of which the node cavity is the starting node in the seepage direction. For example, in case the first node saturation is less than the second threshold value, the second node saturation is determined according to the first node saturation and the flux sum.

[0067] In some embodiments, in case the second node saturation is greater than the second threshold value, the second node saturation is assigned to the second threshold value, and the fluid pressure of any node cavity at the next time step is determined according to the fluid pressure of any node cavity at the current time step and the flux sum of all hydraulic-joint elements to which the node cavity belongs.

[0068] The seepage calculation will be described in detail below in connection with Figure 5 In the embodiments of the present application, as shown in Figure 5 , the three-dimensional discontinuous surface seepage can be simplified as a two-dimensional plane seepage problem. The total pressure of each node cavity can be respectively represented as:

[0069]

[0070] wherein, , and are the fluid pressures in the node cavities 1, 2 and 3 respectively; is the density of the fluid; g is the gravitational acceleration in the direction of z (z is assumed to be the negative direction); , and are the elevations of the node cavities 1, 2 and 3 respectively.

[0071] Then, a new local two-dimensional coordinate system Figure 5 is established on the seepage surface as shown in x 0 y , the component of the pressure gradient can be represented by the divergence theorem as:

[0072]

[0073] wherein, S is the surface area of the seepage channel, and are the x and y components of the outer normal. Assuming that the fluid pressure between the node cavities varies linearly, the equations (10) and (11) can be approximated as:

[0074]

[0075] where, , and are the fluid pressures at the midpoints of the three edges of the seepage triangle.

[0076] Next, assuming that Darcy's law holds, the fluxes in the directions of the seepage plane x and y can be expressed as:

[0077]

[0078] where, μ is the viscosity of the fluid; a is the average aperture of the hydraulic joint element, where, , and are the hydraulic apertures of the node cavities 1, 2 and 3, respectively. It is noted that when the fluid pressures at the three node cavities of the triangular seepage plane are zero, the fluxes calculated from Equations (14) and (15) are not zero due to gravity, which is not reasonable because the flow rate should decrease with the decrease of the saturation, so that the fluid will not flow out of the zero saturation region. Therefore, an empirical function of saturation f s can be expressed as where, is the saturation of the node at the time step t .

[0079] Next, the fluxes of the node cavities 1, 2 and 3 can be calculated using the following equations:

[0080]

[0081] where, and ( i = AB , BC and CA ) are the components of the outward normal of the three edges of the seepage triangle shown in FIG. 1 in the directions of Figure 5 and x and y . For example, taking the node cavity 1 as an example, Figure 5 the total flux of the node cavity 1 For example, if a cavity at a node is connected to two hydraulic joint units (such as hydraulic joint unit 1 and hydraulic joint unit 2), then the sum of the flow rates (or total flow rates) of the cavity at that node = the sum of the flow rates of hydraulic joint unit 1 + the sum of the flow rates of hydraulic joint unit 2.

[0082] For any node cavity, the fluid pressure is calculated using equation (19):

[0083]

[0084] In the formula, and These are the fluid pressures at the nodes at time steps t and t-1, respectively. It is the bulk modulus of the fluid; It is the time step interval; , ,in, and These represent the volumes of the cavity at time steps t and t-1, respectively. t Time node cavity volume For example, it can be derived from , j It is the index of the seepage channel connected to the cavity of this node, where, .if Then, use equation (19) to calculate. ;if If the fluid is insufficient to fill the node cavity, then Then, the saturation of the node cavity is updated using equation (20):

[0085]

[0086] In the formula, and These are the saturation levels of the node cavity at time steps t and t-1, respectively. If ,but At time step t+1, the saturation of the nodal cavity is calculated using equation (20), but not using equation (19). ;if ,but At time step t+1, the fluid pressure in the node cavity is calculated using equation (19).

[0087] In step 104, based on the stress information of each joint element, the total equivalent nodal force on each node in the multiple joint elements is determined.

[0088] In embodiments of this application, for example at the current time step, the total equivalent nodal force on each node in multiple joint elements of the rock 3D model region can be determined based on the stress information of each joint element. For example, for each joint element, the stress information of that joint element can be converted into the equivalent nodal force of each node within that joint element. Since some nodes may be shared by multiple joint elements (i.e., some nodes may belong to multiple joint elements), for such shared nodes, the equivalent nodal force of that node within each joint element can be calculated first, and then all the equivalent nodal forces of that node can be summed to obtain the total equivalent nodal force of that node. For example, taking a node associated with only one joint element (of type bonded joint element) as an example, the stress state of that bonded joint element can be converted into the equivalent nodal force of that node, which is the total equivalent nodal force of that node. For example, taking a node shared by joint element 1 (a bonded joint element), joint element 2 (a fractured joint element), and joint element 3 (a hydraulic joint element) as an example, the stress state of joint element 1 can be converted into the equivalent nodal force N1 of that node; the contact force (including normal repulsion and tangential friction) of joint element 2 can be converted into the equivalent nodal force N2 of that node; and the contact force (including normal repulsion and tangential friction) of joint element 3 and the first fluid pressure ( If the forces are converted into equivalent nodal forces N3 and N4 respectively, then the total equivalent nodal force of the node = equivalent nodal force N1 + equivalent nodal force N2 + equivalent nodal force N3 + equivalent nodal force N4.

[0089] As an example, the fluid pressure in each nodal cavity of a hydraulic joint element acts as surface pressure on the two walls of the hydraulic joint element, and the fluid pressure can be converted into equivalent nodal forces using equations (21) and (22):

[0090]

[0091] In the formula, p It is the average fluid pressure on the two walls of the hydraulic joint unit, and its value is ;( ), ( ), ( ), ( ), ( )and( They are respectively Figure 6 The spatial coordinates of the six nodes of the hydraulic joint unit (N1N2N3N4N5N6) are shown.

[0092] In step 105, based on the total equivalent node force on each node in the plurality of joint elements, the type identifier of each joint element in the plurality of joint elements is updated, and the total equivalent node force of each node in the next time step is updated to update the type identifier of each joint element.

[0093] For any node in the plurality of joint elements of the region of the three-dimensional model of rock, the acceleration can be determined according to the total equivalent node force of the node and the mass of the node, the velocity can be determined according to the acceleration and the time step interval, the displacement amount can be determined according to the velocity and the time step interval, the coordinate of the node is updated according to the displacement amount to obtain the updated coordinate of the node, and the type identifier of the corresponding joint element is updated according to the updated coordinate of the node (for example, the type identifier of the joint element associated with the hydraulic element is updated). The mass of the node can be equal to one-third of the mass of a triangular pyramid solid element.

[0094] For example, the updated coordinates of other nodes in the associated joint element can be determined in the same way, and the type identifier of the joint element associated with the hydraulic element is updated according to the updated coordinates of the nodes in the associated joint element.

[0095] For example, if the associated joint element includes a cohesive joint element, whether the cohesive joint element is damaged can be determined according to the updated coordinates of the nodes in the cohesive joint element, and if the cohesive joint element is damaged, the type identifier of the cohesive joint element can be updated to the type identifier of a fractured joint element. Then, for the fractured joint element, it can be determined whether the fractured joint element is a hydraulic joint element, and if so, the type identifier of the fractured joint element is updated to the type identifier of the hydraulic joint element, and then the softening process calculation and seepage calculation of the hydraulic joint element can be performed.

[0096] For example, if the normal distance between the midpoints of the two edges of the cohesive joint element is greater than the normal distance threshold, it can be considered that the cohesive joint element is damaged, so that the cohesive joint element is determined to be a fractured joint element.

[0097] For example, if the lateral slip distance of the cohesive joint element is greater than the lateral slip distance threshold, it can be considered that the cohesive joint element is damaged, so that the cohesive joint element is determined to be a fractured joint element.

[0098] For example, if one-half of the sum of the normal distance and the square of the lateral slip distance of the cohesive joint element is outside the preset envelope region, it can be considered that the cohesive joint element is damaged, so that the cohesive joint element is determined to be a fractured joint element.

[0099] For example, whether the fractured joint element is a hydraulic joint element can be determined by the following method: whether the fractured joint element is connected with the hydraulic boundary can be determined, if the fractured joint element is connected with the hydraulic boundary, it can be determined that the fractured joint element becomes a hydraulic joint element.

[0100] For example, whether any node in the fractured joint element belongs to any hydraulic joint element can be determined, if any node in the fractured joint element belongs to any hydraulic joint element, it can be considered that the fractured joint element is connected with the hydraulic boundary, and it can be determined that the fractured joint element is a hydraulic joint element.

[0101] It can be understood that for each node of the joint element in the three-dimensional model region, the acceleration can be determined according to the sum of all equivalent node forces acting on the node and the mass of the node, the velocity can be determined according to the acceleration and the time step interval, the displacement amount can be determined according to the velocity and the time step interval, and the coordinates of the node can be updated according to the displacement amount, so that for the cohesive joint element in the three-dimensional model region, whether the cohesive joint element is damaged can be determined according to the updated coordinates of the nodes of the cohesive joint element, if the cohesive joint element is damaged, the type identifier of the cohesive joint element can be updated to the type identifier of the fractured joint element. Then for the fractured joint element, whether the fractured joint element is a hydraulic joint element can be determined, if yes, the type identifier of the fractured joint element is updated to the type identifier of the hydraulic joint element, and the softening process calculation and the seepage calculation of the hydraulic joint element are performed.

[0102] For example, the type identifier of the cohesive joint element is 0, the type identifier of the fractured joint element is 1, and the type identifier of the hydraulic joint element is marked as 2, if a cohesive joint element is damaged, that is, becomes a fractured joint element, the type identifier of the joint element is updated from 0 to 1, if the fractured joint element is a hydraulic joint element, the type identifier of the joint element can be updated from 1 to 2.

[0103] The above calculation steps of stress information of the joint element, the calculation steps of the total equivalent node force of the node, and the type identifier updating steps of the joint element are repeated for each time step until the calculation time step reaches the target time step, that is, the three-dimensional numerical simulation of the seepage softening of the rock mass discontinuous surface is completed.

[0104] In the embodiments of the present application, the time-dependent seepage-softening model considering the geometric and strength characteristics of the discontinuous surface can intuitively and ingeniously reflect the full-process shear mechanical properties of the discontinuous surface in the tangential sliding of the rock mass, can improve the accuracy of the total equivalent nodal force acting on the nodes of the joint element, improve the prediction accuracy of the rock mass disaster, and then effective prevention and control strategies can be taken, which is suitable for the fluid-structure coupling calculation of the discontinuous surface of various rocks, and improves the robustness and applicability. In addition, when the softening process of the hydraulic joint element is calculated, the decay of the wall strength and the decay of the friction angle of the discontinuous surface with the immersion time are considered, the corrected peak shear strength and the corrected residual friction angle are obtained based on the interface softening function and the friction angle decay function, and the tangential friction force of the hydraulic joint element is calculated based on the corrected peak shear strength and the corrected residual friction angle, which further improves the accuracy of the tangential friction force and the accuracy of the equivalent nodal force converted from the tangential friction force.

[0105] Figure 7 A block diagram of a three-dimensional numerical simulation system for seepage softening of a discontinuous surface of a rock mass is provided in the embodiments of the present application. As shown in the figure, Figure 7 the three-dimensional numerical simulation system for seepage softening of the discontinuous surface of the rock mass can include a preprocessing module 701, a node topology updating module 702, a calculation module 703, and a type updating module 704.

[0106] The preprocessing module 701 is configured to generate a surface domain in the three-dimensional rock model region to represent various discontinuous surfaces in the rock mass as a grid surface, and divide the grid surface into tetrahedral solid element grids to obtain tetrahedral solid elements in the three-dimensional rock model region.

[0107] The node topology updating module 702 is configured to update the node topology of the tetrahedral solid elements to obtain a plurality of joint elements in the three-dimensional rock model region.

[0108] The calculation module 703 is configured to, at an initial time step, determine the stress information of each joint element according to the type identifier of each joint element in the plurality of joint elements by using a stress algorithm matched with the type identifier; wherein the type of the joint element is a hydraulic joint element, and the stress information of the hydraulic joint element includes the contact force of the two walls of the hydraulic joint element calculated based on the softening process and the first fluid pressure of each node of the hydraulic joint element calculated based on the seepage.

[0109] In the embodiments of the present application, the calculation module 703 is further configured to determine the total equivalent nodal force on each node in the plurality of joint elements based on the stress information of each joint element.

[0110] The type updating module 704 is configured to update the type identification of each joint element in the plurality of joint elements based on the total equivalent node force on each node in the plurality of joint elements, and to update the type identification of each joint element in the next time step based on the total equivalent node force on each node in the next time step.

[0111] In some embodiments, the contact force includes a tangential friction force. The calculation module 703 is configured to: obtain an interface softening coefficient and a friction angle decay function corresponding to the rock; the interface softening coefficient is a decay function of a parameter measuring the strength of the discontinuous surface wall with the immersion time, and the friction angle decay function is a decay function of the friction angle with the immersion time; correct the basic friction angle according to the interface softening coefficient to obtain a corrected basic friction angle; correct the rock tensile strength parameter according to the friction angle decay function to obtain a corrected rock tensile strength parameter; determine the corrected peak shear strength according to the corrected basic friction angle and the corrected rock tensile strength parameter; and determine the tangential friction force according to the corrected peak shear strength.

[0112] In some embodiments, the calculation module 703 is configured to: determine a sliding distance change amount and a sliding distance corresponding to the peak shear strength on the interaction edge of the hydraulic joint element; in a case where the absolute value of the sliding distance change amount is less than or equal to the absolute value of the sliding distance, determine the tangential friction force according to the corrected peak shear strength and the sliding distance change amount; in a case where the absolute value of the sliding distance change amount is greater than the sliding distance, correct the residual friction angle according to the friction angle decay function, and determine the corrected residual shear strength according to the corrected residual friction angle; and determine the tangential friction force according to the corrected peak shear strength, the sliding distance change amount, and the corrected residual shear strength.

[0113] In some embodiments, the calculation module 703 is configured to: for any node cavity in each node of the hydraulic joint element, determine the first fluid pressure of the any node cavity in the current time step according to the second fluid pressure of the any node cavity in the previous time step.

[0114] In some embodiments, the calculation module 703 is configured to: in a case where the second fluid pressure is greater than the first threshold value, determine the first fluid pressure according to the second fluid pressure and the flow sum of all hydraulic joint elements to which the any node cavity belongs.

[0115] In some embodiments, the calculating module 703 is further configured to, in a case that the second fluid pressure is less than the first threshold value, assign the second fluid pressure as the first threshold value, and determine the second node saturation of any node cavity at the current time step according to the first node saturation of any node cavity at the previous time step; wherein the second node saturation is used to determine the flow of the hydraulic-joint unit taking any node cavity as a starting node of a seepage direction. For example, the calculating module 703 is configured to, in a case that the first node saturation is less than the second threshold value, determine the second node saturation according to the first node saturation and the flow sum.

[0116] In some embodiments, the calculating module 703 is further configured to, in a case that the second node saturation is greater than the second threshold value, assign the second node saturation as the second threshold value, and determine the fluid pressure of any node cavity at the next time step according to the fluid pressure of any node cavity at the current time step and the flow sum of all hydraulic-joint units to which any node cavity belongs.

[0117] It should be noted that the above-mentioned explanation and description of the embodiments of the three-dimensional numerical simulation method for seepage softening of rock mass discontinuities are also applicable to the three-dimensional numerical simulation system for seepage softening of rock mass discontinuities of the embodiments, which will not be repeated here.

[0118] In order to achieve the above-mentioned embodiments, the present application further provides an electronic device, comprising: at least one processor; and a memory connected in communication with the at least one processor; wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to execute the method provided by the above-mentioned embodiments.

[0119] In order to achieve the above-mentioned embodiments, the present application further provides a storage medium, which stores instructions, and when the instructions run on an electronic device, the electronic device executes the method provided by the above-mentioned embodiments.

[0120] In order to achieve the above-mentioned embodiments, the present application further provides a program product, which comprises at least one of a program and instructions, and when the at least one of the program and instructions is executed by a processor, the steps of the method provided by the above-mentioned embodiments are implemented.

[0121] The collection, storage, use, processing, transmission, provision and disclosure of user personal information involved in the present application comply with relevant laws and regulations and do not violate public order and good customs.

[0122] In the foregoing detailed description, reference is made to descriptive terms such as "one embodiment", "some embodiments", "an example", "a specific example" or "some examples" etc. for describing various embodiments of the application. These descriptive terms are used for the purpose of the description and are not meant to limit or restrict the scope of the application. The use of these terms does not imply that the described feature is essential to the application or that it will be combined with other features as described to implement the application. The scope of the application is defined only by the claims.

[0123] Further, the terms "first", "second", etc. are used herein only to describe various embodiments and do not imply either a relative importance or a specific order of the features being described. Thus, a feature defined with "first" or "second" can implicitly or explicitly include at least one of the feature. The meaning of "a plurality" is at least two, for example, two, three, etc., unless specifically defined otherwise.

[0124] Any process or method descriptions or blocks in flow charts or otherwise described herein represent embodiments of the application that can be managed as one or more modules, segments, or portions of code that include one or more steps for implementing specific logic functions or steps, and the various embodiments of the application can include additional or fewer steps performing the same or equivalent functions. In some embodiments, the blocks can be combined into a software or firmware routine.

[0125] The logic and / or steps represented in flow diagrams or otherwise described herein, for example, can be considered as a sequence of instructions to implement logic functions, and can be embodied in any computer-readable medium for use by an instruction execution system, apparatus, or device, such as a computer-based system, processor- containing system, or other system that can fetch the instructions from the instruction execution system, apparatus, or device and execute the instructions. In the context of this specification, a "computer-readable medium" can be any means that can contain, store, communicate, propagate or transport the program for use by or in connection with the instruction execution system, apparatus, or device. The computer-readable medium can be a machine-readable storage device (e.g., magnetic, optical or other) a machine-readable storage diskette (e.g., floppy disk, optical disk, etc.), a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), optical fibers, and a portable compact disc read-only memory (CDROM). Further, the computer-readable medium can even be paper or other suitable medium upon which the program is printed, as the program can be electronically captured, for example via the optical scanner of a device or other electronic capture device, and then compiled, interpreted, or otherwise processed in a suitable manner, if necessary, and stored in a computer memory.

[0126] It should be understood that aspects of the application can be implemented in hardware, software, firmware or combinations thereof. In the above embodiments, various steps or methods can be implemented in software or firmware that is stored in memory and executed by a suitable instruction execution system. As such, in some embodiments, the steps or methods can be implemented in a combination of hardware and software. If implemented in hardware, as in another embodiment, any of the above techniques can be implemented with or without the use of the following technologies, which technologies are well known in the art: discrete logic circuitry having logic gates for implementing logic functions upon an application of data signals, application specific integrated circuits having appropriate combinational logic gates, programmable gate arrays (PGA), field programmable gate arrays (FPGA), and the like.

[0127] Those of skill in the art would understand that the steps or methods carried out in the above-described embodiments can be carried out by program instructions executed by a processor, and that the program instructions can be stored in a computer readable storage medium. The program instructions, when executed by the processor, can cause the processor to carry out the steps or methods of the embodiments.

[0128] In addition, each of the functional units in the various embodiments of the present application can be integrated in one processing module, or each of the units can be physically present separately, or two or more units can be integrated in one module. The integrated module can be implemented in the form of hardware or in the form of a software functional module. When the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it can also be stored in a computer readable storage medium.

[0129] The storage medium mentioned above can be a read-only memory, a magnetic disk or an optical disk, etc. Although the embodiments of the present application have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present application, and those skilled in the art can make changes, modifications, replacements and variations to the above embodiments within the scope of the present application.

Claims

1. A three-dimensional numerical simulation method for seepage softening of rock mass discontinuities, characterized in that, The method comprises the following steps: generating a surface domain used to represent various types of discontinuous surfaces in a rock mass as a grid surface in a rock three-dimensional model region, and dividing the grid surface into a tetrahedral element grid to obtain a tetrahedral element of the rock three-dimensional model region; updating the node topology of the tetrahedral element to obtain a plurality of joint elements in the rock three-dimensional model region; at a current time step, determining the stress information of each joint element according to the type identification of each joint element in the plurality of joint elements by using a stress algorithm matched with the type identification; wherein the type of the joint element is a hydraulic joint element, the stress information of the hydraulic joint element includes the contact force of the two walls of the hydraulic joint element calculated based on a softening process and the first fluid pressure of each node of the hydraulic joint element calculated based on seepage; determining the total equivalent node force on each node in the plurality of joint elements based on the stress information of each joint element; updating the type identification of each joint element in the plurality of joint elements based on the total equivalent node force on each node in the plurality of joint elements, and updating the total equivalent node force of each node in the next time step to update the type identification of each joint element; wherein the type of the joint element includes a cohesive joint element, a fracture joint element and the hydraulic joint element; the type identification updating mode of the joint element includes updating the type identification of the cohesive joint element to the type identification of the fracture joint element, and updating the type identification of the fracture joint element to the type identification of the hydraulic joint element.

2. The method of claim 1, wherein, The contact force includes a tangential friction force; the softening process calculation of the hydraulic joint element to determine the contact force acting on the two walls of the hydraulic joint element comprises the following steps: obtaining an interface softening coefficient and a friction angle attenuation function corresponding to the rock; wherein the interface softening coefficient is an attenuation function of the parameter measuring the strength of the wall surface of the discontinuous surface with the immersion time, and the friction angle attenuation function is an attenuation function of the friction angle with the immersion time; correcting the basic friction angle according to the interface softening coefficient to obtain a corrected basic friction angle; correcting the rock tensile strength parameter according to the friction angle attenuation function to obtain a corrected rock tensile strength parameter; determining a corrected peak shear strength according to the corrected basic friction angle and the corrected rock tensile strength parameter; determining the tangential friction force according to the corrected peak shear strength.

3. The method of claim 2, wherein, The determination of the tangential friction force according to the corrected peak shear strength comprises the following steps: determining the sliding distance change amount of the interaction edge of the hydraulic joint element and the sliding distance corresponding to the peak shear strength; in the case that the absolute value of the sliding distance change amount is less than or equal to the absolute value of the sliding distance, determining the tangential friction force according to the corrected peak shear strength and the sliding distance change amount; In a case where an absolute value of the sliding distance variable is greater than the sliding distance, a residual friction angle is corrected according to the friction angle attenuation function, and a corrected residual shear strength is determined according to the corrected residual friction angle; The tangential friction force is determined according to the corrected peak shear strength, the sliding distance variable and the corrected residual shear strength.

4. The method of claim 1, wherein, The permeation calculation of the hydraulic joint unit is performed to determine the first fluid pressure acting on each node of the hydraulic joint unit, including: For any node cavity in each node of the hydraulic joint unit, the first fluid pressure of the any node cavity at a current time step is determined according to the second fluid pressure of the any node cavity at a previous time step.

5. The method of claim 4, wherein, The first fluid pressure of the any node cavity at the current time step is determined according to the second fluid pressure of the any node cavity at the previous time step, including: In a case where the second fluid pressure is greater than a first threshold value, the first fluid pressure is determined according to the second fluid pressure and a flow sum of all hydraulic joint units to which the any node cavity belongs.

6. The method of claim 5, wherein, The method further includes: In a case where the second fluid pressure is less than the first threshold value, the second fluid pressure is assigned as the first threshold value, and a second node saturation of the any node cavity at the current time step is determined according to a first node saturation of the any node cavity at the previous time step; The second node saturation is used to determine a flow of a hydraulic joint unit having the any node cavity as a starting node in a seepage direction.

7. The method of claim 6, wherein, The second node saturation of the any node cavity at the current time step is determined according to the first node saturation of the any node cavity at the previous time step, including: In a case where the first node saturation is less than a second threshold value, the second node saturation is determined according to the first node saturation and the flow sum.

8. The method according to claim 6 or 7, characterized in that, The method further includes: In a case where the second node saturation is greater than the second threshold value, the second node saturation is assigned as the second threshold value, and a fluid pressure of the any node cavity at a next time step is determined according to the fluid pressure of the any node cavity at the current time step and the flow sum of all hydraulic joint units to which the any node cavity belongs.

9. A three-dimensional numerical simulation system for seepage softening of rock mass discontinuities, characterized in that, Including: A pretreatment module is configured to generate a surface domain representing various discontinuous surfaces in a rock mass in a rock three-dimensional model region as a grid surface, and divide the grid surface into a tetrahedral element grid to obtain a tetrahedral element of the rock three-dimensional model region; A node topology updating module is configured to update a node topology of the tetrahedral element to obtain a plurality of joint units in the rock three-dimensional model region; A node topology updating module is configured to update a node topology of the tetrahedral element to obtain a plurality of joint units in the rock three-dimensional model region; The computing module is configured to, at a current time step, determine stress information of each of the joint elements according to a type identifier of each of the joint elements, by using a stress algorithm matched with the type identifier; wherein the type of the joint element is a hydraulic joint element, and the stress information of the hydraulic joint element includes contact force of two walls of the hydraulic joint element calculated based on a softening process and first fluid pressure of each node of the hydraulic joint element calculated based on seepage calculation; The computing module is further configured to determine total equivalent node force on each node in the plurality of joint elements based on the stress information of each of the joint elements; The type updating module is configured to update the type identifier of each of the joint elements in the plurality of joint elements based on the total equivalent node force on each node in the plurality of joint elements, and to update the type identifier of each of the joint elements in the next time step based on the total equivalent node force on each node in the next time step. The type of the joint element includes a cohesive joint element, a fracture joint element and the hydraulic joint element, and the type identifier updating manner of the joint element includes updating the type identifier of the cohesive joint element to the type identifier of the fracture joint element, and updating the type identifier of the fracture joint element to the type identifier of the hydraulic joint element.

10. An electronic device, comprising: The method comprises the following steps: at least one processor; and a memory connected with the at least one processor in communication; wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the method in any one of claims 1-8.

Images