A three-dimensional rough fracture nonlinear flow field simulation method and system

Abstract

The present application provides a method and system for simulating a nonlinear flow field of a three-dimensional rough fracture, which belongs to the field of electrical digital data processing technology. The method includes: generating an initial geometric opening field of a three-dimensional rough fracture based on the three-dimensional coordinates of the upper wall of the three-dimensional rough fracture obtained by three-dimensional laser scanning technology; using the Reynolds equation with a time-varying term of the opening to describe the head distribution of the linear flow field in the fracture, embedding the Forchheimer equation into the Reynolds equation with a time-varying term of the opening, obtaining an improved Reynolds equation for describing the nonlinear flow field, and solving the head distribution of the nonlinear flow field under the control of the improved Reynolds equation based on the preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture. By embedding the Forchheimer equation into the Reynolds equation, the description capability of the nonlinear flow field is improved, and the fluid behavior in the three-dimensional rough fracture can be described and predicted more accurately, which not only broadens its scope of application, but also lays the foundation for conducting research on the long-term dissolution mechanism of three-dimensional rough fractures.

Description

Technical Field

[0001] The present application relates to the technical field of electrical digital data processing, and in particular to a three-dimensional rough fracture nonlinear flow field simulation method and system. Background Art

[0002] The cracks in the dam foundation rock mass may be both potential sliding surfaces and the main interface of water-rock interaction. Usually, low-solubility minerals in the cracks in the dam foundation rock mass will dissolve under the action of the seepage field and move to the crack outlet through convection and diffusion. Under short-term seepage conditions, due to the extremely low solubility of minerals, the seepage volume of the dam foundation is small, the amount of dissolution of the crack wall caused by water-rock interaction is small, and the change of the crack wall structure is not obvious; but under long-term seepage conditions, the rough wall surface undergoes qualitative changes due to dissolution accumulation, and the solid wall morphology structure, flow field material concentration, pore wall water pressure, etc. undergo significant changes, causing macroscopic dissolution of the bedrock cracks, further aggravating the leakage of the dam foundation cracks. For dam foundation cracks with strong solubility, the problem of dam foundation crack dissolution is more prominent. Therefore, studying the dissolution characteristics and dissolution mechanism of dam foundation cracks over a longer time scale is of great significance for the long-term safe and efficient operation of the dam.

[0003] At present, most of the studies on water-rock interaction in the literature focus on the influence of water immersion on the macroscopic mechanical properties of joints, or the influence of seepage on mass transfer in the aperture field formed by fracture shear slip, but rarely consider the structure of bedrock fractures under long-term water-rock interaction and the gradual evolution characteristics of various physical fields over time. Some studies have analyzed the changes in pore or fracture structure caused by the action of reactive liquids on fractures from the perspective of chemical engineering, but the time scale is short (tens of seconds or hours), and the fractures are dissolved in a relatively fast time, which is inconsistent with the leakage characteristics of dam foundation fractures (low reactivity, low material concentration, less or no reactive fluid, and mainly clear water dissolution). Although there are a few studies that analyze the fracture dissolution characteristics over a longer time scale, most of these coupling models use the Reynolds equation with time-varying aperture to calculate the flow field. However, the Reynolds equation assumes that the flow rate at each point in the rough fracture satisfies the cubic law (so the Reynolds equation is also called the local cubic law), and the flow capacity calculated by this equation is linear. A large number of experimental and numerical studies have shown that the pressure gradient and volume flow in rough fractures are nonlinear; therefore, the Reynolds equation with time-varying terms of opening is not sufficient to reflect the true flow capacity of rough fractures. In addition, the long-term dissolution mechanism and dissolution characteristics of rough fractures are still unclear in the literature.

[0004] Therefore, it is necessary to provide an improved technical solution to address the above-mentioned deficiencies in the prior art. Summary of the invention

[0005] The purpose of the present application is to provide a three-dimensional rough fracture nonlinear flow field simulation method and system to solve or alleviate the problems existing in the above-mentioned prior art.

[0006] In order to achieve the above objectives, this application provides the following technical solutions:

[0007] In a first aspect, the present application provides a three-dimensional rough fracture nonlinear flow field simulation method, comprising:

[0008] Based on the three-dimensional coordinates of the upper wall of the three-dimensional rough crack obtained by three-dimensional laser scanning technology, the initial geometric opening field of the three-dimensional rough crack is generated;

[0009] The Reynolds equation with a time-varying term of the opening is used to describe the head distribution of the linear flow field in the fracture. The Forchheimer equation is embedded in the Reynolds equation with a time-varying term of the opening to obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows:

[0010] ,

[0011] In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = T s ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )for t time(x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 ;

[0012] Based on the preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture, the water head distribution of the nonlinear flow field under the control of the improved Reynolds equation is solved.

[0013] In the embodiment of the present application, a three-dimensional rough surface is obtained by a three-dimensional laser scanning method to generate a rough fracture space, and the nonlinear Forchheimer equation is embedded in the Reynolds equation with a time-varying term of the opening to achieve the simulation of nonlinear flow in rough fractures under transient conditions. By embedding the Forchheimer equation into the Reynolds equation, the ability to describe the nonlinear flow field is improved, and common nonlinear flow phenomena in rough fractures, such as turbulence, local eddies, etc., can be simulated. The enhanced nonlinear processing capability enables the model to be applicable to a wider range of flow conditions. By combining high-precision geometric data acquisition, dynamic nonlinear flow models, and flexible boundary condition settings, it is possible to more accurately describe and predict the fluid behavior in three-dimensional rough fractures, broadening its scope of application. It is an ideal tool for studying the flow process in complex fracture media, and also lays the foundation for conducting research on the long-term dissolution mechanism of three-dimensional rough fractures.

[0014] In combination with the first aspect, in some optional implementations, the Forchheimer equation is embedded into the Reynolds equation with the time-varying term of the opening, including:

[0015] Use fracture geometry b Alternative fracture hydraulic aperture e h , used to solve the linear coefficients in the Forchheimer equation A and nonlinear coefficientB , and then obtain the linear coefficient after replacement A and nonlinear coefficient B Expression of

[0016] The linear coefficients after replacement A and nonlinear coefficient B Substituting into the Forchheimer equation, we get the fluid volume flow rate Q The first expression of

[0017] Combined with existing research, the fluid volume flow rate is obtained Q The second expression of

[0018] Using the fluid volume flow Q The first expression and the fluid volume flow rate Q The second expression of t time( x , y ) T s ;

[0019] According to the local cubic law, the flow velocity in the opening direction is constant, and the flow velocity in the flow field is determined U =( u , v )

[0020] In the crack t time( x , y ) T s Substitute the velocity into the flow field U =( u , v ) to obtain the improved Reynolds equation.

[0021] In combination with the first aspect, in some optional implementations, the initial geometric opening field of the three-dimensional rough crack is generated based on the three-dimensional coordinates of the upper side wall of the three-dimensional rough crack obtained by the three-dimensional laser scanning technology, including:

[0022] Importing the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture into MATLAB as data points to construct an initial data matrix;

[0023] Generating a three-dimensional simulation grid according to the size of the three-dimensional rough crack, and performing interpolation processing on the initial data matrix to generate a geometric opening simulation surface of the upper side wall of the three-dimensional rough crack;

[0024] The geometric opening simulation surface of the upper side wall of the three-dimensional rough fracture is moved horizontally and vertically to simulate fracture shear sliding, thereby obtaining the geometric opening simulation surface of the lower side wall of the three-dimensional rough fracture;

[0025] The geometric opening of the upper and lower walls is calculated to simulate the opening difference at each grid in the curved surface, and the initial geometric opening field of the three-dimensional rough crack is obtained.

[0026] In a second aspect, this embodiment provides a method for analyzing the long-term dissolution mechanism of three-dimensional rough fractures, comprising:

[0027] The three-dimensional rough fracture nonlinear flow field simulation method provided in any of the above embodiments is used to calculate the flow field of the three-dimensional rough fracture under steady-state conditions, and the concentration field is calculated based on the depth-averaged solute transport equation to obtain the initial distribution of the flow field and the concentration field;

[0028] Through the one-way cycle coupling of flow field, aperture field and concentration field, a fracture dissolution coupling model is established;

[0029] Analyze the difference Da Number and sum P The evolution characteristics of the aperture field, flow field, concentration field, and roughness of three-dimensional rough fractures under the action of number and long-term dissolution are analyzed to reveal the long-term dissolution characteristics and dissolution mechanism of three-dimensional rough fractures.

[0030] In this embodiment, the linear flow field simulation method is used to calculate the three-dimensional flow field under steady-state conditions and the depth-averaged solute transport equation is used to accurately simulate the complex movement process of the fluid in the rough fracture, and to capture the diffusion and convection behavior of the solute in the flow field, which helps to more accurately predict the development trend of the dissolution process in the fracture, reduce the cumulative error, and lay the foundation for analyzing the long-term evolution mechanism of the fracture. Through the unidirectional cycle coupling flow field, aperture field, and concentration field, the complex physical process can be simplified to a certain extent, the complex feedback relationship in the numerical calculation can be reduced within a certain range, and the calculation stability and efficiency can be improved. On this basis, by analyzing different Da Number and sum P The evolution characteristics of the three-dimensional rough fracture aperture field, flow field, concentration field, and roughness under the action of long-term dissolution can provide a deep understanding of the effects of chemical reaction rate, convective transport, and diffusion transport on the fracture dissolution process, which is conducive to predicting the long-term evolution of fractures under different environmental conditions. By studying the evolution characteristics of the aperture field, flow field, concentration field, and roughness under long-term dissolution, it helps to reveal the complex dynamic mechanism in the fracture dissolution process. In addition, this method takes into account both computational efficiency and accuracy, and is an ideal technical path for studying fracture dissolution mechanisms and long-term evolution.

[0031] In a third aspect, this embodiment provides a three-dimensional rough fracture nonlinear flow field simulation system, comprising:

[0032] A simulation unit configured to generate an initial geometric opening field of the three-dimensional rough fracture based on the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture obtained by the three-dimensional laser scanning technology;

[0033] The improved unit is configured to use the Reynolds equation with a time-varying term of opening to describe the water head distribution in the fracture, and the Forchheimer equation is embedded in the Reynolds equation with a time-varying term of opening to obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows:

[0034] ,

[0035] In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = T s ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )for t time( x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 ;

[0036] The solving unit is configured to solve the head distribution of the nonlinear flow field under the control of the improved Reynolds equation based on preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture.

[0037] In conjunction with the third aspect, in some optional implementation manners, the improvement unit is further configured as follows:

[0038] Use fracture geometry b Alternative fracture hydraulic aperture e h , used to solve the linear coefficients in the Forchheimer equation A and nonlinear coefficient B , and then get the linear coefficient after replacement A and nonlinear coefficient B Expression of

[0039] The linear coefficients after replacement A and nonlinear coefficient B Substituting into the Forchheimer equation, we get the fluid volume flow rate Q The first expression of

[0040] Combined with existing research, fluid volume flow rate Q The second expression of

[0041] Using the fluid volume flow Q The first expression and the fluid volume flow rate Q The second expression of t time( x , y ) T s ;

[0042] According to the local cubic law, the flow velocity in the opening direction is constant, and the flow velocity in the flow field is determined U =( u , v )

[0043] In the crack t time( x , y ) T s Substitute the velocity into the flow field U =( u , v ) is used to obtain the improved Reynolds equation for describing nonlinear flow field.

[0044] In conjunction with the third aspect, in some optional implementation manners, the simulation unit is further configured as follows:

[0045] Importing the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture into MATLAB as data points to construct an initial data matrix;

[0046] Generating a three-dimensional simulation grid according to the size of the three-dimensional rough crack, and performing interpolation processing on the initial data matrix to generate a geometric opening simulation surface of the upper side wall of the three-dimensional rough crack;

[0047] The geometric opening simulation surface of the upper side wall of the three-dimensional rough fracture is moved horizontally and vertically to simulate fracture shear sliding, thereby obtaining the geometric opening simulation surface of the lower side wall of the three-dimensional rough fracture;

[0048] The geometric opening of the upper and lower walls is calculated to simulate the opening difference at each grid in the curved surface, and the initial geometric opening field of the three-dimensional rough crack is obtained.

[0049] In a fourth aspect, an embodiment of the present application provides an electronic device, comprising: a memory, a processor, and a program stored in the memory and executable on the processor, wherein when the processor executes the program, any one of the methods provided in the first or second aspect is implemented.

[0050] In a fifth aspect, an embodiment of the present application provides a computer-readable storage medium storing computer-executable instructions, wherein the computer-executable instructions are used to enable a computer to execute any of the methods provided in the first aspect or the second aspect above.

[0051] In a sixth aspect, an embodiment of the present application provides a computer program product comprising computer executable instructions, wherein the computer executable instructions are used to enable a computer to execute any of the methods provided in the first aspect or the second aspect.

[0052] The technical effects that can be achieved in any of the third to sixth aspects mentioned above can refer to the description of the beneficial effects in the first and second aspects mentioned above, and will not be repeated here. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1A schematic diagram of an electronic device provided in this embodiment.

[0054] Figure 2 A schematic flow chart of a three-dimensional rough fracture nonlinear flow field simulation method provided in this embodiment.

[0055] Figure 3 Schematic diagram for modeling a rough fracture in 3D.

[0056] Figure 4 Schematic diagram of the head distribution characteristics of the flow field under the control of the Reynolds equation and the improved Reynolds equation.

[0057] Figure 5 Schematic diagram of the head distribution of the flow field under the control of Reynolds equation and improved Reynolds equation on the monitoring surface.

[0058] Figure 6 This is the velocity distribution cloud diagram of the flow field under the control of Reynolds equation and improved Reynolds equation.

[0059] Figure 7 It is the characteristic diagram of flow velocity distribution in the flow field controlled by Reynolds equation and improved Reynolds equation on the monitoring surface.

[0060] Figure 8 Schematic diagram for comparing the flow forms of fractures under the control of Reynolds equation and improved Reynolds equation.

[0061] Fig. 9 Schematic diagram of the coupling relationship between various physical fields in the fracture dissolution model. DETAILED DESCRIPTION

[0062] In order to make the purpose, technical solution and advantages of the embodiments of the present application clearer, the embodiments of the present application will be described in detail below in conjunction with the accompanying drawings. The terms used in the implementation method section of the present application are only used to explain the specific embodiments of the present application, and are not intended to limit the present application.

[0063] Before introducing the specific solutions provided by the embodiments of the present application, some terms in the present application are explained to facilitate understanding by those skilled in the art, and the terms in the present application are not limited.

[0064] Three-dimensional rough fractures refer to cracks or fissures described and simulated in three-dimensional space, whose surface has a complex and irregular rough structure and has more significant seepage characteristics compared with two-dimensional fractures or smooth fractures.

[0065] Nonlinear flow field means that the fluid flow no longer conforms to a simple linear relationship, that is, the relationship between flow velocity and pressure gradient is no longer linear. In rough fractures, fluid flow is affected by surface roughness, local eddies, complex boundary conditions, etc., and the flow performance may be complex and changeable, which cannot be simply described by classic linear fluid mechanics formulas (such as Darcy's law).

[0066] Unless otherwise specified, the ordinal numbers such as "first" and "second" mentioned in the embodiments of the present application are used to distinguish multiple objects and are not used to limit the order, timing, priority or importance of the multiple objects.

[0067] The embodiments of the present application are described below in conjunction with the accompanying drawings.

[0068] The embodiments of the present application can be applied to Figure 1 Among the electronic devices shown, the electronic device can be but is not limited to mobile terminals such as mobile phones, tablet computers, handheld computers, personal digital assistants (PDAs), smart home devices such as smart TVs and smart cameras, wearable devices such as smart bracelets, smart watches, smart glasses, or other computer devices such as desktops, laptops, notebook computers, ultra-mobile personal computers (UMPCs), netbooks, smart screens, etc.

[0069] like Figure 1 As shown, the electronic device 200 may include one or more of the following components: a processor 201, a memory 203, a communication interface 202, and a communication bus 204. The memory 203 may be connected to the processor 201 via the bus 204. The bus may transmit data between the processor 201 and the memory 203. The bus may be divided into an address bus, a data bus, a control bus, and the like.

[0070] The processor 201 may include one or more processing cores, and the processor 201 may use various interfaces and lines to connect various parts of the entire electronic device 200, and execute various functions and process data of the electronic device 200 by running or executing instructions, programs, code sets or instruction sets stored in the memory 203, and calling data stored in the memory 203. Exemplarily, the processor 201 may include an application processor (application processor, AP), a modem processor, a CPU, a graphics processor (graphics processing unit, GPU), an image signal processor (image signal processor, ISP), a controller, a video codec, a digital signal processor (digital signal processor, DSP), a field programmable gate array (field-programmable gate array, FPGA), a programmable logic array (programmable logic array, PLA) and / or a neural network processor (neural-network processing unit, NPU), etc. Among them, the CPU mainly processes the operating system, user interface and application programs, etc.; the GPU is responsible for rendering and drawing the content to be displayed; the NPU is used to implement artificial intelligence (artificial intelligence, AI) functions; and the modem is used to process wireless communications. Different processing units may be independent devices or integrated into one or more processors. For example, the multiple processing units shown above are integrated into a SoC, or the AP is a separate semiconductor chip and other processing units are integrated into a SoC, which is not limited in this application.

[0071] The memory 203 may include a random access memory (RAM), a read-only memory (ROM), or a non-transitory computer-readable storage medium. The memory 203 may be used to store instructions, programs, codes, code sets, or instruction sets. The memory 203 may include a program storage area and a data storage area, wherein the program storage area may store instructions for implementing an operating system and instructions for at least one function, such as a three-dimensional rough fracture nonlinear flow field simulation method, a three-dimensional rough fracture long-term dissolution mechanism analysis method, etc.; the data storage area may store data created according to the use of the electronic device 200, such as input data for flow field simulation, etc.

[0072] In addition, those skilled in the art can understand that the structure of the electronic device 200 shown in the above drawings does not constitute a limitation on the electronic device 200, and the electronic device may include more or fewer components than shown in the drawings, or combine certain components, or arrange the components differently. For example, the electronic device 200 also includes components such as a microphone, a speaker, a radio frequency circuit, a sensor, an audio circuit, a power supply, and a Bluetooth module, which will not be described in detail here.

[0073] The embodiment of the present application provides a three-dimensional rough fracture nonlinear flow field simulation method for realizing the simulation of nonlinear flow through rough fractures under transient conditions and accurately reflecting the actual flow capacity of rough fractures.

[0074] Figure 2 A schematic flow chart of a three-dimensional rough fracture nonlinear flow field simulation method provided in this embodiment.

[0075] like Figure 2 As shown in Figure 2, the three-dimensional rough fracture nonlinear flow field simulation method includes:

[0076] Step S101: generating an initial geometric opening field of the three-dimensional rough crack based on the three-dimensional coordinates of the upper wall surface of the three-dimensional rough crack obtained by three-dimensional laser scanning technology.

[0077] Step S102: Use the Reynolds equation with time-varying terms of opening to describe the hydraulic head distribution of the linear flow field in the fracture, embed the Forchheimer equation into the Reynolds equation with time-varying terms of opening, and obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows:

[0078] (1)

[0079] In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = Ts ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )for t time( x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 .

[0080] To facilitate understanding of step S102, the Reynolds equation and the Forchheimer equation considering the time-varying property of the opening are respectively described in detail below.

[0081] 1) Flow field control equation - Reynolds equation considering the time-varying opening

[0082] In existing research, the Navier-Stokes equations (NS equations) are usually used to describe the flow of fluids in fractures. However, the NS equations are relatively complex and difficult to couple with other physical fields. In order to simplify the flow field control equations, the convection terms in the NS equations are usually ignored, resulting in a simpler form of the Stokes equations; further, by ignoring the momentum component perpendicular to the fracture surface, the Stokes equations can be simplified to the Reynolds equations. The Reynolds equation is a "depth averaging" method that can reduce the dimension of the solution model and is often used for approximate calculations of the flow field. By adding the time-varying term of the opening as a source term on the right side of the Reynolds equation, the Reynolds equation that takes into account the time-varying nature of the opening can be obtained:

[0083] (2)

[0084] In the formula, T s = T s ( x , y , t )for t time( x , y ) at the hydraulic conductivity (m 2 / s), H = H ( x , y , t )for t time( x , y ) at the water head (m), b = b ( x , y , t )for t time( x , y ) is the crack opening (m).

[0085] The Reynolds equation is also called the local cubic law. The permeability of each crack is T s Satisfies the cubic law:

[0086] (3)

[0087] In the formula, ρ is the fluid density (kg / m 3 ), g is the acceleration due to gravity (m / s 2 ), μis the dynamic viscosity of the fluid (Pa·s).

[0088] The essence of the Reynolds equation (Equations (2) and (3), ignoring the source term on the right side of equation (2)) is to treat each point in the fracture as a parallel plate. The relationship between the fluid volume flow rate and the pressure gradient in each local fracture satisfies the cubic law. The head is assigned to each point in the fracture by assembling the overall hydraulic conductivity matrix. Therefore, the Reynolds equation is also called the local cubic law. Since the flow capacity of each fracture satisfies a linear relationship, the overall flow capacity of the flow field obtained by solving the Reynolds equation is still a linear relationship. Many studies have shown that when the fluid flow rate in the fracture increases, the inertia effect increases, and the fluid volume flow rate and pressure gradient change from a linear relationship to a nonlinear relationship. Therefore, the Reynolds equation is not sufficient to fully describe the actual distribution characteristics of the flow field in the fracture.

[0089] 2) Forchheimer equation

[0090] The nonlinear relationship between pressure gradient and volume flow rate can be described by the Forchheimer equation:

[0091] (4)

[0092] In the formula P is the pressure (Pa), Q is the volume flow rate (m 3 / s), A is the linear coefficient characterizing the crack flow capacity (kg / (m 5 s)), B is a nonlinear term mainly related to the crack surface roughness (kg / m 8 ).

[0093] P They can be expressed as:

[0094] (5)

[0095] A and B They can be expressed as:

[0096] (6)

[0097] (7)

[0098] In the formula β is the coefficient related to roughness (m -1 ), the greater the roughness β The bigger; w is the crack width (m), eh is the fracture hydraulic aperture (m). Many studies have established the fracture hydraulic aperture e h and mechanical opening e m experience relationship.

[0099] In order to fully describe the actual distribution characteristics of the flow field in the fracture, this embodiment embeds the Forchheimer equation into the Reynolds equation with a time-varying term of the opening, and obtains the improved Reynolds equation (i.e., equation (1)) for describing the nonlinear flow field.

[0100] Step S103: Based on the preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture, solve the water head distribution of the nonlinear flow field under the control of the improved Reynolds equation.

[0101] In the embodiment of the present application, a three-dimensional rough surface is acquired by a three-dimensional laser scanning method to generate a rough crack space. By embedding the nonlinear Forchheimer equation into the Reynolds equation with a time-varying term of the opening, the simulation of nonlinear flow in rough cracks under transient conditions is achieved, laying a foundation for conducting research on the long-term dissolution mechanism of three-dimensional rough cracks.

[0102] Optionally, in step S101, based on the three-dimensional coordinates of the upper side wall of the three-dimensional rough crack obtained by three-dimensional laser scanning technology, an initial geometric opening field of the three-dimensional rough crack is generated, including:

[0103] Step S111: Import the three-dimensional coordinates of the upper wall surface of the three-dimensional rough fracture as data points into MATLAB to construct an initial data matrix.

[0104] Step S121: Generate a three-dimensional simulation grid according to the size of the three-dimensional rough fracture, and perform interpolation processing on the initial data matrix to generate a geometric opening simulation surface of the upper side wall of the three-dimensional rough fracture.

[0105] Step S131: the geometric opening simulation surface of the upper wall of the three-dimensional rough fracture is moved horizontally and vertically to simulate fracture shear sliding, so as to obtain the geometric opening simulation surface of the lower wall of the three-dimensional rough fracture.

[0106] Step S141: Calculate the geometric opening difference of the upper and lower wall surfaces and the opening difference of each grid in the simulated curved surface to obtain the initial geometric opening field of the three-dimensional rough fracture.

[0107] Analyzing equations (4), (6) and (7), we can see that due to the mechanical opening of the crack e m It can be determined by three-dimensional laser scanning and other methods. Therefore, according to equations (6) and (7), A andB is a constant and has nothing to do with the size of each local opening in the fracture. Therefore, although equations (6) and (7) can satisfy the nonlinear relationship between the overall volume flow rate and the average pressure gradient, the head and flow rate distribution in the fracture do not take into account the actual distribution characteristics of each opening.

[0108] In view of this, this embodiment modifies the "local cubic law" of rough cracks into the "local Forchheimer" law, that is, the Forchheimer equation is embedded in the Reynolds equation with a time-varying term of the aperture to obtain the improved Reynolds equation coupled with the Forchheimer equation (abbreviated as: improved Reynolds equation).

[0109] The specific steps of embedding the Forchheimer equation into the Reynolds equation with time-varying terms of the opening are described in detail below, including:

[0110] Step S112: Use the crack geometry opening b Alternative fracture hydraulic aperture e h , used to solve the linear coefficients in the Forchheimer equation A and nonlinear coefficient B , after replacement, the linear coefficient A and nonlinear coefficient B The expression is as follows:

[0111] (8)

[0112] (9)

[0113] In the formula, β is a coefficient related to roughness, the unit is m -1 ; μ is the fluid dynamic viscosity, in Pa·s; due to the geometric opening of the crack b changes over time, so the coefficient A , B also changes over time, i.e., A = A ( x , y , t ), B = B ( x , y , t ), the pressure distribution in the cracks is related to the size of each opening and changes continuously with dissolution.

[0114] Step S122: Replace the linear coefficients Aand nonlinear coefficient B Substituting into the Forchheimer equation, that is, combining equations (4), (5), (8), (9) and the fluid volume flow rate being greater than 0, the fluid volume flow rate is obtained Q The first expression is:

[0115] (10)

[0116] Step S132: Combined with existing research, fluid volume flow rate Q The second expression is as follows:

[0117] (11)

[0118] Step S142: Using the fluid volume flow rate Q The first expression and fluid volume flow rate Q The second expression of t time( x , y ) T s ;

[0119] Step S152: The direction perpendicular to the fracture surface is the “depth average” direction. According to the local cubic law (Equation (3)), the flow velocity in the opening direction is constant, and the flow velocity in the flow field is U =( u , v ) is as follows:

[0120] (12)

[0121] In the formula, u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s;

[0122] Step S162: t time( x , y ) T s Substitute the velocity into the flow field U =(u , v ), we can obtain the velocity distribution of the flow field under the influence of nonlinear seepage in fractures, that is, the improved Reynolds equation (Equation (1)) for describing the nonlinear flow field.

[0123] In order to verify the rationality of the improved Reynolds equation proposed in this application, the flow fields controlled by the Reynolds equation under steady-state conditions and the improved Reynolds equation are compared. The specific steps and analysis results are as follows:

[0124] First, steps S111 to S141 are written into computer code to complete the three-dimensional rough crack modeling. For example, the crack size can be set to 80 mm × 80 mm, and the three-dimensional morphology structure of the upper side wall of the crack can be obtained by three-dimensional laser scanning technology. The result is as follows: Figure 3 As shown in (a), the upper wall is then moved horizontally and vertically to simulate the fracture shear sliding, thereby forming the other side wall (i.e., the lower wall). The results are shown in Figure 3 As shown in (b), the opening field (i.e., the crack opening cloud map) is obtained by the difference in opening between the upper and lower walls. The result is shown in Figure 3 As shown in (c), the average opening is calculated to be about 1.97 mm.

[0125] against Figure 3 The three-dimensional rough fracture model shown in the figure has the following boundary conditions: the upper and lower sides are impermeable boundaries; the left side is the Dirichlet boundary condition, and the water head H 0 =0.02m; the right side is also the Dirichlet boundary condition, water head H 0 =0.0m. The density and dynamic viscosity of water are 1000kg / m 3 and 0.001Pa·s.

[0126] In addition, in order to quantitatively compare the distribution characteristics of field variables under the two conditions, we set x =40mm, y =40mm two monitoring surfaces.

[0127] Based on the above settings, the improved Reynolds equation under steady-state conditions and the head distribution characteristics of the flow field controlled by the Reynolds equation are solved respectively.

[0128] Exemplarily, the model can be solved using multi-physics coupling software COMSOL Multiphysics (version 5.6), or other software or code can be used to solve it. This application does not limit the specific solution method.

[0129] Figure 4 The hydraulic head distribution characteristics of the flow field under the control of the Reynolds equation and the improved Reynolds equation are described. (a) is the hydraulic head distribution characteristic diagram of the flow field under the control of the Reynolds equation (unit: mm); (b) is the hydraulic head distribution characteristic diagram of the flow field under the control of the improved Reynolds equation (unit: mm). It can be seen that the hydraulic head distribution in the flow field under the two situations is relatively close.

[0130] Figure 5 This is a schematic diagram of the water head distribution of the flow field under the control of the Reynolds equation and the improved Reynolds equation on the monitoring surface. In the figure, the red line is x =40mm monitoring surface, blue line is y =40mm monitoring surface, from Figure 5 The results on the monitoring surface also show that the water head distribution in the flow direction and perpendicular to the flow direction is relatively close. Therefore, the modified Reynolds equation has little effect on the distribution of water head in the fracture.

[0131] Figure 6 The velocity distribution cloud diagram of the flow field under the control of Reynolds equation and improved Reynolds equation. In the figure, (a) is the velocity distribution cloud diagram of the flow field under the control of Reynolds equation (unit: m / s); (b) is the velocity distribution cloud diagram of the flow field under the control of improved Reynolds equation (unit: m / s). Figure 6 It can be seen that the overall characteristics of the velocity distribution cloud maps of the two are relatively similar, the relative size distribution areas of the velocity are basically consistent, and the maximum velocity of the flow field controlled by the improved Reynolds equation is slightly smaller than the maximum velocity in the flow field controlled by the Reynolds equation.

[0132] Figure 7 The velocity distribution characteristics of the flow field under the control of the Reynolds equation and the improved Reynolds equation on the monitoring surface. Among them, (a) is x =40mm cross-sectional flow velocity comparison diagram, (b) y =40mm cross-sectional flow velocity comparison diagram. x =40mm cross section and y =40mm section flow rate monitoring results show that the improved Reynolds equation considering the influence of nonlinear seepage reduces the fluid flow rate in the fracture, and the flow capacity of the fracture decreases, which is also in line with the general law of three-dimensional rough fracture seepage. Therefore, the improved Reynolds equation proposed in this application effectively reflects the influence of rough protrusions in the fracture on the seepage process.

[0133] Figure 8 Schematic diagram for comparing the flow forms of fractures under the control of Reynolds equation and improved Reynolds equation. Figure 8 (a) shows that the pressure gradient of the flow field under the control of the Reynolds equation is linearly related to the volume flow rate, while Figure 8 (b) shows that the pressure gradient and volume flow rate of the flow field under the control of the improved Reynolds equation are nonlinear. Figure 8 (b) also shows that when the crack is rougher, the nonlinear coefficient β The larger the value, the stronger the nonlinearity of pressure gradient and volume flow rate, which is consistent with the general characteristics of rough fracture seepage. Therefore, the improved Reynolds equation of this application can better reflect the nonlinear seepage characteristics of fluid in three-dimensional rough rock fractures.

[0134] The present application provides a method for analyzing the long-term dissolution mechanism of three-dimensional rough cracks, comprising the following steps:

[0135] Calculating the flow field of the three-dimensional rough fracture under steady-state conditions using the linear flow field simulation method provided in any of the above embodiments, and calculating the concentration field based on the depth-averaged solute transport equation to obtain the initial distribution of the flow field and the concentration field;

[0136] Through the one-way cycle coupling of flow field, aperture field and concentration field, a fracture dissolution coupling model is established;

[0137] Analyze the difference Da Number and sum P The evolution characteristics of the aperture field, flow field, concentration field, and roughness of three-dimensional rough fractures under the action of number and long-term dissolution are analyzed to reveal the long-term dissolution characteristics and dissolution mechanism of three-dimensional rough fractures.

[0138] The process of using the linear flow field simulation method to calculate the flow field of a three-dimensional rough fracture under steady-state conditions can be performed with reference to the steps of the linear flow field simulation method provided in any of the above embodiments, which will not be described in detail here.

[0139] The calculation process of the concentration field is introduced below.

[0140] The solute in the crack produces a convection effect under the flow field; at the same time, the dissolved minerals diffuse. Some scholars assume that the concentration does not change much along the opening direction and propose a solute transport equation based on depth average:

[0141] (13)

[0142] In the formula, c = c ( x , y , t )for t time( x , y ) solute concentration (mol / m 3 ),D is the local diffusion tensor (m 2 / s), under isotropic conditions D No longer a matrix, use D express; R ( c )= R ( c ( x , y , t )) is the source term, representing the local reaction rate (mol / (m 2 ·s)), which is related to the concentration of the solution. R ( c ) is as follows:

[0143] (14)

[0144] In the formula, k eff is the first-order dynamic response coefficient (m / s), C s is the solute saturation concentration (kg / m 3 ).

[0145] At the same time, when considering surface reactions and solute migration, k eff It can be expressed as:

[0146] (15)

[0147] In the formula, k is the surface reaction rate coefficient (m / s); Sh is the dimensionless Sherwood number, which is 7.56 for parallel plates.

[0148] After the minerals on the surface of the fracture dissolve and migrate, the fracture opening gradually increases due to dissolution. The opening increment can be expressed by the following control equation:

[0149] (16)

[0150] In the formula, M is the solid phase molar mass (kg / mol), ρ s is the solid phase density (kg / m 3 ). Formula (16) is also the control equation of the opening change caused by crack dissolution.

[0151] The fracture dissolution coupling model was solved using the custom partial differential equation interface in the multi-physics coupling software COMSOL Multiphysics (version 5.6).

[0152] First, the initial geometric opening field is imported through the interpolation function. b 0 It can be expressed as:

[0153] (17)

[0154] In the formula, data( x , y ) is the three-dimensional roughness field interpolation function (mm), representing the coordinates ( x , y ) at the initial local opening. Each initial opening b 0 is a constant.

[0155] The change of each opening with time is recorded as b pp = b pp ( x , y , t )(mm), then the opening at each moment b = b ( x , y , t ) can be written as:

[0156] (18)

[0157] Formula (1) gives an improved form of the flow field control equation, namely the improved Reynolds equation. The opening change control equation shown in formula (16) is also relatively easy to solve, so the key to the problem lies in solving the concentration field solute transport control equation (formula (13)). This embodiment uses the general form of the stable convection-diffusion equation to solve formula (13). The general form of the stable convection-diffusion equation is:

[0158] (19)

[0159] In the formula d 1 , d 5 and d 6 are the damping coefficient (m), absorption coefficient (m / s) and source term (mol / (m 2 ·s)); d 2 is a vector representing the diffusion coefficient (assuming isotropic diffusion, m 3 / s); d 3 is a vector representing the conserved flux convection coefficient (m2 / s); d 4 is a vector representing the convection coefficient (m 2 / s).

[0160] It can be seen that there is a big difference between formula (13) and formula (19). Therefore, it is necessary to sort out formula (13) first so that it is consistent with formula (19) in form:

[0161] (20)

[0162] Comparing equation (19) and equation (20), we have:

[0163] (twenty one)

[0164] The calculation of the fracture dissolution coupling model is divided into two steps: first, the flow field and concentration field under steady-state conditions are calculated to obtain the initial distribution of the two fields; then, the coupling of the three fields under transient conditions is calculated. Among them, the flow field, aperture field, and concentration field are one-way coupling relationships, such as Fig. 9 As shown. The field variables (or intermediate variables) obtained by solving each field are substituted into the next field to solve the unknown field variables, thereby realizing the coupling of the model. Specifically, using the given opening field b and the total head difference, the distribution of the head H and the flow velocity U in the flow field can be calculated; substituting the U solved in the flow field into the concentration field, the field variable c can be calculated; then substituting the concentration field variable c into the opening field b, the change in the solid phase structure can be calculated. At the same time, the change in b will further affect the distribution of the head H and the flow velocity U in the flow field, thus forming a one-way cyclic coupling process.

[0165] Furthermore, in order to further verify the rationality of the fissure dissolution model, a three-dimensional parallel plate is first used for verification. By constructing a geometric model of a three-dimensional parallel plate fissure, setting relevant parameters, respectively solving the analytical solution and numerical solution of the concentration field, and performing comparative analysis, the results show that under different dissolution conditions, the solute concentration distribution calculated by the analytical solution is in good agreement with the numerical solution, which indicates that the fissure dissolution model proposed in this embodiment is reasonable. Secondly, the numerical solution and analytical solution of the fissure dissolution model under one-dimensional conditions are compared, and the curves of the time-varying characteristics of the opening of each fissure are analyzed and drawn. The results show that the distance from the entrance at different times and the opening curve, the theoretical solution of the entrance opening and the numerical solution of the entrance opening at different times are in good agreement, which also indicates that the fissure dissolution model proposed in this embodiment is reasonable.

[0166] On the basis of the rationality of the fracture dissolution model, different Da Number and sum PThe evolution characteristics of the aperture field, flow field, concentration field, and roughness of three-dimensional rough fractures under the action of several years and long-term dissolution (for example, up to 12 years) are studied to reveal the long-term dissolution characteristics and dissolution mechanism of three-dimensional rough fractures. This fills the gap in the existing research on the long-term dissolution mechanism and dissolution characteristics of rough fractures to a certain extent.

[0167] The present application provides a three-dimensional rough fracture nonlinear flow field simulation system, comprising:

[0168] A simulation unit configured to generate an initial geometric opening field of the three-dimensional rough fracture based on the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture obtained by the three-dimensional laser scanning technology;

[0169] The improved unit is configured to use the Reynolds equation with a time-varying term of opening to describe the water head distribution in the fracture, and the Forchheimer equation is embedded in the Reynolds equation with a time-varying term of opening to obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows:

[0170] ,

[0171] In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = T s ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )fort time( x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 ;

[0172] The solving unit is configured to solve the head distribution of the nonlinear flow field under the control of the improved Reynolds equation based on preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture.

[0173] The present application also provides a computer program product, which includes computer executable instructions. In one embodiment, the computer executable instructions are used to enable a computer to execute the functions of the above method embodiment.

[0174] The computer executable instructions may be stored in a computer readable storage medium. The present application also provides a computer readable storage medium, wherein the computer readable storage medium stores the executable instructions. In one embodiment, the computer executable instructions are used to enable a computer to execute the functions in the above method embodiment.

[0175] The computer-readable storage medium provided in the embodiments of the present application may be a random access memory (RAM), a flash memory, a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), a register, a hard disk, a mobile hard disk, a CD-ROM, or any other form of computer-readable storage medium known in the art.

[0176] Computer executable instructions can be stored in a computer readable storage medium or transferred from one computer readable storage medium to another computer readable storage medium. For example, the computer program or instructions can be transferred from one website, computer, server or data center to another website, computer, server or data center by wired or wireless means. The computer readable storage medium can be any available medium that can be accessed by a computer or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium, such as a floppy disk, a hard disk, or a magnetic tape; it can also be an optical medium, such as a digital video disc (DVD); it can also be a semiconductor medium, such as a solid state drive.

[0177] The above description is only a preferred embodiment of the present application and is not intended to limit the present application. For those skilled in the art, the present application may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A three-dimensional rough fracture nonlinear flow field simulation method, characterized in that: include: Based on the three-dimensional coordinates of the upper wall of the three-dimensional rough crack obtained by three-dimensional laser scanning technology, the initial geometric opening field of the three-dimensional rough crack is generated; The Reynolds equation with a time-varying term of the opening is used to describe the head distribution of the linear flow field in the fracture. The Forchheimer equation is embedded in the Reynolds equation with a time-varying term of the opening to obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows: , In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = T s ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )for t time( x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 ; Based on the preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture, solving the head distribution of the nonlinear flow field under the control of the improved Reynolds equation; The Forchheimer equation is embedded into the Reynolds equation with the time-varying term of the opening, including: Use fracture geometry b Alternative fracture hydraulic aperture e h , used to solve the linear coefficients in the Forchheimer equation A and nonlinear coefficient B , and then obtain the linear coefficient after replacement A and nonlinear coefficient B Expression of The linear coefficients after replacement A and nonlinear coefficient B Substituting into the Forchheimer equation, we get the fluid volume flow rate Q The first expression of Combined with existing research, the fluid volume flow rate is obtained Q The second expression of Using the fluid volume flow Q The first expression and the fluid volume flow rate Q The second expression of t time( x , y ) T s ; According to the local cubic law, the flow velocity in the opening direction is constant, and the flow velocity in the flow field is determined U =( u , v ) In the crack t time( x , y ) T s Substitute the velocity into the flow field U =( u , v ) to obtain the improved Reynolds equation.

2. The three-dimensional rough fracture nonlinear flow field simulation method according to claim 1, characterized in that: Based on the three-dimensional coordinates of the upper wall of the three-dimensional rough fracture obtained by three-dimensional laser scanning technology, the initial geometric opening field of the three-dimensional rough fracture is generated, including: Importing the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture into MATLAB as data points to construct an initial data matrix; Generating a three-dimensional simulation grid according to the size of the three-dimensional rough crack, and performing interpolation processing on the initial data matrix to generate a geometric opening simulation surface of the upper side wall of the three-dimensional rough crack; The geometric opening simulation surface of the upper side wall of the three-dimensional rough fracture is moved horizontally and vertically to simulate fracture shear sliding, thereby obtaining the geometric opening simulation surface of the lower side wall of the three-dimensional rough fracture; The geometric opening of the upper and lower walls is calculated to simulate the opening difference at each grid in the curved surface, and the initial geometric opening field of the three-dimensional rough crack is obtained.

3. A method for analyzing the long-term dissolution mechanism of three-dimensional rough fractures, characterized in that: include: The three-dimensional rough fracture nonlinear flow field simulation method according to any one of claims 1 to 2 is used to calculate the flow field of the three-dimensional rough fracture under steady-state conditions, and the concentration field is calculated based on the depth-averaged solute transport equation to obtain the initial distribution of the flow field and the concentration field; Through the one-way cycle coupling of flow field, aperture field and concentration field, a fracture dissolution coupling model is established; Analyze the difference Da Number and P The evolution characteristics of the aperture field, flow field, concentration field, and roughness of three-dimensional rough fractures under the action of number and long-term dissolution are analyzed to reveal the long-term dissolution characteristics and dissolution mechanism of three-dimensional rough fractures.

4. A three-dimensional rough fracture nonlinear flow field simulation system, characterized in that: include: A simulation unit configured to generate an initial geometric opening field of the three-dimensional rough fracture based on the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture obtained by the three-dimensional laser scanning technology; The improved unit is configured to use the Reynolds equation with a time-varying term of opening to describe the water head distribution in the fracture, and the Forchheimer equation is embedded in the Reynolds equation with a time-varying term of opening to obtain an improved Reynolds equation for describing the nonlinear flow field. The expression of the improved Reynolds equation is as follows: , In the formula, U is the flow velocity in the flow field, U =( u , v ), u = u ( x , y , t )and v = v ( x , y , t ) are respectively t time( x , y ) X Axis direction and Y Axial velocity component, in m / s; T s = T s ( x , y , t )for t time( x , y ) at the point where the hydraulic conductivity is in m 2 / s; H = H ( x , y , t )for t time( x , y ) at the water head, in meters, b = b ( x , y , t )for t time( x , y ) is the geometric opening of the crack at the position of m; ρ is the fluid density in kg / m 3 , g is the acceleration due to gravity, in m / s 2 ; w is the crack width, in m; A is the linear coefficient, which represents the flow capacity of the crack, and the unit is kg / (m 5 s); B is the nonlinear coefficient related to the crack surface roughness, the unit is kg / m 8 ; A solving unit configured to solve the head distribution of the nonlinear flow field under the control of the improved Reynolds equation based on preset boundary conditions and the initial geometric opening field of the three-dimensional rough fracture; The Forchheimer equation is embedded into the Reynolds equation with the time-varying term of the opening, including: Use fracture geometry b Alternative fracture hydraulic aperture e h , used to solve the linear coefficients in the Forchheimer equation A and nonlinear coefficient B , and then obtain the linear coefficient after replacement A and nonlinear coefficient B Expression of The linear coefficients after replacement A and nonlinear coefficient B Substituting into the Forchheimer equation, we get the fluid volume flow rate Q The first expression of Combined with existing research, the fluid volume flow rate is obtained Q The second expression of Using the fluid volume flow Q The first expression and the fluid volume flow rate Q The second expression of t time( x , y ) T s ; According to the local cubic law, the flow velocity in the opening direction is constant, and the flow velocity in the flow field is determined U =( u , v ) In the crack t time( x , y ) T s Substitute the velocity into the flow field U =( u , v ) to obtain the improved Reynolds equation.

5. The system according to claim 4, characterized in that The simulation unit is further configured as follows: Importing the three-dimensional coordinates of the upper side wall of the three-dimensional rough fracture into MATLAB as data points to construct an initial data matrix; Generating a three-dimensional simulation grid according to the size of the three-dimensional rough crack, and performing interpolation processing on the initial data matrix to generate a geometric opening simulation surface of the upper side wall of the three-dimensional rough crack; The geometric opening simulation surface of the upper side wall of the three-dimensional rough fracture is moved horizontally and vertically to simulate fracture shear sliding, thereby obtaining the geometric opening simulation surface of the lower side wall of the three-dimensional rough fracture; The geometric opening of the upper and lower walls is calculated to simulate the opening difference at each grid in the curved surface, and the initial geometric opening field of the three-dimensional rough crack is obtained.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the three-dimensional rough fracture nonlinear flow field simulation method as described in any one of claims 1 to 2 or the three-dimensional rough fracture long-term dissolution mechanism analysis method as described in claim 3 is implemented.

7. An electronic device, characterized in that: include: A memory, a processor, and a program stored in the memory and executable on the processor, wherein when the processor executes the program, the three-dimensional rough fracture nonlinear flow field simulation method as described in any one of claims 1 to 2 or the three-dimensional rough fracture long-term dissolution mechanism analysis method as described in claim 3 is implemented.

8. A computer program product, characterized in that The method comprises computer executable instructions for causing a computer to execute the three-dimensional rough fracture nonlinear flow field simulation method as described in any one of claims 1 to 2 or the three-dimensional rough fracture long-term dissolution mechanism analysis method as described in claim 3.

Images