Fracture force-seepage-chemical coupling simulation method, medium and equipment

Abstract

The invention provides a fracture force-seepage-chemical coupling simulation method, and relates to the technical field of earth science, and the method comprises the following steps: establishing a rough fracture force-seepage-chemical coupling simulation method for fracture evolution according to a force field model, a seepage field model and a chemical field model, initializing three-field model calculation parameters, and obtaining opening according to coordinates of upper and lower joints; unit contact judgment is conducted according to the opening degree, seepage field and chemical field analysis is conducted on non-contact units, seepage pressure and H + ion concentration of all the units are obtained, the opening degree of the non-contact units is updated, unit contact is judged, and the seepage pressure of the non-contact units is considered; obtaining the normal force, the shearing force, the normal displacement and the shearing displacement of the contact unit by adopting a force field model, and updating the opening degree based on the normal displacement and the shearing displacement; and according to the updated opening degree, new contact and non-contact units are divided until the three-field unknown quantity converges and is in an end time step. The method can describe fracture evolution behaviors under actual working conditions.

Description

Technical Field

[0001] The present invention relates to the field of earth science technology, and in particular to a fracture force-seepage-chemistry coupling simulation method, medium, and equipment. Background Art

[0002] In underground rock engineering and energy development, rock fracture systems largely control the mechanical response, seepage behavior, and material migration pathways of subsurface media. In particular, in projects involving deep geological systems, such as oil and gas extraction, geothermal development, and nuclear waste disposal, variations in permeability within fractures directly impact resource extraction efficiency and engineering stability. Fracture systems, however, are not ideally smooth surfaces but rather possess complex mesostructures, including roughness distribution, irregular contact surfaces, and discontinuous flow channel structures. These mesostructures not only influence seepage pathways and velocity distribution but also induce nonlinear deformation behaviors such as fracture closure and shear slip under external stress, significantly impacting the overall seepage-mechanical coupling relationship.

[0003] At the same time, chemical reactions within the fracture system are also crucial factors that cannot be ignored. This is especially true when injected fluids or subsurface media undergo dissolution and precipitation reactions with the fracture wall rock, further altering the fracture's channel structure and porosity, forming a dynamic feedback mechanism. This complex multi-field interaction of stress, seepage, and chemistry makes it difficult for traditional simplified models to accurately describe fracture evolution under actual operating conditions. Therefore, there is an urgent need to develop a multi-field coupled simulation method that comprehensively considers the influence of fracture mesostructure and integrates mechanical, seepage, and chemical processes. Summary of the Invention

[0004] The purpose of the present invention is to solve the problem that traditional simplified models cannot accurately describe the evolution behavior of fractures under actual working conditions. A fracture force-seepage-chemistry coupling simulation method is proposed, which includes: The rough fracture force field model is constructed based on Hertz contact theory; the rough fracture seepage field model is constructed using the cubic law and mass conservation equation; + The migration process of ion concentration in fluid is used to construct a rough fracture chemical field model; Based on the three-field model, a rough fracture force-seepage-chemistry coupling simulation method is established to simulate fracture evolution. The steps are as follows: Step 1: Initialize the calculation parameters of the rough fracture force field model, seepage field model, and chemical field model, as well as the external force conditions, seepage pressure, and injection fluid concentration, and input the coordinates of the upper and lower joints to obtain the aperture based on the coordinates of the upper and lower joints; Step 2: The contact of the fracture units is judged according to the opening, and the rough fracture seepage field model is used to analyze the seepage field of the uncontacted units to obtain the seepage pressure of each unit; the rough fracture chemical field model is used to analyze the concentration field of H+ ion migration to obtain the H + ion concentration and update the opening of the uncontacted unit; Step 3: Update the uncontacted unit and the contacted unit according to the updated uncontacted unit opening. Consider the seepage pressure of the uncontacted unit and use the rough fracture force field model to obtain the normal force, shear force, normal displacement, and shear displacement of the contacted unit. Update the opening based on the normal displacement and shear displacement. Step 4: Divide the new contact and non-contact units according to the updated opening, and determine whether the three-field unknown quantities have converged. If not, return to step S2; if converged, determine whether it is at the end time step. If not, return to step S2. If it is at the end time step, end the calculation.

[0005] Furthermore, the rough fracture force field includes normal force and shear force. The contact element describes the contact force by two-ball contact. The rough fracture force field model includes the following formula: , , in, represents the normal force of the contact element, represents the shear force of the contact element, R represents the equivalent radius of the sphere, E represents the elastic modulus of the sphere, v represents Poisson's ratio, represents the normal deformation, and θ represents the angle between the two spheres and the horizontal.

[0006] Furthermore, the rough fracture seepage field model includes the following formula: The seepage control equation is: , , Where x and y represent the horizontal and vertical coordinates, respectively, K represents the permeability, b represents the opening of random rough fractures, and μ is the dynamic viscosity of the fluid; The finite difference method is used to discretize the seepage control equation and obtain: , Where h represents the water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Indicates the spacing of grid nodes in the horizontal direction, Indicates the spacing of grid nodes in the vertical direction; Mesh Node The conductivity in the four adjacent directions is defined as: , in, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node penetration rate; node Water head The discrete form of is: .

[0007] Furthermore, the rough fracture chemical field model includes the following formula: H + The transport equation of ion concentration in fluid is: , Where C represents H + Ion concentration; u represents velocity vector; t represents time; In , ▽ represents the convection derivative operator; represents the divergence; D represents the concentration diffusion coefficient; represents the gradient of C; The finite difference method was used to + The ion concentration is discretely expressed in the fluid transport equation, and we get: , in, Indicates that the n+1th time layer is at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at 、 Respectively represent the grid in x Direction and y The grid spacing in the direction, represents the time step, 、 Represents the fluid in x and y The velocity component in the direction; The change in crack opening is expressed as: , in, is the change in crack opening, is the initial moment H + Ion concentration, A represents concentration B represents the exponential sensitivity of the crack opening response to concentration changes.

[0008] Furthermore, the upper joint coordinates minus the lower joint coordinates can be used to obtain the opening, which is expressed as: , in, Indicates the opening degree, represents the upper joint coordinates, Indicates the lower joint coordinates.

[0009] Furthermore, an opening greater than 0 is a non-contact unit; an opening less than or equal to 0 is a contact unit.

[0010] Furthermore, considering the seepage pressure of the non-contact unit, the total normal force of the contact unit is the sum of the normal forces of all contact points and the seepage pressure of the non-contact unit; the total shear force of the contact unit is the sum of the shear forces of all contact points.

[0011] Furthermore, the following formula is used to determine whether the three unknowns converge: , , , in, represents the fluid seepage field pressure distribution obtained in the K-1th coupling iteration, represents the fluid seepage field pressure distribution obtained in the K-th coupling iteration, represents the chemical field concentration distribution obtained in the K-1th coupling iteration, represents the chemical field concentration distribution obtained in the K-th coupling iteration, represents the normal stress on the contact surface obtained in the K-1th coupling iteration, represents the normal stress on the contact surface obtained in the K-th coupling iteration, represents the contact surface shear force obtained in the K-1th coupling iteration, represents the shear force on the contact surface obtained in the Kth coupling iteration.

[0012] The present invention also provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned fracture force-permeation-chemistry coupling simulation method is implemented.

[0013] The present invention also proposes an electronic device, comprising a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program includes computer-readable instructions, and the processor is configured to call the computer-readable instructions to execute the above-mentioned fracture force-seepage-chemistry coupling simulation method.

[0014] The beneficial effects brought about by the technical solution provided by the present invention are: This method achieves dynamic simulation and accurate prediction of fracture aperture evolution by constructing a coupled mechanics-seepage-chemistry model. By discretely processing joint surfaces and updating contact states in real time, combined with a three-field convergence mechanism to control iterative solutions and aperture-concentration responses, it effectively reflects the true evolution of fractures under chemical dissolution and mechanical action. This method can be used to assess rock stability in complex environments such as chemical erosion and pressure seepage, providing a reliable basis for safe design and long-term service prediction of underground projects. It possesses strong engineering applicability and theoretical value. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 This is a flow chart of a fracture force-seepage-chemistry coupling simulation method according to an embodiment of the present invention; Figure 2 Schematic diagram of crack contact mechanics analysis according to an embodiment of the present invention, wherein: Figure 2 (a) is the mechanical analysis of the contact surface. Figure 2 (b) is the mechanical analysis of the contact element; Figure 3 This is a schematic diagram of seepage flow according to an embodiment of the present invention; Figure 4 This is the entire rough fracture force-seepage-chemistry coupling calculation process established in the embodiment of the present invention; Figure 5 It is a block diagram of an electronic device in an exemplary embodiment of the present invention. DETAILED DESCRIPTION

[0016] To make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0017] The flow chart of the fracture force-seepage-chemistry coupling simulation method according to the embodiment of the present invention is as follows: Figure 1 , specifically including: (1) The rough crack force field model is constructed based on the Hertz contact theory. The force field solution is divided into two parts: normal force and shear force. The calculation theory of both is based on the Hertz contact theory, which uses the contact of two points to simplify the contact of two balls (contact units) to describe the contact force, such as Figure 2 As shown. First, the normal force can be calculated using the following expression, , in, represents the normal force of the contact element, R represents the equivalent radius of the sphere, ; is the radius of the upper joint, is the radius of the lower joint; E is the elastic modulus of the sphere, v represents Poisson's ratio, Indicates the normal deformation.

[0018] Consider the equilibrium of forces on a single sphere, such as Figure 2 As shown in (b), then: , Here, θ is the angle between the two spheres and the horizontal.

[0019] Thus we can get: , in, Represents the shear force of the contact element.

[0020] At the same time, assuming that the shear deformation and shear force satisfy the linear relationship, then: , in, represents the deformation stiffness of the two balls, Represents the tangential relative displacement.

[0021] For all contact points i, the following static equilibrium relationship exists: , , in, represents the normal force at contact point i on the contact element, represents the shear force at contact point i on the contact element, represents the normal stress acting on the contact surface, represents the tangential stress acting on the contact surface.

[0022] According to the above formula, the normal force on each contact body can be calculated and shear force .

[0023] (2) The cubic law and mass conservation equation are used to construct a rough fracture seepage field model.

[0024] The non-contact zone is where the fluid flows. Its flow behavior is described using the cubic law and the mass conservation equation. Darcy's law reveals the law of fluid transfer in the medium through the relationship between permeability and pressure gradient. The mass conservation equation ensures the conservation of mass throughout the system. After discretizing the non-contact zone, the seepage control equation for the non-contact unit is obtained as follows: , , Where x and y represent the horizontal and vertical coordinates respectively; h is the water head in m; K represents the permeability in m 2 ; b represents the opening of random rough cracks, unit is m; μ is the dynamic viscosity of the fluid, unit is Pa·s.

[0025] In order to solve numerically, the above equations are discretized using the finite difference method. The seepage solution area is divided into grids with a grid node spacing of and For nodes , the second-order partial derivative of the hydraulic head h is expressed as the following central difference: , in, h Indicates water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Indicates the spacing of grid nodes in the horizontal direction, Indicates the spacing of grid nodes in the vertical direction.

[0026] Since the cracks have different openings in each unit after discretization, and the crack opening has a significant impact on the permeability of the cracks, the conductivity factor T is proposed to comprehensively consider the differences in the actual conductivity calculation caused by the different openings of each crack unit. , and the grid nodes The nodes in the four adjacent directions are: The node in the east is: ; The nodes in the west are: ; The nodes in the north are: ; The nodes in the south are: . Grid Node The conductivity of the nodes in the four adjacent directions is defined as: , in, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node penetration rate; node Water head The discrete form of is: .

[0027] For the seepage unit in the uncontacted area, the above solution method has the following characteristics: ① Compared with traditional finite element methods, this solution method converges significantly faster. When coupling reaches a certain level, a large number of scattered contact zones exist within the crack. In the finite element solution, an internal boundary must be established for each internal contact zone to solve the problem. When the number of internal boundaries in the solution domain is large, the finite element solution will converge slowly or even fail to converge, resulting in an inability to solve the problem. This method avoids this problem, ensuring solution speed and accuracy, and fully considering the microscopic contact state within the crack.

[0028] ② The finite difference method has a local structure and clear physical meaning, which makes it easy to embed the dynamic changes of anisotropic permeability and fracture aperture, and can flexibly adapt to the evolution process of heterogeneous rock mass and fracture structure.

[0029] ③ Its calculation format is regular and the coefficient matrix is ​​sparse, which is suitable for parallel acceleration and grid adaptive update. It has good numerical stability and scalability, and can efficiently handle the seepage solution of large-scale complex fracture networks in coupled simulations. Figure 3 shown. Figure 3 The crack shown in is composed of upper and lower joint surfaces. The crack has a seepage inlet and a seepage outlet, and the boundaries on both sides are no-flow boundaries, that is, there is no fluid exchange between the boundary and the outside world. After discretization, the j1 group of sub-cracks is magnified for illustration. There are three types of units in the crack, (a) is a contact unit, no fluid passes through the unit, (b) is a seepage unit, this unit can normally undertake the seepage task, (c) is an adjacent boundary unit, the special feature of this unit is that there is a no-flow boundary with a crack boundary, and there is no liquid exchange in this direction. In the enlarged view of the three units, there are the following symbols: double arrows represent that there is fluid exchange between the unit and the adjacent unit, and × represents no fluid exchange. Grid unit The adjacent grid cells above, below, left and right are 、 、 、 , 、 、 Represents grid cells 、 、 The opening degree.

[0030] (3) According to H + The migration process of ion concentration in fluid is used to construct a rough fracture chemical field model.

[0031] Fluid H + Ions, which produce chemical corrosion on the surface of rock cracks. This chemical corrosion is related to H + The ion concentration is related to H +The migration of ion concentration in the fluid satisfies the following equation: , Where C represents H + Ion concentration, u represents velocity vector, t represents time, In , ▽ represents the convection derivative operator, represents the divergence, D represents the concentration diffusion coefficient, represents the gradient of C.

[0032] The finite difference method was used to + The ion concentration is discretely expressed in the fluid transport equation, and we get: , in, Indicates that the n+1th time layer is at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at 、 Respectively represent the grid in x The grid spacing in the a and y directions, represents the time step, 、 Represents the fluid in x and the velocity component in the y direction.

[0033] The H of each unit can be obtained by the above formula + Ion concentration. + The degree of corrosion of the crack surface due to the ion concentration is different, resulting in changes in the crack opening. The change in crack opening is related to the consumption of H+ ions in the solution. The change in crack opening is expressed as: , in, is the change in crack opening, is the initial moment H + Ion concentration, A represents concentration The change in crack opening under the condition of , B represents the exponential sensitivity of concentration change to crack opening response, and both are determined by least squares fitting of the logarithmic linearization of the experimental data.

[0034] (4) Based on the three-field model, a rough fracture force-seepage-chemistry coupling simulation method is established to simulate fracture evolution, referring to Figure 4 , Figure 4 The entire rough fracture force-permeability-chemistry coupling calculation process established in the embodiment of the present invention has the following steps: Step 1: Initialize the calculation parameters of the rough fracture force field model, seepage field model, and chemical field model, as well as the external force conditions, seepage pressure, and injection fluid concentration. Enter the coordinates of the upper and lower joints, and obtain the aperture by subtracting the lower joint coordinates from the upper joint coordinates: ,in, Indicates the opening degree, represents the upper joint coordinates, Indicates the lower joint coordinates.

[0035] Step 2: Determine contact between fracture elements based on their openings. Elements with openings greater than 0 are considered non-contact elements; elements with openings less than or equal to 0 are considered contact elements. A rough fracture seepage field model is used to analyze the seepage field of non-contact elements to obtain the seepage pressure of each element. Considering the input seepage pressure boundary, the Newton-Raphson iterative method is used to solve the problem until convergence, obtaining the hydraulic head and flow velocity of each element. Based on the relationship between hydraulic head and seepage pressure: , calculate the seepage pressure of each unit , is the unit head, ρ is the fluid density, and g is the acceleration due to gravity.

[0036] Based on the velocity field of each unit, the rough fracture chemical field model is used to analyze the H + The concentration field of ion migration is analyzed, and the concentration discrete format of each unit is integrated using the finite difference method. Considering the input concentration boundary, the Newton-Raphson iterative method is used to solve until convergence, and the H of each untouched unit is obtained. + Ion concentration, according to the formula Update the opening of each untouched element.

[0037] Step 3: Update the uncontacted unit and the contacted unit according to the updated uncontacted unit opening. Consider the seepage pressure of the uncontacted area in the force balance equation according to the new contacted unit and uncontacted unit. , assuming that all contact normal displacements and shear displacements are equal and that the non-contact area does not bear shear force, the new force balance equation is: , , The new mechanical equilibrium equations are solved by Gauss-Seidel iteration, and the normal force, shear force, normal displacement and shear displacement of the contact element are obtained based on the rough crack force field model. The opening is updated based on the normal displacement and shear displacement. The method for updating the opening is to first subtract the normal displacement from the opening of the non-contact element. , get the new opening, record the new coordinates of the upper and lower joints, and then use shear displacement , determine the upper joint movement distance, taking the x direction as an example, and calculate it according to the following formula: , Among them, Δx is the length of the unit in the x direction, int() is the rounding function, represents the shear displacement, and substituting the coordinates of the newly obtained upper and lower joints into the opening value is: , Where N represents the number of translation grid points.

[0038] Step 4: Divide the new contact and non-contact units according to the updated opening, and determine whether the unknown quantities of the force field, seepage field, and chemical field have converged. If not, return to step S2; if converged, determine whether it is at the end time step. If not, return to step S2. If it is at the end time step, end the calculation.

[0039] Use the following formula to determine whether the three unknowns converge: , , , in, represents the fluid seepage field pressure distribution obtained in the K-1th coupling iteration, represents the fluid seepage field pressure distribution obtained in the K-th coupling iteration, represents the chemical field concentration distribution obtained in the K-1th coupling iteration, represents the chemical field concentration distribution obtained in the K-th coupling iteration, represents the normal stress on the contact surface obtained in the K-1th coupling iteration, represents the normal stress on the contact surface obtained in the K-th coupling iteration, represents the contact surface shear force obtained in the K-1th coupling iteration, represents the shear force on the contact surface obtained in the Kth coupling iteration.

[0040] In an exemplary embodiment, a computer-readable storage medium is included, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned fracture force-permeability-chemistry coupling simulation method is implemented.

[0041] See also Figure 5 In an exemplary embodiment, an electronic device is also included, including at least one processor, at least one memory, and at least one communication bus.

[0042] Wherein, a computer program is stored in the memory, and the computer program includes computer-readable instructions. The processor calls the computer-readable instructions stored in the memory through a communication bus to execute the above-mentioned fracture force-seepage-chemistry coupling simulation method.

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

Claims

1. A fracture force-seepage-chemistry coupling simulation method, characterized in that: include: A rough crack force field model is constructed based on Hertz contact theory; The cubic law and mass conservation equation are used to construct the rough fracture seepage field model; according to H + The migration process of ion concentration in fluid is used to construct a rough fracture chemical field model; Based on the three-field model, a rough fracture force-seepage-chemistry coupling simulation method is established to simulate fracture evolution. The steps are as follows: Step 1: Initialize the calculation parameters of the rough fracture force field model, seepage field model, and chemical field model, as well as the external force conditions, seepage pressure, and injection fluid concentration, and input the coordinates of the upper and lower joints to obtain the aperture based on the coordinates of the upper and lower joints; Step 2: According to the opening, the contact of the fracture units is judged, and the rough fracture seepage field model is used to analyze the seepage field of the non-contact units to obtain the seepage pressure of each unit; the rough fracture chemical field model is used to analyze the H + The concentration field of ion migration is analyzed to obtain the H of each unit. + ion concentration and update the opening of the uncontacted unit; Step 3: Update the uncontacted unit and the contacted unit according to the updated uncontacted unit opening. Consider the seepage pressure of the uncontacted unit and use the rough fracture force field model to obtain the normal force, shear force, normal displacement, and shear displacement of the contacted unit. Update the opening based on the normal displacement and shear displacement. Step 4: Divide the new contact and non-contact units according to the updated opening, and determine whether the three-field unknown quantities have converged. If not, return to step S2; if converged, determine whether it is at the end time step. If not, return to step S2. If it is at the end time step, end the calculation.

2. A fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: The rough fracture force field includes normal force and shear force. The contact element describes the contact force using two-ball contact. The rough fracture force field model includes the following formulas: ， ， in, represents the normal force of the contact element, represents the shear force of the contact element, R represents the equivalent radius of the sphere, E represents the elastic modulus of the sphere, v represents Poisson's ratio, represents the normal deformation, θ Represents the angle between the two spheres and the horizontal.

3. A fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: The rough fracture seepage field model includes the following formulas: The seepage control equation is: ， ， Where x and y represent the horizontal and vertical coordinates, respectively, K represents the permeability, b represents the opening of random rough fractures, and μ is the dynamic viscosity of the fluid; The finite difference method is used to discretize the seepage control equation and obtain: ， Where h represents the water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Representation node The water head, Indicates the spacing of grid nodes in the horizontal direction, Indicates the spacing of grid nodes in the vertical direction; Mesh Node The conductivity in the four adjacent directions is defined as: ， in, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The conductivity, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node The penetration rate, Representation node penetration rate; node Water head The discrete form of is: 。 4. A fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: The rough fracture chemical field model includes the following formulas: H + The transport equation of ion concentration in fluid is: ， Where C represents H + Ion concentration; u represents the velocity vector; t represents time; In , ▽ represents the convection derivative operator; represents the divergence; D represents the concentration diffusion coefficient; represents the gradient of C; The finite difference method was used to + The ion concentration is discretely expressed in the fluid transport equation, and we get: ， in, Indicates that the n+1th time layer is at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at Indicates the nth time layer at the grid point The concentration value at 、 Respectively represent the grid in x Direction and y The grid spacing in the direction, represents the time step, 、 Represents the fluid in x and y The velocity component in the direction; The change in crack opening is expressed as: ， in, is the change in crack opening, is the initial moment H + Ion concentration, A represents concentration B represents the exponential sensitivity of the crack opening response to concentration changes.

5. The fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: Subtract the lower joint coordinates from the upper joint coordinates to get the opening, which is expressed as: ， in, Indicates the opening degree, represents the upper joint coordinates, Indicates the lower joint coordinates.

6. A fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: If the opening is greater than 0, it is a non-contact element; if the opening is less than or equal to 0, it is a contact element.

7. The fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: Considering the seepage pressure of the non-contact unit, the total normal force of the contact unit is the sum of the normal forces of all contact points and the seepage pressure of the non-contact unit; the total shear force of the contact unit is the sum of the shear forces of all contact points.

8. The fracture force-seepage-chemistry coupling simulation method according to claim 1, characterized in that: Use the following formula to determine whether the three unknowns converge: ， ， ， in, represents the fluid seepage field pressure distribution obtained in the K-1th coupling iteration, represents the fluid seepage field pressure distribution obtained in the K-th coupling iteration, represents the chemical field concentration distribution obtained in the K-1th coupling iteration, represents the chemical field concentration distribution obtained in the K-th coupling iteration, represents the normal stress on the contact surface obtained in the K-1th coupling iteration, represents the normal stress on the contact surface obtained in the K-th coupling iteration, represents the contact surface shear force obtained in the K-1th coupling iteration, represents the shear force on the contact surface obtained in the Kth coupling iteration.

9. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 8 is implemented.

10. An electronic device, characterized in that: The method comprises a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program includes computer-readable instructions, and the processor is configured to call the computer-readable instructions to execute the method according to any one of claims 1 to 8.

Images