Fracture and matrix grid-based non-darcy seepage stress coupling simulation method

Abstract

The invention relates to the technical field of engineering design, in particular to a fracture and matrix grid-based non-darcy seepage stress coupling simulation method, which comprises the following steps of: inputting information such as fracture and matrix unit grids, boundary conditions, material parameters and control parameters; giving an initial value and an initial stress of the basic unknown quantity; sequentially carrying out matrix seepage and embedded fracture non-darcy flow iterative solution and matrix stress balance equation solution at each time step; judging whether convergence is carried out or not, and if the time step calculation is not converged, continuing to carry out the next iteration solution; and after a convergence condition is reached, entering a next time step for solving until the time is over. The non-Darcy seepage stress coupling simulation method provided by the invention is high in precision and high in efficiency, the problems of high local structural plane modeling difficulty, low efficiency, poor effect and the like during traditional finite element solution are solved, and the non-Darcy seepage characteristics and mechanical deformation behaviors of the rock mass under high seepage pressure can be effectively simulated.

Description

Technical Field

[0001] The present application relates to the field of engineering design technology, and in particular to a non-Darcy seepage stress coupling simulation method based on fracture and matrix grids. Background Art

[0002] In water conservancy and hydropower projects, the rock mass in the dam area exists under unique hydrogeological conditions before water is stored. After the reservoir is filled with water, the high head has caused severe disturbances to the original hydrogeological environment due to the mutual feedback and coupling between water and rock, causing stress redistribution in the dam area. The dam foundation rock is a multiphase medium system with a complex structure, consisting of rock blocks and structural surfaces such as fractures, weak interlayers, and dislocation zones in the rock strata. The rock blocks are mainly rock matrix and primary or secondary pores, and the weak structural surfaces in the rock strata are mainly cracks. Under the action of high osmotic pressure, the flow of water in the rock mass cracks in the dam area no longer satisfies Darcy's seepage law, but presents non-Darcy seepage characteristics.

[0003] Traditional studies on seepage-stress coupling in fractured rock mass are mostly based on the equivalent continuum model. The equivalent continuum model has great advantages in macroscopic seepage analysis and is relatively simple to calculate. However, when the fracture scale is large, the fractured rock mass should be considered as a discontinuous medium. For example, in the patent CN118095122A, a model construction based on real rough fractures and a non-Darcy seepage simulation method are only for nonlinear seepage indicators of fractured rock mass with different rock types and fracture roughness, without considering discrete fractures and continuous matrix content.

[0004] The discrete fracture network model can better describe the heterogeneity and anisotropy of fractured rock mass. Its disadvantage is that it does not take into account the permeability of rock blocks, and it is difficult to fully obtain parameters such as the distribution, orientation and geometric characteristics of fractures in the rock mass, resulting in a very large workload for fracture generation, unit segmentation and numerical calculation. In fact, fractured rock mass is composed of discrete fractures and continuous matrix. The study of seepage stress coupling in fractured rock mass needs to consider both the hydraulic coupling process of rock blocks and the hydraulic properties of fractures, especially under the action of high head, the seepage of rock mass fractures will show obvious nonlinear characteristics. The advantage of the embedded fracture network model is that it realizes coupled iterative calculation by defining the flow exchange between the calculation fractures and the matrix, which greatly improves the calculation efficiency, can better describe the seepage and mechanical behavior of fractured rock mass, and more conveniently characterizes the non-Darcy preferential flow phenomenon of rock mass fractures. It has developed rapidly in recent years and has broad application prospects.

[0005] The present invention regards the fractured rock mass as a mixed system consisting of matrix and natural fractures, and the fracture water flow and matrix seepage as relatively independent systems. The coupled motion of the matrix and fracture water flow is realized through matrix-fracture flow exchange. Based on the basic principle of the finite element method and using an embedded fracture-matrix system model, a numerical solution model and method for fracture-matrix mixed media are developed. Summary of the invention

[0006] To achieve the above objectives, this application provides the following technical solutions:

[0007] According to a first aspect of the present invention, the present invention claims a non-Darcy seepage stress coupling simulation method based on fracture and matrix grid, comprising:

[0008] Inputting a simulated matrix grid, establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining the intersection reference information of the mesh units of the simulated fracture grid and the simulated matrix grid;

[0009] Selecting initial water pressure and initial displacement values ​​of grid nodes according to the simulated fracture grid and the simulated matrix grid;

[0010] The time-step iterative solution of fracture non-Darcy seepage and matrix seepage is used to obtain the water pressure tolerance and judge the convergence of the time-step iterative calculation;

[0011] According to the convergence and time step information, the solution of non-Darcy seepage stress coupling is obtained.

[0012] Furthermore, the inputting of the simulated matrix grid, establishing the simulated fracture grid according to the attribute information of the fracture, and obtaining the intersection reference information of the mesh units of the simulated fracture grid and the simulated matrix grid also includes:

[0013] Inputting a simulation matrix grid, and obtaining a matrix grid unit and a matrix grid node according to the simulation matrix grid;

[0014] Establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining fracture grid units and fracture grid nodes according to the simulated fracture grid;

[0015] Obtaining intersection reference information of the matrix grid unit and the fracture grid unit;

[0016] The intersection reference information at least includes:

[0017] Boundary conditions, material parameters and control parameters.

[0018] Furthermore, the selecting of initial water pressure and initial displacement values ​​of grid nodes according to the simulated fracture grid and the simulated matrix grid also includes:

[0019] Selecting the initial fracture water pressure of the fracture grid node according to the simulated fracture grid;

[0020] The initial matrix water pressure and the initial matrix displacement of the matrix grid nodes are selected according to the simulation matrix grid.

[0021] Furthermore, the method of using the time-step iterative solution of fracture non-Darcy seepage and matrix seepage to obtain the water pressure tolerance and judge the convergence of the time-step iterative calculation also includes:

[0022] The seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system are solved in sequence;

[0023] By calculating the matrix and fracture flow exchange terms, the matrix and fracture seepage coupling solution is realized. By calculating the matrix water pressure and fracture water pressure tolerance, it is determined whether the time step iterative calculation converges.

[0024] Furthermore, the method of obtaining the solution result of non-Darcy seepage stress coupling based on the convergence and time step information also includes:

[0025] Use the displacement tolerance to determine whether the time step iterative calculation has converged. If so, proceed to the next time step calculation; otherwise, proceed to the previous step to continue the iterative calculation.

[0026] Adopt adaptive time step to determine whether the calculation of the previous time step needs to be repeated;

[0027] When the simulation calculation time is reached, the solution of non-Darcy seepage stress coupling is obtained.

[0028] Furthermore, when obtaining the intersection reference information of the matrix grid unit and the crack grid unit, it also includes:

[0029] Given a simulation calculation time, set the initial time step, first iteration tolerance, second iteration tolerance, third iteration tolerance, and the maximum number of iterations between two iteration steps.

[0030] Furthermore, when selecting the initial fracture water pressure of the fracture mesh nodes according to the simulated fracture mesh and selecting the initial matrix water pressure and initial matrix displacement of the matrix mesh nodes according to the simulated matrix mesh, it also includes:

[0031] Let the time step n = 0, select the matrix water pressure of the matrix grid nodes, the fracture water pressure of the fracture grid nodes and the displacement of the matrix grid nodes as the basic unknowns, and give the initial values ​​of the displacement, matrix water pressure, fracture water pressure and stress time step as the initial matrix displacement, initial matrix water pressure, initial fracture water pressure and initial stress time step.

[0032] Furthermore, when solving the seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system in sequence, it also includes:

[0033] Solve the fracture non-Darcy seepage equation, matrix Darcy seepage equation and matrix stress balance equation in sequence;

[0034] Set the iteration step to 0, and give the iteration step displacement, initial iteration matrix water pressure, initial iteration fracture water pressure and initial stress iteration value;

[0035] The method further includes: calculating the matrix and fracture flow exchange term to achieve matrix and fracture seepage coupling solution, and calculating the matrix water pressure and fracture water pressure tolerance to determine whether the time step iterative calculation converges.

[0036] Solve the matrix Darcy flow equation and calculate the pore water pressure;

[0037] Solve the non-Darcy flow equation in fractures and calculate the fracture water pressure;

[0038] Calculate matrix water pressure tolerance and fracture water pressure tolerance;

[0039] According to the relationship between the tolerance of the matrix water pressure and the fracture water pressure and the tolerance of the first iteration and the tolerance of the second iteration, it is selected to solve the stress balance equation or to continue to iteratively solve the fracture non-Darcy seepage and matrix seepage equations;

[0040] When solving the stress equilibrium equation, the displacement increment, strain increment and stress increment are calculated, the displacement vector, strain and stress tensor of the iteration step are updated, and the displacement tolerance is calculated.

[0041] Furthermore, the displacement tolerance is used to judge whether the time step iterative calculation has converged. If it has converged, it will go to the next time step for calculation. Otherwise, it will go to the previous step for iterative calculation. It also includes:

[0042] If the displacement tolerance is less than or equal to the third iteration tolerance, proceed to the next time step calculation;

[0043] If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has not reached the maximum value, the number of iteration steps is incremented by 1, and the process goes back to the previous step to continue the iterative calculation.

[0044] Furthermore, the method of using an adaptive time step to determine whether the calculation of the previous time step needs to be repeated also includes:

[0045] If the number of iterations is less than half of the maximum number of iterations, increase the time step;

[0046] If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has exceeded the maximum value, the time step is reduced and the calculation of the previous time step is repeated.

[0047] The present application relates to the field of engineering design technology, and in particular to a non-Darcy seepage stress coupling simulation method based on fracture and matrix grids, inputting information such as fracture and matrix unit grids, boundary conditions, material parameters and control parameters; giving the initial value and initial stress of the basic unknown quantity; performing iterative solutions of matrix seepage and embedded fracture non-Darcy flow in turn at each time step, and solving the matrix stress equilibrium equation; judging whether it converges, and continuing to the next iterative solution if the time step calculation does not converge; after reaching the convergence condition, entering the next time step solution until the time ends. The present invention provides a high-precision and high-efficiency non-Darcy seepage stress coupling simulation method, which improves the problems of difficulty, low efficiency and poor effect in modeling local structural surfaces during traditional finite element solution, and can effectively simulate the non-Darcy seepage characteristics and mechanical deformation behavior of rock mass under high seepage pressure. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a workflow diagram of the non-Darcy seepage stress coupling simulation method based on fracture and matrix grid claimed in the embodiment of the present application;

[0049] Figure 2 It is a schematic diagram of an embedded fracture network model of a non-Darcy seepage stress coupling simulation method based on fractures and matrix grids claimed in the embodiments of the present application;

[0050] Figure 3 It is a matrix and fracture grid diagram of the non-Darcy seepage stress coupling simulation method based on fracture and matrix grids claimed in the embodiment of the present application;

[0051] Figure 4 It is a water pressure distribution diagram of a fractured rock mass according to a non-Darcy seepage stress coupling simulation method based on fracture and matrix grids claimed in the embodiments of the present application;

[0052] Figure 5 It is a horizontal displacement distribution diagram of a fractured rock mass according to a non-Darcy seepage stress coupling simulation method based on fracture and matrix grids as claimed in the embodiments of the present application. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0054] The terms "first", "second" and "third" in this application are only used for descriptive purposes, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of indicated technical features. Thus, the features defined as "first", "second" and "third" can explicitly or implicitly include at least one of the features. In the description of this application, the meaning of "multiple" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined. All directional indications in the embodiments of the present application (such as up, down, left, right, front, back...) are only used to explain the relative position relationship, movement, etc. between the components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication also changes accordingly. In addition, the terms "including" and "having" and any of their variations are intended to cover non-exclusive inclusions. For example, a process, method, system, product or device that includes a series of steps or units is not limited to the listed steps or units, but optionally also includes steps or units that are not listed, or optionally also includes other steps or units inherent to these processes, methods, products or devices.

[0055] Reference to "embodiments" herein means that a particular feature, structure, or characteristic described in conjunction with the embodiments may be included in at least one embodiment of the present application. The appearance of the phrase in various locations in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment that is mutually exclusive with other embodiments. It is explicitly and implicitly understood by those skilled in the art that the embodiments described herein may be combined with other embodiments.

[0056] According to the first embodiment of the present invention, a high-precision and high-efficiency non-Darcy seepage stress coupling simulation method based on fracture and matrix grids is provided, which improves the problems of difficulty, low efficiency and poor effect in modeling local structural surfaces in traditional finite element solution, and can effectively simulate the non-Darcy seepage characteristics and mechanical deformation behavior of rock mass under high permeability pressure.

[0057] Reference Figure 1 , the method comprising:

[0058] Inputting a simulated matrix grid, establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining the intersection reference information of the mesh units of the simulated fracture grid and the simulated matrix grid;

[0059] Selecting initial water pressure and initial displacement values ​​of grid nodes according to the simulated fracture grid and the simulated matrix grid;

[0060] The time-step iterative solution of fracture non-Darcy seepage and matrix seepage is used to obtain the water pressure tolerance and judge the convergence of the time-step iterative calculation;

[0061] According to the convergence and time step information, the solution of non-Darcy seepage stress coupling is obtained.

[0062] Furthermore, the inputting of the simulated matrix grid, establishing the simulated fracture grid according to the attribute information of the fracture, and obtaining the intersection reference information of the mesh units of the simulated fracture grid and the simulated matrix grid also includes:

[0063] Inputting a simulation matrix grid, and obtaining a matrix grid unit and a matrix grid node according to the simulation matrix grid;

[0064] Establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining fracture grid units and fracture grid nodes according to the simulated fracture grid;

[0065] Obtaining intersection reference information of the matrix grid unit and the fracture grid unit;

[0066] The intersection reference information at least includes:

[0067] Boundary conditions, material parameters and control parameters.

[0068] Furthermore, the selecting of initial water pressure and initial displacement values ​​of grid nodes according to the simulated fracture grid and the simulated matrix grid also includes:

[0069] Selecting the initial fracture water pressure of the fracture grid node according to the simulated fracture grid;

[0070] The initial matrix water pressure and the initial matrix displacement of the matrix grid nodes are selected according to the simulation matrix grid.

[0071] Furthermore, the method of using the time-step iterative solution of fracture non-Darcy seepage and matrix seepage to obtain the water pressure tolerance and judge the convergence of the time-step iterative calculation also includes:

[0072] The seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system are solved in sequence;

[0073] By calculating the matrix and fracture flow exchange terms, the matrix and fracture seepage coupling solution is realized. By calculating the matrix water pressure and fracture water pressure tolerance, it is determined whether the time step iterative calculation converges.

[0074] Furthermore, the method of obtaining the solution result of non-Darcy seepage stress coupling based on the convergence and time step information also includes:

[0075] Use the displacement tolerance to determine whether the time step iterative calculation has converged. If so, proceed to the next time step calculation; otherwise, proceed to the previous step to continue the iterative calculation.

[0076] Adopt adaptive time step to determine whether the calculation of the previous time step needs to be repeated;

[0077] When the simulation calculation time is reached, the solution of non-Darcy seepage stress coupling is obtained.

[0078] Furthermore, when obtaining the intersection reference information of the matrix grid unit and the crack grid unit, it also includes:

[0079] Given a simulation calculation time, set the initial time step, first iteration tolerance, second iteration tolerance, third iteration tolerance, and the maximum number of iterations between two iteration steps.

[0080] In this embodiment, the rock mass is meshed by finite element, the unit and node information of the matrix mesh is input, the intersection of the fracture and the matrix mesh is found according to the spatial position of the fracture, an embedded fracture finite element calculation mesh is established, the unit and node information of the fracture mesh is input, and the unit number of the intersection with the matrix mesh is recorded, such as Figure 2 shown.

[0081] Given the simulation calculation time t max , set the initial time step Δt, iteration tolerance Tol1, Tol2, Tol3 and the maximum number of iterations l for two iterations max , k max .

[0082] Furthermore, when selecting the initial fracture water pressure of the fracture mesh nodes according to the simulated fracture mesh and selecting the initial matrix water pressure and initial matrix displacement of the matrix mesh nodes according to the simulated matrix mesh, it also includes:

[0083] Let the time step n = 0, select the matrix water pressure of the matrix grid nodes, the fracture water pressure of the fracture grid nodes and the displacement of the matrix grid nodes as the basic unknowns, and give the initial values ​​of the displacement, matrix water pressure, fracture water pressure and stress time step as the initial matrix displacement, initial matrix water pressure, initial fracture water pressure and initial stress time step.

[0084] In this embodiment, let the time step n = 0, select the matrix water pressure of the matrix grid node, the fracture water pressure of the fracture grid node and the displacement of the matrix grid node as the basic unknown quantities, and give the initial matrix displacement, initial matrix water pressure, initial fracture water pressure and initial stress time step as u0, p respectively. m,0 、p f,0 ,σ′0;

[0085] Among them, m represents matrix and f represents fracture.

[0086] Furthermore, when solving the seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system in sequence, it also includes:

[0087] Solve the fracture non-Darcy seepage equation, matrix Darcy seepage equation and matrix stress balance equation in sequence;

[0088] Set the iteration step to 0, and give the iteration step displacement, initial iteration matrix water pressure, initial iteration fracture water pressure and initial stress iteration value;

[0089] The method further includes: calculating the matrix and fracture flow exchange term to achieve matrix and fracture seepage coupling solution, and calculating the matrix water pressure and fracture water pressure tolerance to determine whether the time step iterative calculation converges.

[0090] Solve the matrix Darcy flow equation and calculate the pore water pressure;

[0091] Solve the non-Darcy flow equation in fractures and calculate the fracture water pressure;

[0092] Calculate matrix water pressure tolerance and fracture water pressure tolerance;

[0093] According to the relationship between the tolerance of the matrix water pressure and the fracture water pressure and the tolerance of the first iteration and the tolerance of the second iteration, it is selected to solve the stress balance equation or to continue to iteratively solve the fracture non-Darcy seepage and matrix seepage equations;

[0094] When solving the stress equilibrium equation, the displacement increment, strain increment and stress increment are calculated, the displacement vector, strain and stress tensor of the iteration step are updated, and the displacement tolerance is calculated.

[0095] In this embodiment, the fracture-matrix system non-Darcy seepage stress coupling control equation includes the fracture non-Darcy seepage equation (1), the matrix seepage equation (2) and the matrix stress balance equation (3).

[0096]

[0097] In the formula, a f is the crack opening, k f is the matrix intrinsic permeability, J is the non-Darcy correction term, k m is the matrix intrinsic permeability, g is the gravitational acceleration, ρ w is the density of water, μ w is the viscosity coefficient of water, φ is the matrix porosity, K s is the compression modulus of rock particles, u is displacement, t is time, div is the divergence sign, c w is the volume compression coefficient of water, α is the Biot constant, k n is the intrinsic permeability of the matrix along the fracture normal direction, is the pressure gradient in the normal direction of the fracture, Among them, p m is the matrix pore water pressure, p f is the fracture water pressure, Ω mis the matrix domain, Γ f is the crack domain, σ is the total stress, ρ is the total density, ρ=φρ w +(1-φ)ρ s , ρ s is the solid particle density, b is the body force vector, n f is the crack normal direction vector.

[0098] 1) Let the iteration step l = 0, and give the iteration step displacement, initial iteration matrix water pressure, initial iteration fracture water pressure and initial stress iteration value;

[0099]

[0100] 2) The fracture non-Darcy seepage and matrix seepage are iteratively solved. The matrix seepage equation is solved first, and then the fracture non-Darcy seepage equation is solved. By accurately calculating the matrix and fracture flow exchange terms, the coupled solution of matrix and fracture seepage is achieved.

[0101] First solve the matrix Darcy seepage equation to calculate the matrix pore water pressure

[0102] Then solve the fracture non-Darcy seepage equation and calculate the fracture water pressure

[0103] Among them, n represents the time step calculation number, and k represents the iteration step calculation number.

[0104] Calculating matrix water pressure tolerance and fracture water pressure tolerance

[0105]

[0106] Among them, the left side of formula (4) and (5) is the symbol for matrix water pressure tolerance and fracture water pressure tolerance. The absolute value has no practical meaning and only indicates that the number is a non-negative number. The subscripts 1 and 2 are used to distinguish the two symbols. The double absolute value symbol plus the subscript 2 on the right side of formula (4) and (5) indicates the 2-norm of the column vector, which is actually the 2-norm of the difference of the pore water pressure array calculated in two iteration steps divided by the 2-norm of the pore water pressure array in the previous iteration step.

[0107] if and make Go to 3) to solve the stress equilibrium equation;

[0108] if or And k <k max , let k = k + 1, go to 2), continue to iteratively solve the fracture non-Darcy seepage and matrix seepage equations;

[0109] 3) Solve the stress equilibrium equation and calculate the displacement increment Strain increment and stress increment Update the displacement vector, strain and stress tensor of the l+1 iteration step:

[0110]

[0111] Calculate displacement tolerance

[0112]

[0113] Furthermore, the displacement tolerance is used to judge whether the time step iterative calculation has converged. If it has converged, it will go to the next time step for calculation. Otherwise, it will go to the previous step for iterative calculation. It also includes:

[0114] If the displacement tolerance is less than or equal to the third iteration tolerance, proceed to the next time step calculation;

[0115] If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has not reached the maximum value, the number of iteration steps is incremented by 1, and the process goes back to the previous step to continue the iterative calculation.

[0116] Furthermore, the method of using an adaptive time step to determine whether the calculation of the previous time step needs to be repeated also includes:

[0117] If the number of iterations is less than half of the maximum number of iterations, increase the time step;

[0118] If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has exceeded the maximum value, the time step is reduced and the calculation of the previous time step is repeated.

[0119] In this embodiment, it also includes:

[0120] 4) Convergence judgment:

[0121] The displacement tolerance is used to determine whether the time step iteration calculation has converged. Go to the next time step calculation, t n =t n +Δt, n=n+1; if And l <l max , let l = l + 1, go to the previous step to continue iterative calculation;

[0122] Adaptive time step is used. The specific implementation method is as follows: if k < 0.5 × k max , Δt=Δt×1.1, if And l>l max , let Δt = Δt × 0.5 and repeat the calculation of the previous time step.

[0123] 5) End judgment

[0124] If t n <t max , go to the third step, otherwise, the calculation ends.

[0125] The embodiment of the present invention is subjected to actual simulation experiments, and combined with a typical fractured rock mass seepage stress coupling problem, a fractured rock mass numerical calculation model and method based on an embedded fracture network proposed by the present invention are further described, including the following steps:

[0126] 1) Input information such as fracture and matrix unit mesh, boundary conditions, material parameters and control parameters

[0127] The calculation model of the embodiment adopts a rectangular area with a height of 50m and a width of 11m. There is a crack with a length of 20m in the middle. The pore matrix calculation grid and the crack calculation grid of the model are relatively independent. The total number of matrix calculation grid units is 550, the total number of nodes is 612, the total number of crack grid units is 21, and the total number of nodes is 20. Figure 3 shown.

[0128] The boundary conditions are as follows: 1) The upper boundary is a permeable boundary, the pore pressure is equal to the atmospheric pressure, and the other boundaries are impermeable boundaries; 2) The upper boundary is subjected to a uniform load of 20 MPa during the entire calculation time, and the left, right and bottom boundaries are all fixed support constraints. The fracture permeability is 1.0×10 -10 m 2 The horizontal permeability of the matrix is ​​2.0×10 -15 m 2 , the vertical permeability is 1.0×10 -15 m 2 , the density of water is 1000kg / m 3 , the compression factor is 5.0×10 -9 Pa -1 , the dynamic viscosity coefficient is 1.0×10 -3 Pa·s, the rock mass elastic modulus is 10 GPa, and the Poisson's ratio is 0.25. The calculation adopts a variable step size, with an initial step size of 0.01h, a maximum step size of 10h, and a simulation time of 30d.

[0129] 2) Given the initial values ​​and initial stresses of the basic unknown quantities

[0130] The initial pore water pressure is 0 (atmospheric pressure), the initial fracture water pressure is 0, the initial displacement is 0, the influence of gravity is not considered, and the initial stress is 0.

[0131] 3) Results of non-Darcy seepage stress coupling solution in fracture-matrix system

[0132] Figure 4The distribution of rock matrix water pressure calculated by this model at the end of the calculation time (30 days) is plotted. It can be seen that the permeability of the fracture is good. In the middle of the calculation area, the matrix water pressure gradient is small, and the water in the medium can be discharged through the fractures quickly, indicating that this calculation can better simulate the seepage coupling process of the fracture-matrix system. Figure 5 The horizontal and vertical displacement distributions at the end of the calculation time (30 days) were plotted respectively. It can be seen that the medium deformation is mainly manifested as settlement deformation, the horizontal displacement is small, and it shows obvious discontinuous characteristics, which is related to the deformation of the cracks.

[0133] In the several embodiments provided in the present application, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of units is only a logical function division. There may be other division methods in actual implementation. For example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be an indirect coupling or communication connection through some interfaces, devices or units, which can be electrical, mechanical or other forms.

[0134] In addition, each functional unit in each embodiment of the present application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The above integrated unit can be implemented in the form of hardware or in the form of software functional units. The above is only an implementation method of the present application, and does not limit the patent scope of the present application. Any equivalent structure or equivalent process transformation made using the contents of the specification and drawings of this application, or directly or indirectly used in other related technical fields, is also included in the patent protection scope of the present application.

[0135] The specific implementation methods of the invention are described in detail above, but they are only examples, and the present application is not limited to the specific implementation methods described above. For those skilled in the art, any equivalent modification or substitution of the invention is also within the scope of the present application, and therefore, the equalization, modification, and improvement made without departing from the spirit and principle of the present application should be included in the scope of the present application.

Claims

1. A non-Darcy seepage stress coupling simulation method based on fracture and matrix grid, characterized in that: include: Inputting a simulated matrix grid, establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining the intersection reference information of the mesh units of the simulated fracture grid and the simulated matrix grid; Selecting initial water pressure and initial displacement values ​​of grid nodes according to the simulated fracture grid and the simulated matrix grid; The time-step iterative solution of fracture non-Darcy seepage and matrix seepage is used to obtain the water pressure tolerance and judge the convergence of the time-step iterative calculation; According to the convergence and time step information, the solution of non-Darcy seepage stress coupling is obtained.

2. A non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 1, characterized in that: The input simulation matrix grid, establish the simulation fracture grid according to the attribute information of the fracture, obtain the intersection reference information of the mesh unit of the simulation fracture grid and the simulation matrix grid, and also include: Inputting a simulation matrix grid, and obtaining a matrix grid unit and a matrix grid node according to the simulation matrix grid; Establishing a simulated fracture grid according to the attribute information of the fracture, and obtaining fracture grid units and fracture grid nodes according to the simulated fracture grid; Obtaining intersection reference information of the matrix grid unit and the fracture grid unit; The intersection reference information at least includes: Boundary conditions, material parameters and control parameters.

3. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 1, characterized in that: The selecting of the initial water pressure and initial displacement value of the grid nodes according to the simulated fracture grid and the simulated matrix grid also includes: Selecting the initial fracture water pressure of the fracture grid node according to the simulated fracture grid; The initial matrix water pressure and the initial matrix displacement of the matrix grid nodes are selected according to the simulation matrix grid.

4. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 1, characterized in that: The method adopts the time step iteration to solve the fracture non-Darcy seepage and matrix seepage to obtain the water pressure tolerance and judge the convergence of the time step iteration calculation, and also includes: The seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system are solved in sequence; By calculating the matrix and fracture flow exchange terms, the matrix and fracture seepage coupling solution is realized. By calculating the matrix water pressure and fracture water pressure tolerance, it is determined whether the time step iterative calculation converges.

5. A non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 4, characterized in that: The method of obtaining the solution result of non-Darcy seepage stress coupling based on the convergence and time step information also includes: Use the displacement tolerance to determine whether the time step iterative calculation has converged. If so, proceed to the next time step calculation; otherwise, proceed to the previous step to continue the iterative calculation. Adopt adaptive time step to determine whether the calculation of the previous time step needs to be repeated; When the simulation calculation time is reached, the solution of non-Darcy seepage stress coupling is obtained.

6. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 2, characterized in that: The step of obtaining the intersection reference information of the matrix grid unit and the crack grid unit further includes: Given a simulation calculation time, set the initial time step, first iteration tolerance, second iteration tolerance, third iteration tolerance, and the maximum number of iterations between two iteration steps.

7. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 3, characterized in that: When selecting the initial fracture water pressure of the fracture grid nodes according to the simulated fracture grid and selecting the initial matrix water pressure and initial matrix displacement of the matrix grid nodes according to the simulated matrix grid, it also includes: Let the time step n = 0, select the matrix water pressure of the matrix grid nodes, the fracture water pressure of the fracture grid nodes and the displacement of the matrix grid nodes as the basic unknowns, and give the initial values ​​of the displacement, matrix water pressure, fracture water pressure and stress time step as the initial matrix displacement, initial matrix water pressure, initial fracture water pressure and initial stress time step.

8. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 4, characterized in that: When solving the seepage equation of the matrix grid system and the non-Darcy seepage equation of the fracture grid system in sequence, it also includes: Solve the fracture non-Darcy seepage equation, matrix Darcy seepage equation and matrix stress balance equation in sequence; Set the iteration step to 0, and give the iteration step displacement, initial iteration matrix water pressure, initial iteration fracture water pressure and initial stress iteration value; The method further includes: calculating the matrix and fracture flow exchange term to achieve matrix and fracture seepage coupling solution, and calculating the matrix water pressure and fracture water pressure tolerance to determine whether the time step iterative calculation converges. Solve the matrix Darcy flow equation and calculate the pore water pressure; Solve the non-Darcy flow equation in fractures and calculate the fracture water pressure; Calculate matrix water pressure tolerance and fracture water pressure tolerance; According to the relationship between the tolerance of the matrix water pressure and the fracture water pressure and the tolerance of the first iteration and the tolerance of the second iteration, it is selected to solve the stress balance equation or to continue to iteratively solve the fracture non-Darcy seepage and matrix seepage equations; When solving the stress equilibrium equation, the displacement increment, strain increment and stress increment are calculated, the displacement vector, strain and stress tensor of the iteration step are updated, and the displacement tolerance is calculated.

9. A non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 8, characterized in that: The displacement tolerance is used to determine whether the time step iterative calculation has converged. If it has converged, it will proceed to the next time step calculation. Otherwise, it will proceed to the previous step to continue the iterative calculation. It also includes: If the displacement tolerance is less than or equal to the third iteration tolerance, proceed to the next time step calculation; If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has not reached the maximum value, the number of iteration steps is incremented by 1, and the process goes back to the previous step to continue the iterative calculation.

10. The non-Darcy seepage stress coupling simulation method based on fracture and matrix grid according to claim 8, characterized in that: The method of using an adaptive time step to determine whether the calculation of the previous time step needs to be repeated also includes: If the number of iterations is less than half of the maximum number of iterations, increase the time step; If the displacement tolerance is greater than the third iteration tolerance and the number of iteration steps has exceeded the maximum value, the time step is reduced and the calculation of the previous time step is repeated.

Images