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

Abstract

The invention provides a three-dimensional numerical simulation method and system for seepage softening of a discontinuous surface of a rock mass. The method comprises the following steps: generating a surface domain used for representing various discontinuous surfaces in a rock mass in a rock three-dimensional model region as a grid surface, and dividing the grid surface into tetrahedron entity unit grids to obtain tetrahedron entity units of the rock three-dimensional model region; updating a node topological structure of the tetrahedral entity unit to obtain a plurality of joint units in the rock three-dimensional model region; at the initial time step, according to the type identifier of each joint unit, determining stress information of each joint unit by adopting a stress algorithm matched with the type identifier; based on the stress information of each joint unit, determining a total equivalent node force on each node in the plurality of joint units; and based on the total equivalent node force, updating the type identifier of each joint unit in the plurality of joint units, and performing total equivalent node force of each node in the next time step so as to update the type identifier of each joint unit.

Description

Technical Field

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

[0002] Discontinuity surface systems such as fractures, bedding planes, and joints are widely developed in natural rock masses. Their seepage softening effect not only directly affects engineering stability problems caused by the attenuation of the self-supporting force of rock masses (such as underground cavern collapses and slope instabilities), but is also a key factor restricting the efficiency of energy development such as shale gas reservoir stimulation and deep mining. This softening process is essentially a fluid-solid coupling process of the interaction between the deterioration of the mechanical properties of fracture surfaces and the deformation of the medium driven by fluid pressure.

[0003] However, the existing theoretical models and technologies of fluid-solid coupling of rock masses are not universal, and their completeness urgently needs to be optimized and supplemented. Summary of the Invention

[0004] Embodiments of the present application provide a three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuity surfaces.

[0005] In a first aspect, embodiments of the present application provide a three-dimensional numerical simulation method for seepage softening of rock mass discontinuity surfaces, including: Generating a domain used to represent various discontinuity surfaces in the rock mass as a mesh surface within the three-dimensional model region of the rock, and dividing the mesh surface into tetrahedral solid element meshes to obtain the tetrahedral solid elements in the three-dimensional model region of the rock; Updating the node topological structure of the tetrahedral solid elements to obtain multiple joint elements in the three-dimensional model region of the rock; At an initial time step, according to the type identifier of each joint element among the multiple joint elements, determining the stress information of each joint element by using a stress algorithm matching 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 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 seepage; Based on the stress information of each joint element, determining the total equivalent nodal force on each node among the multiple joint elements; Based on the total equivalent nodal force on each node among the multiple joint elements, updating the type identifier of each joint element among the multiple joint elements, and performing the total equivalent nodal force of each node in the next time step to update the type identifier of each joint element.

[0006] In a second aspect, embodiments of the present application provide a three-dimensional numerical simulation system for seepage softening of rock mass discontinuity surfaces, including: A preprocessing module, configured to generate a region for characterizing various discontinuity surfaces in a rock mass within a three-dimensional rock model region as a mesh surface, and divide the mesh surface into a tetrahedral solid element mesh to obtain the tetrahedral solid elements in the three-dimensional rock model region; A node topological structure updating module, configured to update the node topological structure of the tetrahedral solid elements to obtain a plurality of joint elements in the three-dimensional rock model region; A calculation module, configured to, at an initial time step, according to the type identifier of each joint element in the plurality of joint elements, determine the stress information of each joint element by using a stress algorithm matching 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 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 seepage; The calculation module 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; A type updating module, configured to update the type identifier of each joint element in the plurality of joint elements based on the total equivalent nodal force on each node in the plurality of joint elements, and perform the total equivalent nodal force of each node in the next time step to update the type identifier of each joint element.

[0007] In a third aspect, an embodiment of the present application provides an electronic device, including: At least one processor; and A memory communicatively connected to 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 described in the foregoing first aspect.

[0008] In a fourth aspect, an embodiment of the present application provides a storage medium, where the storage medium stores instructions, and when the instructions run on an electronic device, the electronic device is enabled to execute the method described in the foregoing first aspect.

[0009] In a fifth aspect, an embodiment of the present application provides a program product, where the program product includes at least one of a program and instructions, and when at least one of the program and instructions is executed by a processor, the steps of the method described in the foregoing first aspect are implemented.

[0010] According to the technical solution of the present application, a time-dependent seepage-softening model that takes into account the geometric morphology characteristics and strength characteristics of discontinuity surfaces can ingeniously and intuitively reflect the full-process shear mechanical properties of rock mass discontinuity surfaces during tangential sliding, broadening the applicable range of the rock mass fluid-solid coupling model.

[0011] Additional aspects and advantages of the present application will be given in part in the following description, become apparent in part from the following description, or be understood through the practice of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The above and / or additional aspects and advantages of the present application will become apparent and understandable from the following description of embodiments in conjunction with the accompanying drawings, in which: Figure 1 is a schematic flow chart of a three-dimensional numerical simulation method for seepage softening of rock mass discontinuity surfaces provided by an embodiment of the present application; Figure 2 is a schematic diagram of the effect of tetrahedral mesh generation provided by an embodiment of the present application; Figure 3 is an example diagram of the three-dimensional joint element topological structure provided by an embodiment of the present application; Figure 4 is a schematic diagram of the definition of tangential frictional force provided by an embodiment of the present application; Figure 5 is an example diagram of seepage channels of triangular surfaces in the global and local coordinate systems provided by an embodiment of the present application; Figure 6 is a schematic diagram of the equivalent nodal force for fluid pressure calculation provided by an embodiment of the present application; Figure 7 is a block diagram of a three-dimensional numerical simulation system for seepage softening of rock mass discontinuity surfaces provided by an embodiment of the present application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0013] Embodiments of the present application will be described in detail below. Examples of the embodiments are shown in the accompanying drawings, where the same or similar reference numerals denote 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 should not be construed as limiting the present application.

[0014] It should be noted that discontinuum systems such as fractures, bedding planes, and joints are widely developed in natural rock masses. Their seepage softening effect not only directly affects engineering stability problems caused by the attenuation of the self-supporting force of rock masses (such as underground cavity collapses and slope instabilities), but is also a key factor restricting the efficiency of energy development such as shale gas reservoir stimulation and deep mining. This softening process is essentially a fluid-solid coupling process of the interaction between the deterioration of the mechanical properties of fracture surfaces and the deformation of the medium driven by fluid pressure. With the progress of high-performance computing technology, numerical methods such as the discrete element method (DEM) and discontinuous deformation analysis (DDA) have become important tools for studying this problem. Existing methods usually simplify the discontinuum surface into a parallel plate model, establish the relationship between the aperture and the flow rate based on the cubic law, and achieve fluid-solid coupling calculations by iteratively updating the fluid pressure and the fracture aperture. However, there are still significant gaps in the refined characterization of the seepage softening mechanism, especially in the analysis of the softening effect in a three-dimensional complex environment. Existing numerical models mainly target jointed rock masses with weak water sensitivity. Their core limitations are as follows: (1) only considering the hydrodynamics of fracture surfaces and ignoring the physicochemical softening mechanisms such as mineral swelling and cement weakening that occur in highly water-sensitive rock masses during the seepage process; (2) using two-dimensional or quasi-three-dimensional simplified models and being unable to truly reflect the coupling effect of the three-dimensional seepage field and stress field in complex fracture networks; (3) lacking the ability to dynamically track the localization of rock mass deformation and the propagation of fractures during the softening process. When facing typical highly water-sensitive rock masses such as water-rich fault zones and mudstone interbeds, existing methods can neither accurately describe the law of fracture surface strength attenuation caused by seepage nor simulate the multi-scale fracture evolution process induced by softening, resulting in deviations in engineering stability assessment.

[0015] In view of the above problems, there is an urgent practical need to construct a refined model that can consider the coupling of the softening characteristics of rock mass materials and the three-dimensional seepage field. The core advantages of the three-dimensional seepage softening model are as follows: by introducing softening parameters to characterize the law of fracture surface strength attenuation under seepage action and combining the pore-fracture dual-medium seepage theory to describe the three-dimensional flow field distribution, it can accurately depict the spatial differences in the softening effect of rock masses under different seepage directions and stress states. This model not only breaks through the limitations of the traditional parallel plate hypothesis, but can also reveal the complex mechanical behaviors of fracture opening-sliding-penetration caused by seepage softening under three-dimensional stress states, providing a more reliable theoretical support for the support design of deep underground engineering and the prediction of hydraulic fracture propagation. Therefore, carrying out the theoretical construction and numerical implementation of the three-dimensional seepage softening model is not only a key link to fill the current gap in fluid-solid coupling research, but also an important technical path to promote the transformation of rock mass engineering stability analysis from qualitative assessment to quantitative prediction.

[0016] Fluids vary in their degree of altering the physical and mechanical properties of different rock types. Existing fluid-solid coupling models for rock mass discontinuities do not focus on the modification process of fluids on different types of rocks. Instead, the effect of fluids on rocks is uniformly characterized by applying fluid pressure to the boundary of the rock skeleton. Although the above approach has a certain degree of applicability to some rocks with low or even non-water sensitivity (such as granite, marble, etc.), for cemented rocks (such as mudstone, shale, argillaceous sandstone, etc.) widely developed in some mountainous areas and coal-bearing strata, using existing technologies often results in calculated results being better than actual on-site results, which has a direct negative impact on the prediction accuracy of engineering soft rock mass disasters and the selection of prevention and control strategies. Therefore, the existing fluid-solid coupling theoretical models and technologies for rock masses are not universal, and their completeness urgently needs to be optimized and supplemented. Aiming at the shortcomings of the existing technologies, this application proposes a three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuities. This method can accurately match the corresponding rock types through a unified form of softening equation, thereby improving the robustness and applicability of the fluid-solid coupling model in the fields of rock mechanics and rock engineering.

[0017] The modification effect of fluids on the physical and mechanical properties of rocks has significant lithological differences, while existing fluid-solid coupling models for rock mass discontinuities generally ignore this core feature. Their simplified calculation method of uniformly applying fluid pressure boundary conditions to the rock skeleton has a certain degree of applicability to low-water-sensitivity rocks such as granite and marble, but significant defects are exposed in water-sensitive cemented rocks such as mudstone, shale, and argillaceous sandstone widely distributed in some mountainous areas. During the seepage process, such rocks will undergo physicochemical softening reactions such as clay mineral swelling and cement dissolution, resulting in a non-linear attenuation of mechanical strength. However, the "homogenized" treatment method of existing models often makes the calculated results overly optimistic, directly affecting the prediction accuracy of soft rock engineering disasters (such as large tunnel deformations and water inrush from coal seam floors) and the effectiveness of prevention and control strategies.

[0018] There are three major key gaps in current research: First, a multi-mechanism coupling model of fluid-rock interaction has not been established, lacking differential characterization of softening paths (expansion type, dissolution type, ion exchange type) dominated by different mineral compositions (such as montmorillonite, illite); second, two-dimensional plane or quasi-three-dimensional simplified assumptions are adopted, unable to capture the anisotropic influence of the three-dimensional seepage field (such as seepage direction, velocity gradient, pressure distribution) in complex fracture networks on the softening process; third, a full coupling framework of "seepage-softening-deformation" has not been constructed, making it difficult to describe the dynamic feedback process of softening-induced fracture propagation, stress redistribution, and seepage channel evolution. This theoretical deficiency is particularly prominent in projects such as gas drainage in coal-bearing strata and hydraulic fracturing of shale gas horizontal wells, often leading to the failure of early warnings for accidents such as wellbore instability and fracture cross-layer induced by softening.

[0019] In response to the above challenges, constructing a three-dimensional seepage softening model has irreplaceable engineering value. By introducing a lithology-sensitive softening constitutive equation, the model can accurately characterize the strength attenuation law of different rock types in a three-dimensional seepage environment, and by coupling the discrete fracture network (DFN) modeling technology, it can realize the three-dimensional dynamic simulation of crack initiation, expansion, and penetration during the softening process. This model can not only provide accurate tools for the stability analysis of surrounding rocks of hydropower projects under certain complex geological conditions, but also support the construction of a "precise prediction-targeted prevention and control" technical system for softening disasters in coal-bearing strata, fundamentally solving the problem of insufficient universality of existing models, and promoting the rock mass fluid-solid coupling theory from "simplified simulation" to "real characterization". Specifically, the three-dimensional numerical simulation method and system for seepage softening of rock mass discontinuities in the embodiments of the present application are described below with reference to the accompanying drawings.

[0020] It should be noted that the execution subject of the three-dimensional numerical simulation method of seepage softening of rock mass discontinuity surface in the embodiment of the present application can be a three-dimensional numerical simulation system of seepage softening of rock mass discontinuity surface, which can be implemented by software and / or hardware, and the system can be configured in an electronic device. Exemplarily, the electronic device can include but is not limited to a terminal, a server, etc.

[0021] Figure 1 The following is a flow chart of a three-dimensional numerical simulation method for seepage softening of rock mass discontinuities provided in the embodiment of the present application. Figure 1 As shown, the three-dimensional numerical simulation method for seepage softening of the discontinuity surface of the rock mass may include but is not limited to the following steps.

[0022] In step 101, a surface domain used to characterize various discontinuities in the rock mass is generated as a mesh surface in the rock three-dimensional model region, and the mesh surface is divided into tetrahedral solid unit meshes to obtain tetrahedral solid units in the rock three-dimensional model region.

[0023] For example, Figure 2 As shown, taking a certain type of rock 3D model as an example, a surface domain used to characterize various discontinuities such as cracks, bedding, and joints in the rock mass can be generated in the rock 3D model area, and the surface domains of various discontinuities are used as mesh surfaces and divided into tetrahedral solid unit meshes, so that the tetrahedral solid units in the rock 3D model area can be obtained. The rock type can be any rock type. For example, the rock type can be any rock type among mudstone, shale, argillaceous sandstone, granite, and marble.

[0024] In step 102, the node topology of the tetrahedral solid element is updated to obtain a plurality of joint elements in the rock three-dimensional model region.

[0025] For example, the node topology of the tetrahedral solid elements in the three-dimensional rock model area can be updated so that all tetrahedral solid elements have independent node numbers, and then, as shown in Figure 3 , a six-node joint element (N1N2N3N4N5N6) that can simulate the continuous-discontinuous process can be embedded, and the joint elements are divided into fracture joint elements representing discontinuous surfaces and bonded joint elements representing continuous surfaces. Among them, the bonded joint elements can be transformed into fracture joint elements after calculation. Figure 3 The boundary area between the triangular pyramids in

[0026] can be regarded as a joint element. For example, the boundary area between the surface composed of a1, b1, c1 and the surface composed of a2, b2, c2 can be regarded as a joint element. If these two surfaces are characterized as continuous surfaces, the type of this joint element is a bonded joint element. If the bonded joint element is softened by fluid seepage to form crack initiation, propagation, and penetration, resulting in these two surfaces being characterized as discontinuous surfaces, the type of this joint element is updated from a bonded joint element to a fracture joint element. Exemplarily, the bonded joint elements on the hydraulic boundary can be marked as hydraulic joint elements (HMJ), and only the fracture joint elements connected to the HMJ are subjected to discontinuous surface seepage softening calculation.

[0027] In the embodiment of the present application, at the initial time step, the discontinuous surfaces such as cracks, bedding planes, and joints in the three-dimensional rock model area can be marked as fracture joint elements, and then the continuous surfaces between the remaining tetrahedral solid elements are marked as bonded joint elements; the bonded joint elements on the water boundary are marked as hydraulic joint elements (HMJ), and only the fracture joint elements connected to the HMJ are subjected to discontinuous surface seepage softening calculation. The stress information of the corresponding joint elements can be determined by using a stress algorithm matching the type identifier according to the type identifier of each joint element in the three-dimensional rock model area. Among them, in the embodiment of the present application, for the joint element of the type of 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.

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

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

[0030] Continuing with the above embodiments, if the joint element is damaged, it can be marked as a fractured joint element and solved using the discrete contact algorithm in the next time step. However, before the calculation, it is necessary to determine whether the fractured joint element is connected to the HMJ. In one possible implementation, the optional implementation of determining whether the fractured joint element is connected to the HMJ is as follows: Mark the set of newly generated fractures in the current time step as NJ. Select fracture 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 the HMJ, then delete FJ(i) from NJ and incorporate it into the HMJ. Repeat the above process until all FJ(i) connected to the HMJ in NJ are deleted from NJ and incorporated into the HMJ. For the above fractured joint element, the discrete contact algorithm can be used to solve the contact force between its two walls (that is, the stress information of the fractured joint element), that is, it can include the normal repulsive force and the tangential frictional force.

[0031] 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 contact elements (i.e., contact pairs). This method assumes that the contact pairs can penetrate each other, calculates the size and shape of the overlapping area, and introduces the normal contact stiffness pn , the contact force can be calculated.

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

[0033] In the formula, h is the side length of the tetrahedral element; is the tangential contact stiffness; is the peak shear strength. For the three-dimensional rock mass discontinuity surface, taking the three-dimensional shear strength theoretical model considering the statistical roughness parameter as an example, that is , where 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 dip 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 this application; is the residual shear strength, and its expression is , where is the residual friction angle. The sign of is the same as the sign of the change in the sliding distance (also called the slip distance) on the interaction surface, where

[0034] In the formula, is the change in the sliding distance with respect to the time step ; is the relative velocity of the interaction edge. The peak shear strength corresponding to the sliding distance can be expressed as:

[0035] It should be noted that for the fractured joint unit of the non-hydro-joint unit, the normal repulsive force can be solved by the above-mentioned distributed contact force penalty function method at each time step, and the tangential frictional force can be solved by the above-mentioned formulas (3)-(5). In some embodiments, for the hydro-joint unit, the normal repulsive force between the two side units is calculated by the above-mentioned distributed contact force penalty function method at each time step, while the tangential frictional force needs to introduce an interface softening coefficient on the basis of the above-mentioned tangential frictional force calculation method SC ( t ) and the friction angle decay function , and correct the relevant parameters.

[0036] For the tangential frictional force of the hydro-joint unit, at each time step, the softening process of the hydro-joint unit can be calculated to determine the tangential frictional force acting on the two walls of the hydro-joint unit. In one possible implementation manner, the optional implementation manner for calculating the softening process of the hydro-joint unit and determining the tangential frictional force acting on the two walls of the hydro-joint unit is as follows: obtain 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 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 the corrected basic friction angle; correct the rock tensile strength parameter according to the friction angle decay function to obtain the 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; determine the tangential frictional force according to the corrected peak shear strength.

[0037] Exemplarily, the optional implementation manner for determining the tangential frictional force according to the corrected peak shear strength is as follows: determine the change amount of the sliding distance on the interaction side of the hydro-joint unit and the sliding distance corresponding to the peak shear strength; when the absolute value of the change amount of the sliding distance is less than or equal to the absolute value of the sliding distance, determine the tangential frictional force according to the corrected peak shear strength and the change amount of the sliding distance; when the absolute value of the sliding distance variable 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; determine the tangential frictional force according to the corrected peak shear strength, the change amount of the sliding distance and the corrected residual shear strength.

[0038] It should be noted that in the embodiments of the present application, for the tangential frictional force of the hydro-joint unit, at each time step, the tangential frictional force needs to introduce an interface softening coefficient on the basis of the above-mentioned tangential frictional force calculation method (i.e., the above-mentioned formulas (3)-(5)) SC ( t ) and the friction angle decay function , the relevant parameters are corrected, that is

[0039] In the formula, is the attenuation function of the friction angle obtained from the test with respect to the immersion time t ; The fitting relationship between the uniaxial compressive strength and the immersion time obtained from the actual test can be adopted. It can be the attenuation function of the parameter used to measure the strength of the discontinuity surface wall with respect to the immersion time, or can be adopted, where D , r and t 0 are all fitting parameters and are related to the cementation type of the rock. It should be noted that the change degree of the geometric morphology of the rock joints caused by the immersion time is very small. Therefore, the statistical roughness parameters describing the joint morphology in formula (6) remain unchanged, and only the strength parameter and the friction mechanical parameters and change. In the specific calculation process, when a certain fractured joint unit is marked as HMJ at the time step t i , then at the subsequent time step t, the tensile strength of the rock and the friction mechanical parameters are respectively , and .

[0040] When performing the above softening process calculation, in the hydraulic joint unit, the seepage calculation also needs to be carried out. That is to say, the permeability calculation can be performed on the hydraulic joint unit to determine the first fluid pressure acting on each node of the hydraulic joint unit. In some embodiments, the optional implementation manner of the above-mentioned permeability calculation on 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 among the nodes 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. Exemplarily, when the second fluid pressure is greater than the first threshold, the first fluid pressure can be determined according to the second fluid pressure and the sum of the flow rates of all hydraulic joint units to which any node cavity belongs.

[0041] In some embodiments, when the second fluid pressure is less than the first threshold, the second fluid pressure is assigned the first threshold, and the second node saturation of any node cavity at the current time step is determined based on the first node saturation of the any node cavity at the previous time step; wherein, the second node saturation is used to determine the flow rate of the hydraulic joint element with any node cavity as the starting node of the seepage direction. Exemplarily, when the first node saturation is less than the second threshold, the second node saturation is determined based on the first node saturation and the sum of the flow rates.

[0042] In some embodiments, when the second node saturation is greater than the second threshold, the second node saturation is assigned the second threshold, and the fluid pressure of any node cavity at the next time step is determined based on the fluid pressure of the any node cavity at the current time step and the sum of the flow rates of all the hydraulic joint elements to which the any node cavity belongs.

[0043] The following will be combined with Figure 5 to describe the seepage calculation in detail. In the embodiments of the present application, as Figure 5 shown, the three-dimensional discontinuous surface seepage can be simplified to a two-dimensional plane seepage problem. The total pressure of each node cavity can be respectively expressed as:

[0044] In the formula, , and are the fluid pressures in node cavities 1, 2, and 3 respectively; is the density of the fluid; g is z the gravitational acceleration in the direction (assuming its direction is the negative direction of z); and are the elevations of node cavities 1, 2, and 3 respectively.

[0045] Next, a new local two-dimensional coordinate system Figure 5 0 x is established on the seepage surface as shown in y , then the components of the pressure gradient can be expressed by the divergence theorem as:

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

[0047] In the formula, , and are respectively the fluid pressures at the midpoints of the three sides of the seepage triangular surface.

[0048] Next, assuming that Darcy's law holds and adopting the parallel plate seepage model approximated by the cubic law, the flow rates in the x and y directions of the seepage plane can be expressed as:

[0049] In the formula, μ is the viscosity of the fluid; a is the average aperture of the hydraulic joint element, , where , and are respectively the hydraulic apertures of the node cavities 1, 2, and 3. It should be noted that when the fluid pressures at the three node cavities of the triangular seepage surface are zero, due to the action of gravity, the flow rates calculated according to Eqs. (14) and (15) are not zero, which is unreasonable because the flow velocity should decrease with the decrease of saturation, so that the fluid will not flow out of the zero-saturation region. Therefore, the saturation empirical function f s can be introduced and expressed as , where is the saturation of the node at the time step t .

[0050] Next, the flow rates of the node cavities 1, 2, and 3 can be calculated using the following formula:

[0051] In the formula, and ( i = AB , BC and CA ) are respectively the components of the outer normal vectors of the three sides of the seepage triangular surface shown in Figure 5 in the x and y directions. For example, taking the node cavity 1 as an example, Figure 5 in which there is only one hydraulic joint element connected to it, so the total flow rate of the node cavity 1. Another example is that if there are two hydraulic joint elements connected to a certain node cavity (such as hydraulic joint element 1 and hydraulic joint element 2), then the sum of the flow rates of the node cavity (or the total flow rate) = the sum of the flow rates of hydraulic joint element 1 + the sum of the flow rates of hydraulic joint element 2.

[0052] For any node cavity, its fluid pressure is calculated using Equation (19):

[0053] wherein, and are the fluid pressures of the node at time steps t and t - 1, respectively; is the bulk modulus of the fluid; is the time step interval; , , where, and are the volumes of the above node cavity at time steps t and t - 1, respectively. Taking the volume t of the node cavity at time step as an example, it can be obtained from , j is the seepage channel index connected to the node cavity, where, . If , then Equation (19) is used to calculate ; if , the fluid is not sufficient to fill the node cavity, then , and then Equation (20) is used to update the saturation of the node cavity:

[0054] wherein, and are the saturations of the node cavity at time steps t and t - 1, respectively. If , then , and the saturation of the node cavity is calculated using Equation (20) at time step t + 1, and Equation (19) is not used to calculate ; if , then , and the fluid pressure of the node cavity is calculated using Equation (19) at time step t + 1.

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

[0056] In an embodiment of the present application, for example, at the current time step, the total equivalent nodal force at each node in the multiple joint elements in the three-dimensional rock model area can be determined according to the stress information of each joint element. Exemplarily, for each joint element, the stress information of the joint element can be converted into the equivalent nodal forces of the nodes of the 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 forces of the node in each joint element can be calculated first, and then all the equivalent nodal forces of the node are summed up to obtain the total equivalent nodal force of the node. For example, taking the case where a node is only associated with one joint element (type: bonded joint element) as an example, the stress state of the bonded joint element can be converted into the equivalent nodal force of the node, which is the total equivalent nodal force of the node. Another example is that taking the case where a node is shared by joint element 1 (type: bonded joint element), joint element 2 (type: fractured joint element), and joint element 3 (type: hydraulic joint element) as an example, the stress state of joint element 1 can be converted into the equivalent nodal force N1 of the node, the contact force of joint element 2 (including normal repulsive force and tangential frictional force) can be converted into the equivalent nodal force N2 of the node, and the contact force of joint element 3 (including normal repulsive force and tangential frictional force) and the first fluid pressure ( ) are respectively converted into the equivalent nodal forces N3 and N4 of the node, 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.

[0057] As an example, the fluid pressure in each node cavity of the hydraulic joint element acts on the two walls of the hydraulic joint element as the surface pressure, and the fluid pressure can be converted into the equivalent nodal force through Eqs. (21) and (22):

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

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

[0060] For any node in multiple joint elements in the three-dimensional rock model area, the acceleration can be determined based on the total equivalent nodal force of the node and the mass of the node. Based on the acceleration and the time step interval, the velocity can be determined. Based on the velocity and the time step interval, the displacement amount can be determined. Based on the displacement amount, the coordinates of the node are updated to obtain the updated coordinates of the node. Based on the updated coordinates of the node, the type identifier of the corresponding joint element is updated (such as updating the type identifier of the joint element associated with the hydraulic element). Among them, the mass of the node can be equal to one-third of the mass of a triangular pyramid solid element.

[0061] Exemplarily, in the same way, the updated coordinates of other nodes in the associated joint element can be determined, and based on the updated coordinates of each node in the associated joint element, the type identifier of the joint element associated with the hydraulic element is updated.

[0062] For example, if the associated joint element includes a cohesive joint element, based on the updated coordinates of each node of the cohesive joint element, it is determined whether the cohesive joint element has been damaged. If it has been 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, it can be determined whether the fractured joint element is a hydraulic joint element. 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 can be performed on the hydraulic joint element.

[0063] Exemplarily, if the normal distance between the midpoints of two sides of the cohesive joint element is greater than the normal distance threshold, it can be considered that the cohesive joint element is damaged, and thus the cohesive joint element is determined to be a fractured joint element.

[0064] Exemplarily, if the lateral slip distance of the cohesive joint element is greater than the lateral slip distance, it can be considered that the cohesive joint element is damaged, and thus the cohesive joint element is determined to become a fractured joint element.

[0065] Exemplarily, if the square root of the sum of the squares of the normal distance and the lateral slip distance of the cohesive joint element is outside the preset envelope area, it can be considered that the cohesive joint element is damaged, and thus the cohesive joint element is determined to become a fractured joint element.

[0066] Exemplarily, the following method can be used to determine whether a fractured joint element is a hydraulic joint element: It can be determined whether the fractured joint element is connected to the hydraulic boundary. If the fractured joint element is connected to the hydraulic boundary, it can be determined that the fractured joint element becomes a hydraulic joint element.

[0067] As an example, it can be determined whether any node in the fractured joint element belongs to any hydraulic joint element. 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 to the hydraulic boundary, and it can be determined that the fractured joint element is a hydraulic joint element.

[0068] It can be understood that for each node of the joint element in the three-dimensional model area, the acceleration can be determined according to the sum of all equivalent nodal forces acting on the node and the mass of the node. According to the acceleration and the time step interval, the velocity can be determined. According to the velocity and the time step interval, the displacement amount can be determined. According to this displacement amount, the coordinates of the node can be updated. Thus, for the bonded joint element in the three-dimensional model area, it can be determined whether the bonded joint element has been damaged according to the updated coordinates of the nodes of the bonded joint element. If it has been damaged, the type identifier of the bonded joint element can be updated to the type identifier of the fractured joint element. Then, for the fractured joint element, it can be judged whether the fractured joint element is a hydraulic joint element. If so, 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 seepage calculation are performed on the hydraulic joint element.

[0069] For example, the class identifier of the bonded 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 bonded joint element is damaged, that is, it becomes a fractured joint element, then the type identifier of the joint element is updated from 0 to 1. If this fractured joint element is a hydraulic joint element, then the type identifier of the joint element can be updated from 1 to 2.

[0070] Repeat the above calculation steps of the stress information of the joint element, the calculation steps of the total equivalent nodal force of the node, and the update steps of the joint element type identifier for each time step until the calculation time step reaches the target time step and ends, that is, the three-dimensional numerical simulation of the seepage softening of the rock mass discontinuity surface is completed.

[0071] In the embodiments of the present application, a time-dependent seepage-softening model that takes into account the geometric morphology characteristics and strength characteristics of discontinuity surfaces can ingeniously and intuitively reflect the full-process shear mechanical characteristics of rock mass discontinuity surfaces during tangential sliding. It can improve the accuracy of the total equivalent nodal forces acting on the nodes of joint elements, improve the prediction accuracy of rock mass disasters, and then effective prevention and control strategies can be adopted. It is applicable to the fluid-solid coupling calculation of discontinuity surfaces of various rocks, improving the robustness and applicability. In addition, when calculating the softening process of hydraulic joint elements, the attenuation of the strength of the discontinuity surface wall with immersion time and the attenuation of the friction angle with immersion time are considered. Based on the interface softening function and the friction angle attenuation function, the modified peak shear strength and the modified residual friction angle are obtained. Based on the modified peak shear strength and the modified residual friction angle, the tangential frictional force of the hydraulic joint element is calculated, further improving the accuracy of the tangential frictional force and further improving the accuracy of the equivalent nodal force converted from the tangential frictional force.

[0072] Figure 7 It is a block diagram of a three-dimensional numerical simulation system for seepage softening of rock mass discontinuity surfaces provided by an embodiment of the present application. As Figure 7 shown, the three-dimensional numerical simulation system for seepage softening of rock mass discontinuity surfaces may include: a preprocessing module 701, a node topology structure updating module 702, a calculation module 703, and a type updating module 704.

[0073] Among them, the preprocessing module 701 is used to generate a surface region used to represent various discontinuity surfaces in the rock mass as a mesh surface within the three-dimensional rock model region, and divide the mesh surface into a tetrahedral solid element mesh to obtain the tetrahedral solid element of the three-dimensional rock model region.

[0074] The node topology structure updating module 702 is used to update the node topology structure of the tetrahedral solid element to obtain multiple joint elements in the three-dimensional rock model region.

[0075] The calculation module 703 is used to, at the initial time step, determine the stress information of each joint element in multiple joint elements by using a stress algorithm matching 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 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.

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

[0077] A type update module 704 is configured to update the type identifier of each joint element among a plurality of joint elements based on the total equivalent nodal force at each node in the plurality of joint elements, and perform the total equivalent nodal force of each node in the next time step to update the type identifier of each joint element.

[0078] In some embodiments, the contact force includes a tangential frictional force. The calculation module 703 is configured to: obtain an interface softening coefficient and a friction angle decay function corresponding to the rock; wherein, the interface softening coefficient is a decay function of a parameter measuring the strength of the discontinuous surface wall with respect to the immersion time, and the friction angle decay function is a decay function of the friction angle with respect to 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 a corrected peak shear strength according to the corrected basic friction angle and the corrected rock tensile strength parameter; and determine the tangential frictional force according to the corrected peak shear strength.

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

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

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

[0082] In some embodiments, the computing module 703 is further configured to: when the second fluid pressure is less than the first threshold, assign the second fluid pressure to the first threshold, 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 rate of the hydraulic joint unit with any node cavity as the starting node of the seepage direction. Exemplarily, the computing module 703 is configured to: when the first node saturation is less than the second threshold, determine the second node saturation according to the first node saturation and the sum of the flow rates.

[0083] In some embodiments, the computing module 703 is further configured to: when the second node saturation is greater than the second threshold, assign the second node saturation to the second threshold, 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 sum of the flow rates of all the hydraulic joint units to which any node cavity belongs.

[0084] It should be noted that the foregoing explanation of the embodiments of the three-dimensional numerical simulation method for seepage softening of rock mass discontinuities also applies to the three-dimensional numerical simulation system for seepage softening of rock mass discontinuities in this embodiment, and will not be repeated here.

[0085] To implement the above embodiments, the present application also provides an electronic device, including: at least one processor; and a memory communicatively connected to 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 so that the at least one processor can execute the method for implementing the foregoing embodiments.

[0086] To implement the above embodiments, the present application also provides a storage medium storing instructions that, when run on an electronic device, cause the electronic device to execute the method for implementing the foregoing embodiments.

[0087] To implement the above embodiments, the present application also provides a program product including at least one of a program and instructions, and when at least one of the program and instructions is executed by a processor, the steps of the method for implementing the foregoing embodiments are realized.

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

[0089] In the description of the foregoing embodiments, the descriptions with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples", etc. mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, without contradiction, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0090] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include at least one of such features. In the description of the present application, "a plurality of" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0091] Any process or method description in a flowchart or described otherwise herein can be understood to represent a module, segment, or portion of code including one or more executable instructions for implementing a customized logic function or process, and the scope of the preferred embodiments of the present application includes additional implementations, where the functions may be executed in a substantially simultaneous manner or in a reverse order according to the functions involved, rather than in the order shown or discussed, which should be understood by those skilled in the art to which the embodiments of the present application pertain.

[0092] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as a definitional sequence list of executable instructions for implementing logical functions, which can be specifically implemented in any computer-readable medium for use by an instruction execution system, apparatus, or device (such as a computer-based system, a system including a processor, or other systems that can fetch and execute instructions from the instruction execution system, apparatus, or device), or in conjunction with these instruction execution systems, apparatuses, or devices. For the purposes of this specification, a "computer-readable medium" can be any device that can contain, store, communicate, propagate, or transport a program for use by or in conjunction with an instruction execution system, apparatus, or device. More specific examples (non-exhaustive list) of computer-readable media include the following: electrical connection parts with one or more wirings (electronic devices), portable computer disk cartridges (magnetic devices), random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber devices, and portable compact disc read-only memory (CDROM). Additionally, the computer-readable medium can even be paper or other suitable media on which the program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other media, followed by editing, interpretation, or other suitable processing as necessary, and then storing it in a computer memory.

[0093] It should be understood that various parts of the present application can be implemented using hardware, software, firmware, or a combination thereof. In the above-described embodiments, multiple steps or methods can be implemented using software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented using hardware, as in another embodiment, it can be implemented using a combination of any one or more of the following techniques known in the art: discrete logic circuits having logic gate circuits for implementing logical functions on data signals, application-specific integrated circuits having suitable combinational logic gate circuits, programmable gate arrays (PGA), field-programmable gate arrays (FPGA), etc.

[0094] Those of ordinary skill in the art of this technology can understand that all or part of the steps carried by the method of implementing the above embodiments can be completed by instructing relevant hardware through a program, and the program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiments.

[0095] In addition, each functional unit in various embodiments of the present application may be integrated into one processing module, or each unit may exist physically alone, or two or more units may be integrated into one module. The above-mentioned integrated module may 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 may also be stored in a computer-readable storage medium.

[0096] The above-mentioned storage medium may be a read-only memory, a magnetic disk or an optical disc, etc. Although the embodiments of the present application have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present application. Those of ordinary skill in the art can make changes, modifications, substitutions, 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, Including: Generating a region used to represent various discontinuity surfaces in the rock mass within the three-dimensional rock model region as a mesh surface, and dividing the mesh surface into tetrahedral solid element meshes to obtain the tetrahedral solid elements of the three-dimensional rock model region; Updating the node topological structure of the tetrahedral solid elements to obtain multiple joint elements in the three-dimensional rock model region; At the initial time step, according to the type identifier of each joint element among the multiple joint elements, using a stress algorithm matching the type identifier 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 on both walls of the hydraulic joint element calculated based on the softening process and the first fluid pressure at each node of the hydraulic joint element calculated based on seepage; Based on the stress information of each joint element, determining the total equivalent nodal force at each node among the multiple joint elements; Based on the total equivalent nodal force at each node among the multiple joint elements, updating the type identifier of each joint element among the multiple joint elements, and performing the total equivalent nodal force at each node in the next time step to update the type identifier of each joint element.

2. The method according to claim 1, wherein The contact force includes tangential frictional force; performing a softening process calculation on the hydraulic joint element to determine the contact force acting on both walls of the hydraulic joint element, including: Obtaining an interface softening coefficient and a friction angle attenuation function corresponding to the rock; wherein, the interface softening coefficient is a decay function of a parameter measuring the strength of the discontinuity surface wall with respect to the immersion time, and the friction angle attenuation function is a decay function of the friction angle with respect to the immersion time; According to the interface softening coefficient, correcting the basic friction angle to obtain a corrected basic friction angle; According to the friction angle attenuation function, correcting the rock tensile strength parameter to obtain a corrected rock tensile strength parameter; According to the corrected basic friction angle and the corrected rock tensile strength parameter, determining a corrected peak shear strength; According to the corrected peak shear strength, determining the tangential frictional force.

3. The method according to claim 2, wherein The determining the tangential frictional force according to the corrected peak shear strength includes: Determining the change amount of the sliding distance on the interaction edge of the hydraulic joint element and the sliding distance corresponding to the peak shear strength; In the case where the absolute value of the change amount of the sliding distance is less than or equal to the absolute value of the sliding distance, determining the tangential frictional force according to the corrected peak shear strength and the change amount of the sliding distance; In the case where the absolute value of the sliding distance variable is greater than the sliding distance, correcting the residual friction angle according to the friction angle attenuation function, and determining a corrected residual shear strength according to the corrected residual friction angle; Determining the tangential frictional force according to the corrected peak shear strength, the change amount of the sliding distance, and the corrected residual shear strength.

4. The method according to claim 1, characterized in that Performing a seepage calculation on the hydraulic joint element to determine the first fluid pressure at each node of the hydraulic joint element, including: For any node cavity among the nodes of the hydraulic joint element, determine the first fluid pressure of the any node cavity at the current time step according to the second fluid pressure of the any node cavity at the previous time step.

5. The method according to claim 4, wherein The determining the first fluid pressure of the any node cavity at the current time step according to the second fluid pressure of the any node cavity at the previous time step includes: When the second fluid pressure is greater than the first threshold, determine the first fluid pressure according to the second fluid pressure and the sum of the flows of all the hydraulic joint elements to which the any node cavity belongs.

6. The method according to claim 5, wherein The method further includes: When the second fluid pressure is less than the first threshold, assign the second fluid pressure to the first threshold, and determine the second node saturation of the any node cavity at the current time step according to the first node saturation of the any node cavity at the previous time step; wherein, the second node saturation is used to determine the flow of the hydraulic joint element with the any node cavity as the starting node of the seepage direction.

7. The method according to claim 6, wherein The determining the second node saturation of the any node cavity at the current time step according to the first node saturation of the any node cavity at the previous time step includes: When the first node saturation is less than the second threshold, determine the second node saturation according to the first node saturation and the sum of the flows.

8. The method according to claim 6 or 7, characterized in that The method further includes: When the second node saturation is greater than the second threshold, assign the second node saturation to the second threshold, and determine the fluid pressure of the any node cavity at the next time step according to the fluid pressure of the any node cavity at the current time step and the sum of the flows of all the hydraulic joint elements to which the any node cavity belongs.

9. A three-dimensional numerical simulation system for seepage softening of rock mass discontinuities, characterized in that, Includes: A preprocessing module, configured to generate a region representing various discontinuity surfaces in the rock mass as a mesh surface within the three-dimensional rock model region, and divide the mesh surface into a tetrahedral solid element mesh to obtain the tetrahedral solid element of the three-dimensional rock model region; A node topology structure updating module, configured to update the node topology structure of the tetrahedral solid element to obtain a plurality of joint elements in the three-dimensional rock model region; A calculation module, configured to, at the initial time step, according to the type identifier of each joint element among the plurality of joint elements, use a stress algorithm matching the type identifier 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 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; The calculation module is further configured to determine the total equivalent nodal force on each node among the plurality of joint elements based on the stress information of each joint element. A type update module, configured to update the type identifier of each of the plurality of joint units based on the total equivalent nodal force at each node in the plurality of joint units, and perform the total equivalent nodal force at each node in the next time step to update the type identifier of each of the joint units.

10. An electronic device, characterized in that, Comprising: At least one processor; And A memory communicatively connected to the at least one processor; wherein, The memory stores instructions executable by the at least one processor, and when the instructions are executed by the at least one processor, the at least one processor is enabled to execute the method according to any one of claims 1-8.

Images