A numerical simulation method for karst reservoir evolution

Abstract

The present invention discloses a method for numerical simulation of karst evolution, comprising the following steps: step 1, establishing a fracture and matrix coupling geological model; step 2, establishing a multi-type reservoir fluid-solid-chemical coupling mathematical model; step 3, forming a fast and stable fully coupled numerical solution algorithm. By using an embedded discrete fracture method and combining a fracture mechanics constitutive model to process fracture or fault diversion and deformation processes, not only can the calculation efficiency be improved, but also the control effect of fracture dissolution and mechanical deformation on the karst process can be accurately reflected. The present invention can simulate the evolution law of reservoir porosity and permeability under complex karst environments, and provide a theoretical basis for revealing the formation mechanism of carbonate fracture-cavity reservoirs.

Description

Technical Field

[0001] The invention relates to the field of oil and gas field development geology, and in particular to a karst reservoir evolution numerical simulation method. Background Art

[0002] The development plan of carbonate fracture-cavity reservoirs depends largely on the spatial size, spatial morphology, connectivity and combination relationship of fractures and caves in the reservoir. Therefore, clarifying the fracture-cavity combination morphology of carbonate reservoirs and understanding the laws of karst evolution are one of the key technologies for achieving efficient development of carbonate fracture-cavity reservoirs. Inverting the fracture-cavity combination characteristics and laws of carbonate reservoirs through evolutionary simulation of the karst reservoir formation process is an important means and way to describe the reservoir.

[0003] At present, there are two main types of numerical simulation models related to the evolution of carbonate karst reservoirs: one is the pipe flow model for the cave formation process, and the other is the dual-scale wormhole model for the acid fracturing process.

[0004] In 1999, Kaufmann G. and Braun J. published a paper titled Karst Aquifer Evolution in Fractured Rocks in the journal Water Resources Research, establishing a pipe flow model for the cave formation process. In 2009, Kaufmann G. published another paper titled Modelling Karst Geomorphology on Different Time Scales, which improved the model. This model equates karst caves of different shapes and sizes to interconnected cylindrical conduits with different radii, and uses the Hagen-Poiseuille law to calculate the fluid flow in the conduit. By solving the solute transport-diffusion equation with chemical reactions, the fluid flow velocity, solute concentration and chemical reaction rate are obtained. Then, the relationship between the conduit radius and the fluid flow velocity and reaction rate is established, and the shape and size changes and extension laws of the conduit under the current flow velocity and reaction rate are determined, thereby simulating the evolution of the cave in the reservoir under karstification. This method does not consider the flow in the matrix, nor the controlling role of fractures in the entire karst process. In addition, the influence of geostress on the matrix and fractures, as well as the flow of fluid in the cave, is ignored.

[0005] In 2005, Panga M. et al. published a paper titled Two-scale continuum model for simulation of wormholes in carbonate acidification in AIChE Journal, establishing a continuous medium model for carbonate acidification. In 2007, Kalia N. and Balakotaiah V. published a paper titled Modeling and Analysis of Wormhole Formation in Reactive Dissolution of Carbonate Rocks in Chemical Engineering Science, establishing a model for the formation of wormholes in carbonate acidification. This type of model simulates acidification wormholes and studies their improvement of permeability in the near-wellbore zone to achieve the purpose of increasing the production of a single well. First, a method that uniformly treats the matrix, dissolved pores, fractures and caves as a single matrix with a certain porosity and permeability, and obtains the fluid flow velocity, solute concentration and chemical reaction rate by solving the solute transport-diffusion equation with chemical reactions. Then, the mathematical dependence of porosity and permeability on flow rate, solute concentration and chemical reaction rate is established to simulate the evolution of porosity and permeability of the reservoir under the action of high-pressure and highly corrosive fracturing fluid. At present, this method simulates the acidification process without considering the influence of stress on the matrix, fracture porosity and permeability, and also fails to consider the control of fractures on the migration speed of the chemical reaction front. The time and spatial scales of the acidification process are small, and only the near-well area within a hundred meters is considered. The process of reservoir porosity and permeability changes measured in hours cannot be used for the evolution of pores, fractures and caves in reservoirs within a few kilometers or even larger, and over a hundred years or even longer.

[0006] The formation process of dissolved fractures, dissolved pores, dissolved caves and their composite reservoirs in carbonate reservoirs is extremely complex. At present, there is no mathematical model that considers the karst evolution process of fracture-cavity reservoirs at home and abroad. Summary of the invention

[0007] In order to simulate the evolution and formation process of karst pore reservoirs in complex geological environments, realize the dynamic characterization of the formation process of fracture-cavity reservoirs, and reveal the formation mechanism of fracture-cavity reservoirs in carbonate rocks, this paper proposes a karst evolution numerical simulation method based on fluid-solid-chemical full coupling of embedded discrete fractures, which comprehensively considers fluid flow, solute migration and diffusion, chemical reaction and rock deformation process, and realizes the dynamic characterization of the formation process of fracture-cavity reservoirs.

[0008] The present invention discloses a karst reservoir evolution numerical simulation method, comprising the following steps:

[0009] Step 1: Establish a fracture and matrix coupling geological model;

[0010] Step 2: Establish a multi-type reservoir fluid-solid-chemical coupling mathematical model;

[0011] Step 3: Form a fast and stable fully coupled numerical solution algorithm.

[0012] In one embodiment, step 1 includes the following sub-steps:

[0013] Step 1.1, mesh the matrix and fractures according to the geometric boundary data of the simulation area and the vertex data of each natural fracture;

[0014] Step 1.2, establishing geometric connection relationships between matrix and cracks, and between cracks;

[0015] Step 1.3, calculating the length of the intersection line between the cracks, the intersection area between the cracks and the matrix, and the corresponding equivalent distance, and calculating the conductivity coefficient in combination with the permeability of the cracks and the matrix;

[0016] Step 1.4: Establish a connection table, including three connection relationships between matrix and matrix, matrix and fracture, and fracture and fracture, to obtain a fracture and matrix coupling geological model.

[0017] In one embodiment, in step 1.1, a Cartesian orthogonal grid or a corner point grid is used for grid generation.

[0018] In one embodiment, step 2 includes the following sub-steps:

[0019] Step 2.1, based on the main physical and chemical processes of karst evolution and their coupling relationships, derive the control equations of the fluid-solid-chemical coupling system;

[0020] Step 2.2: Based on the control equation in step 2.1 and in combination with chemical reaction kinetic equations of different types of reservoirs, a multi-type reservoir fluid-solid-chemical coupling mathematical model is established.

[0021] In one embodiment, the control equations of the fluid-solid-chemical coupling system include: the fluid mass conservation equation shown in formula (1), the solute convection-diffusion equation with chemical reaction shown in formula (2), and the mechanical equilibrium equation shown in formula (3).

[0022]

[0023]

[0024]

[0025] Among them, the porosity in formula (2) is obtained by evolving formula (4),

[0026]

[0027] Where: p is the fluid pressure, ε is the strain tensor, D is the fourth-order elastic matrix tensor, δ is the Kronecker tensor, ρ f and ρ b are the fluid density and the average volume density of the fluid in the rock skeleton and pores, is the Darcy velocity, C f and C s are the calcium ion concentration in the pores and the calcium ion concentration on the solid-liquid interface, respectively. e is the diffusion coefficient tensor, k c is the local mass exchange coefficient, a v is the surface area per unit volume for chemical reactions, ε V is the body strain, g is the gravitational acceleration, b is the Biot coefficient, K s is the bulk modulus of matrix particles, μ is the viscosity, k is the permeability, B is the fluid volume coefficient, φ and φ0 are the current porosity and initial porosity, t is the simulation time, q f is the mass source / sink term.

[0028] In one embodiment, step 3 includes the following sub-steps:

[0029] Step 3.1, numerically discretizing the multi-type reservoir fluid-solid-chemical coupling mathematical model described in step 2;

[0030] Step 3.2, constructing a Jacobian matrix structure based on the connection table described in step 1.4 to form a nonlinear / linear matrix solution format;

[0031] Step 3.3: Use the finite volume method and finite element method to discretely solve the fluid-solid-chemical coupling mathematical model of multiple types of reservoirs.

[0032] In one embodiment, the finite volume method is used to solve for fluid pressure and solute concentration, and the finite element method is used to solve for stress, strain, and displacement.

[0033] In one embodiment, the numerical values ​​in step 3.1 are discretized to obtain formula (19):

[0034]

[0035] Where: R u , R p and R Care the residuals of the mechanical equilibrium equation, the fluid mass conservation equation and the solute convection-diffusion equation with chemical reaction, K is the stiffness matrix, Δt is the time step, L and M are the poroelastic coupling matrices, The compression matrix related to fluid pressure, is the compression matrix related to fluid stress, F is the conductivity matrix related to fluid pressure, Q is the conductivity matrix related to fluid concentration, superscripts u, p, C, n and k represent displacement, fluid pressure, concentration, time step and iteration step respectively, subscripts i / j represent units, subscripts a / b represent nodes, Δu, Δp and ΔC represent the changes of displacement, fluid pressure and calcium ion concentration in Newton iteration step respectively, C n represents the calcium ion concentration at the initial moment of time step n.

[0036] In one embodiment, the numerical discretization equation (19) in step 3.3 is solved using the Newton-Raphson method to obtain equation (29):

[0037]

[0038]

[0039] Where: X is the variable system to be solved, δ is the variable increment symbol, R is the residual within a certain Newton iteration step, P and S represent the main variable and secondary variable respectively, n and k represent the time step and iteration step respectively.

[0040] In one embodiment, the solution using the Newton-Raphson method comprises the following steps:

[0041] Firstly, the Newton-Raphson method is used to uniformly solve equations (1), (2) and (3) to obtain fluid pressure, displacement, stress, strain and solute concentration;

[0042] Then, the porosity is calculated according to formula (4) and substituted into the Carman-Kozeny equation to solve for the pore throat radius, permeability, and surface area per unit volume for chemical reaction;

[0043] Finally, the diffusion coefficient, flow velocity and density are obtained through the relationship between pore structure and physicochemical parameters.

[0044] The above-mentioned technical features can be combined in various technically feasible ways to produce new embodiments, as long as the purpose of the present invention can be achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The present invention will be described in more detail below based on embodiments and with reference to the accompanying drawings, wherein:

[0046] Figure 1 It is a schematic flow chart of a numerical simulation method for fluid-solid-chemical coupling karst evolution based on embedded discrete fractures in an embodiment of the present invention;

[0047] Figure 2 is a schematic diagram of a fluid-solid-chemical coupling system in an embodiment of the present invention;

[0048] Figure 3 is an interpolation function of the finite volume finite element method in the embodiment of the present invention;

[0049] Figure 4 It is a submodule project for realizing the karst evolution model in the embodiment of the present invention;

[0050] Figure 5 is a karst evolution structural fracture model in an embodiment of the present invention;

[0051] Figure 6 The karst evolution structural fracture model in the embodiment of the present invention simulates the distribution and change of reservoir porosity after 100 years;

[0052] Figure 7 The karst evolution structural fracture model in the embodiment of the present invention simulates the distribution and change of reservoir porosity after 200 years;

[0053] Figure 8 The karst evolution structural fracture model in the embodiment of the present invention simulates the distribution and change of reservoir porosity after 300 years.

[0054] In the drawings, like reference numerals are used for like parts. The drawings are not necessarily to scale. DETAILED DESCRIPTION

[0055] The present invention will be further described below in conjunction with the accompanying drawings.

[0056] like Figure 1 The formula is a schematic flow chart of the numerical simulation method for karst evolution in an example of the present invention. The numerical simulation method for karst evolution in the present invention is specifically introduced below.

[0057] The present invention proposes a karst evolution numerical simulation method based on fluid-solid-chemical full coupling of embedded discrete fractures, which mainly includes three steps: step 1, establishing a fracture and matrix coupling geological model; step 2, establishing a multi-type reservoir fluid-solid-chemical coupling mathematical model; step 3, forming a fast and stable full-coupled numerical solution algorithm.

[0058] The three steps are described in detail below.

[0059] Step 1: Establish a fracture-matrix coupling geological model.

[0060] The fracture system and the matrix system are coupled as two independent systems, and the fracture system is processed by explicit discretization instead of the equivalent processing used by the traditional dual medium method. The advantage of this method is that it does not require equivalent hidden processing of the fractures, but keeps the original geometric morphology of the natural fractures for display processing, and has a computational accuracy very close to that of the discrete fracture method (DFN) after local encryption processing. There is no need to encrypt the area where the fracture is located, and the effect of the high conductivity channel of the fracture can be simulated in a simple structured grid, which has high computational efficiency. Therefore, the advantages of establishing a fracture-matrix coupling geological model based on the embedded discrete fracture method (EDFM) can be summarized as follows: 1) The fracture modeling is convenient, and the study area can be divided into a simple structured grid, and the four vertex coordinates of each fracture are given. The pre-processing program can automatically generate a connection table for subsequent calculations; 2) There is no need to encrypt the grid in the area near the fracture, and the computational efficiency is high; 3) The computational accuracy is high and the chemical reaction front can be accurately simulated.

[0061] The specific implementation steps of the discrete fracture modeling method, that is, the sub-steps included in step 1 are as follows:

[0062] Step 1.1, meshing the matrix and the fractures according to the geometric boundary data of the simulation area and the vertex data of each natural fracture (the vertex may not be within the simulation area, and the non-planar natural fracture may be spliced ​​by multiple planar fractures);

[0063] Step 1.2, establishing the geometric connection relationship between the matrix and the cracks, and between the cracks (i.e. establishing a non-adjacent connection table);

[0064] Step 1.3, calculating the length of the intersection line between the cracks, the intersection area between the cracks and the matrix, and the corresponding equivalent distance, and calculating the conductivity coefficient in combination with the permeability of the cracks and the matrix;

[0065] Step 1.4: Establish a connection table, including three connection relationships between matrix and matrix, matrix and fracture, and fracture and fracture, to obtain a fracture and matrix coupling geological model.

[0066] In a preferred embodiment, in step 1.1, a Cartesian orthogonal grid or a corner point grid is used for grid generation.

[0067] Step 2: Establish a multi-type reservoir fluid-solid-chemical coupling mathematical model.

[0068] The establishment of a fracture and matrix coupling geological model is the basis of the fluid-solid-chemical coupling simulation method. Next, the control equations of the fluid-solid-chemical coupling system (i.e., the mathematical model behind the geological model) need to be derived. First, based on the main physical and chemical processes of karst evolution and their coupling relationships, the control equations of the fluid-solid-chemical coupling system are derived; then, combined with the chemical reaction kinetic equations of different types of reservoirs, a multi-type reservoir fluid-solid-chemical coupling mathematical model is established.

[0069] Karst evolution involves chemical reactions, rock deformation, solute migration, fluid convection and diffusion, etc. Each process can be generalized into a corresponding physical mathematical model based on certain assumptions and simplifications. The basic principle of model generalization is: reasonable assumptions, appropriate simplifications, and accurate reflection of the real physical process. Figure 2 The figure shows the principle diagram of the fluid-solid-chemical coupling system. The chemical reaction of the salt rock mass changes the pore structure and specific surface area of ​​the rock, causing corresponding changes in the porosity, permeability, and reaction rate, and ultimately having a significant impact on fluid flow and solute transport; fluid flow is accompanied by solute transport; solute transport changes the local reactant concentration and affects the reaction rate. Figure 2 The figure shows the various processes involved in karst evolution and the coupling relationship between them: (1) Fresh water and groundwater containing CO2 react chemically with carbonate rocks to form different types of dissolution structures. This process depends on the mineral composition of the salt rock, chemical reaction rate, flow type, reservoir heterogeneity, and the rate at which fluids filter into the matrix. (2) Driven by pressure potential, fluids seep or flow freely in the matrix, fractures, and dissolution caves. (3) Solute (Ca 2+ ,Mg 2+ , H + ,CO3 2- (4) Under conditions such as tectonic stress and self-weight stress, the rock mass deforms, including changes in pore volume and crack opening. The rock mass deforms under stress, which changes the porosity and permeability of the matrix and the crack opening, thereby affecting the flow of fluid and solute transport. The change in crack opening is mainly related to the magnitude of stress, mineral dissolution, transportation, and precipitation on the fracture surface. Among them, chemical dissolution changes the pore structure of the rock mass, affects the specific surface area, permeability, and porosity of the pores, thereby changing the flow rate of the fluid and the distribution of solutes; solute transport changes the solute concentration, thereby affecting the chemical dissolution process; fluid flow causes changes in concentration, and concentration changes in turn affect the chemical dissolution process. Therefore, the entire system jointly affects the karst evolution process through coupling.

[0070] The control equations of the system are specifically introduced below. The control equations of the fluid-solid-chemical coupling system include: the fluid mass conservation equation shown in formula (1), the solute convection-diffusion equation with chemical reaction shown in formula (2), and the mechanical equilibrium equation shown in formula (3):

[0071]

[0072]

[0073]

[0074] Among them, the porosity in formula (2) is evolved from formula (4), which represents the porosity change caused by dissolution and changes in pressure and stress;

[0075]

[0076] In the above formulas: p is the fluid pressure, ε is the strain tensor, D is the fourth-order elastic matrix tensor, δ is the Kronecker tensor, ρ f and ρ b are the fluid density and the average volume density of the fluid in the rock skeleton and pores, is the Darcy velocity, C f and C s are the calcium ion concentration in the pores and the calcium ion concentration on the solid-liquid interface, respectively. e is the diffusion coefficient tensor, k c is the local mass exchange coefficient, a v is the surface area per unit volume for chemical reactions, ε V is the body strain, g is the gravitational acceleration, b is the Biot coefficient, K s is the bulk modulus of matrix particles, μ is the viscosity, k is the permeability, B is the fluid volume coefficient, φ and φ0 are the current porosity and initial porosity, t is the simulation time, q f is the mass source / sink term.

[0077] In addition to the above control equations, there are also some other auxiliary equations, such as strain-displacement equation, stress-strain equation, crack constitutive equation, Darcy seepage equation, etc.

[0078] The mass transport of solutes in the pore network is determined by the following equation:

[0079]

[0080]

[0081]

[0082] In the formula, Sh is the Sherwood number or dimensionless mass transfer coefficient, k c is the local mass exchange coefficient, r P is the average pore radius, D m is the effective molecular expansion coefficient of the acid at any shear rate, is the asymptotic Sherwood number, m p is the ratio of pore length to radius, is the Reynolds number at the pore scale, is the ratio of inertial force to viscous force, υ is the dynamic viscosity, is the Schmidt number, α os is a constant coefficient that depends on pore connectivity, D eX is the longitudinal (X direction) diffusion coefficient, D eT is the lateral (Y or Z direction) diffusion coefficient, λ x and λ T is a constant coefficient that depends on the pore structure, and φ is the current porosity.

[0083] The physicochemical parameters of rocks in the Darcy model (local pore size, pore surface area and permeability) depend on the pore structure. The changes in physicochemical parameters can be considered by establishing the relationship between pore structure and physicochemical parameters (Carman-Kozeny equation), as shown in Equations (8) and (9):

[0084]

[0085]

[0086] Where: I is the second-order unit tensor, k0 is the initial average permeability, φ0 is the initial porosity, r0 is the initial average pore radius, a0 is the initial specific surface area, β is the pore expansion coefficient, k, φ, r P 、a v They are the current permeability, current porosity, average pore radius, and surface area per unit volume available for chemical reactions.

[0087] Step 3: Form a fast and stable fully coupled numerical solution algorithm.

[0088] After obtaining the multi-type reservoir fluid-solid-chemical coupling mathematical model in step 2, it is necessary to solve the mathematical model. Step 3 may include the following sub-steps:

[0089] Step 3.1, numerically discretizing the multi-type reservoir fluid-solid-chemical coupling mathematical model described in step 2;

[0090] Step 3.2, constructing a Jacobian matrix structure based on the connection table described in step 1.4 to form a nonlinear / linear matrix solution format;

[0091] Step 3.3: Use finite volume and finite element methods to discretely solve the fluid-solid-chemical coupling mathematical model of multiple types of reservoirs.

[0092] Preferably, the finite volume method is used to solve the fluid pressure and solute concentration, and the finite element method is used to solve the stress, strain and displacement.

[0093] like Figure 3 As shown in the figure, it is the interpolation function of the finite volume finite element method. It includes the interpolation functions of pore pressure, concentration and displacement, where the displacement variables are located at the block nodes and the pore pressure and concentration are located at the center of the block. Figure 3 In the figure, (a) is the numerical distribution of displacement interpolation function within the unit, (b) is the numerical distribution of pore pressure and concentration interpolation function within the unit; (c) is a schematic diagram of the corresponding position of the main variables within the unit.

[0094] The specific expression of the interpolation function is shown in formula (10):

[0095]

[0096] Among them, the displacement interpolation function ψ i The value is 1 at node i and 0 at all other nodes. It is obtained by interpolation of the hexahedral unit finite element interpolation function shown in formula (11) within the unit:

[0097]

[0098]

[0099]

[0100] Interpolation functions for pore pressure and concentration and η l The interpolation function satisfies the value of 1 at units a and l and 0 at other units. The interpolation function is 1 or 0 everywhere in the unit, which is different from the displacement interpolation function. Figure 3 Substituting the displacement interpolation function into the control equation shown in formula (3) and arranging it into the residual iteration format, the mechanical equilibrium equation shown in formula (3) can be arranged as shown in formula (14),

[0101]

[0102] in,

[0103]

[0104]

[0105] Substitute the interpolation functions of pore pressure and concentration into the control equations shown in equations (1) and (2) respectively, and organize them into residual iteration format. Then, the fluid mass conservation equation shown in equation (1) can be organized as shown in equation (17):

[0106]

[0107] The solute convection-diffusion equation with chemical reaction shown in formula (2) can be organized as shown in formula (18):

[0108]

[0109] Further arranging equations (14), (17) and (18) yields the Newton iteration matrix form of the entire system as shown in equation (19):

[0110]

[0111] Where: R u , R p and R C are the residuals of the mechanical equilibrium equation, the fluid mass conservation equation and the solute convection-diffusion equation with chemical reaction, K is the stiffness matrix, Δt is the time step, L and M are the poroelastic coupling matrices, The compression matrix related to fluid pressure, is the compression matrix related to fluid stress, F is the conductivity matrix related to fluid pressure, Q is the conductivity matrix related to fluid concentration, superscripts u, p, C, n and k represent displacement, fluid pressure, concentration, time step and iteration step respectively, subscripts i / j represent units, subscripts a / b represent nodes, Δu, Δp and ΔC represent the changes of displacement, fluid pressure and calcium ion concentration in Newton iteration step respectively, C n represents the calcium ion concentration at the initial moment of time step n.

[0112] The matrix blocks in equation (19) are shown in equations (20) to (28):

[0113]

[0114]

[0115]

[0116]

[0117]

[0118]

[0119]

[0120]

[0121]

[0122] For the multi-field coupled karst process, the present invention intends to adopt a fully coupled solution strategy to solve it (such as Figure 1 As shown in the figure, that is, in the same time step, the main variables pressure, displacement and concentration are solved at the same time, and other variables (including porosity, permeability, diffusion number, density and flow rate) can be expressed as functions of the main variables. The value of the secondary variable can be obtained through the functional relationship, and updated to enter the calculation of the next time step. The whole system uses the Newton-Raphson method to linearize and solve the partial differential equation:

[0123]

[0124]

[0125] Where X is the variable system to be solved, δ is the variable increment symbol, R is the residual within a certain Newton iteration step, P and S represent the main variable and secondary variable respectively, and n and k represent the time step and iteration step respectively.

[0126] The solution system uses fluid pressure, displacement and solute concentration as the main variables, and porosity, permeability, specific surface area, diffusion coefficient, flow rate, density, etc. as secondary variables.

[0127] Preferably, the Newton-Raphson method comprises the following steps:

[0128] Firstly, the Newton-Raphson method is used to uniformly solve equations (1), (2) and (3) to obtain fluid pressure, displacement, stress, strain and solute concentration;

[0129] Then, the porosity is calculated according to formula (4) and substituted into the Carman-Kozeny equation to solve for the pore throat radius, permeability, and surface area per unit volume for chemical reaction;

[0130] Finally, the diffusion coefficient, flow velocity and density are obtained through the relationship between pore structure and physicochemical parameters.

[0131] According to the characteristics of the mathematical model of karst evolution, a corresponding C++ programming framework is designed. According to the modular design concept, each functional module is independent of each other and follows the characteristics of code sharing, achieving the design goals of high code reuse rate and strong module scalability. Figure 4As shown in the figure, the entire karst evolution simulation system KarstEvolSys is divided into three parts according to its functions: data processing part, simulator driving part and model coupling module part. Among them, the data processing part consists of three sub-modules, including the FileIO module for data file input and output management, the Database module responsible for data storage and query functions, and the 3D_Visualization module responsible for outputting one-dimensional curves, two-dimensional surfaces, and three-dimensional isosurfaces; the simulator driving module includes the linear solver module LinearSolver (this model uses the Intel math core function library MKL), the constitutive model module ConstituModel for processing elastic deformation and elastic-plastic deformation, the module TP-MPFA module that automatically generates a two-point or multi-point flow format connection table based on the grid, and the module ChemProps for calculating chemical properties; the third part is the model coupling module, including the fluid flow module FluidFlow, the rock deformation module GeoMech and the chemical reaction module ChemReac. The three modules are combined into an overall solution system through the coupling matrix, so that the main variables fluid pressure p, displacement u and solute concentration C can be calculated at the same time step. f , and at each iteration step, the secondary variables permeability k, porosity φ, and surface area a per unit volume for chemical reactions are V , diffusion coefficient S h Wait for the update to complete.

[0132] The present invention proposes a fluid-solid-chemical fully coupled karst reservoir evolution numerical simulation method, which uses an embedded discrete fracture model to simulate the control effect of fractures on karst evolution, while considering the influence of stress on the permeability and porosity changes of the rock mass. The mixed finite volume finite element method (mixed FV / FE) is used to discretize the entire system, which can simulate the evolution of fractures and caves in complex environments and provide a theoretical basis for carbonate fracture-cavity reservoir modeling. The method can be implemented by software programming. The software system flow can be described as follows:

[0133] (1) Project settings: set the files and parameters required for simulating the karst evolution process;

[0134] (2) Calculation example settings: setting initial and boundary conditions such as flow rate, fluid pressure, concentration, and stress;

[0135] (3) Establish a coupled geological model of fractures and matrix based on the geometric parameters of structural fractures;

[0136] (4) Establish mathematical models for karst evolution of multiple types of reservoirs;

[0137] (5) Initialize fluid flow, solute concentration, and in situ geostress according to the fluid-solid-chemical coupling method;

[0138] (6) Using a fully coupled solution algorithm to solve the fluid-solid-chemical coupled karst evolution system;

[0139] (7) Output required parameters such as porosity, permeability, etc. as needed;

[0140] (8) Use Tecplot drawing software to process and display the simulation results.

[0141] The method proposed in the present invention is used to simulate the evolution process of a reservoir containing natural structural fractures. The main simulation process is: ① Establish a digital model grid model and coarsen the grid data of the three-dimensional geological model; ② Use the embedded discrete fracture method introduced above to perform preprocessing according to the natural structural fractures, generate a non-adjacent connection table, and construct a Jacobian matrix structure; ③ Prepare an input text file and initialize the fluid-solid-chemical coupling model; ④ Run the simulator to obtain the simulation results and process the results. Figure 5 As shown in the figure, the entire simulation area is 500m×500m×10m (10m in the horizontal direction perpendicular to the paper), the matrix porosity is 5%, the permeability is 5mD, the simulation area contains a set of conjugate natural structural fractures, the fracture permeability is 500mD, the annual rainfall is 20mm, the surface temperature is 25℃, the surface atmospheric pressure is 1 standard atmosphere (i.e. 0.1MPa), the fluid pressure and temperature distribution are calculated according to the hydrostatic pressure and geothermal gradient (according to 25℃ / km), respectively, the carbon dioxide partial pressure is 0.05 standard atmosphere (i.e. 0.005MPa), and the carbon dioxide equilibrium concentration is 2.1mol / m 3 Then, the carbon dioxide concentration distribution was calculated according to the carbon dioxide concentration balance theory, and the simulation finally obtained 100 years ( Figure 6 )、200 years( Figure 7 ) and 300 years ( Figure 8 ). The simulation results show the distribution and changes of reservoir porosity at different times. From the results, it can be seen that the present invention can simulate the karst evolution process of natural fractured reservoirs over time, and has a good display effect, which can provide guidance for the formulation of efficient development plans for carbonate fracture-cavity reservoirs.

[0142] Although the present invention has been described with reference to preferred embodiments, various modifications may be made thereto and parts thereof may be replaced by equivalents without departing from the scope of the present invention. In particular, the various technical features mentioned in the various embodiments may be combined in any manner as long as there are no structural conflicts. The present invention is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A karst reservoir evolution numerical simulation method, comprising the following steps: Step 1: Establish a fracture and matrix coupling geological model; Step 1 includes the following sub-steps: Step 1.1, mesh the matrix and fractures according to the geometric boundary data of the simulation area and the vertex data of each natural fracture; Step 1.2, establishing geometric connection relationships between matrix and cracks, and between cracks; Step 1.3, calculating the length of the intersection line between the cracks, the intersection area between the cracks and the matrix, and the corresponding equivalent distance, and calculating the conductivity coefficient in combination with the permeability of the cracks and the matrix; Step 1.4, establish a connection table, including three connection relationships between matrix and matrix, matrix and fracture, and fracture and fracture, and obtain a fracture and matrix coupling geological model; Step 2: Establish a multi-type reservoir fluid-solid-chemical coupling mathematical model; Step 2 includes the following sub-steps: Step 2.1, based on the main physical and chemical processes of karst evolution and their coupling relationships, derive the control equations of the fluid-solid-chemical coupling system; Step 2.2, according to the control equation in step 2.1, combined with the chemical reaction kinetic equations of different types of reservoirs, a multi-type reservoir fluid-solid-chemical coupling mathematical model is established; Step 3, forming a fast and stable fully coupled numerical solution algorithm; Step 3 includes the following sub-steps: Step 3.1, numerically discretizing the multi-type reservoir fluid-solid-chemical coupling mathematical model described in step 2; Step 3.2, constructing a Jacobian matrix structure based on the connection table described in step 1.4 to form a nonlinear / linear matrix solution format; Step 3.3: Use finite volume and finite element methods to discretely solve the fluid-solid-chemical coupling mathematical model of multiple types of reservoirs.

2. The karst reservoir evolution numerical simulation method according to claim 1, characterized in that: In the step 1.1, a Cartesian orthogonal grid or a corner point grid is used for grid generation.

3. The karst reservoir evolution numerical simulation method according to claim 1, characterized in that: The control equations of the fluid-solid-chemical coupling system include: the fluid mass conservation equation shown in formula (1), the solute convection-diffusion equation with chemical reaction shown in formula (2), and the mechanical equilibrium equation shown in formula (3). Among them, the porosity in formula (2) is obtained by evolving formula (4), Where: p is the fluid pressure, ε is the strain tensor, D is the fourth-order elastic matrix tensor, v is the Kronecker tensor, ρ f and ρ b are the fluid density and the average volume density of the fluid in the rock skeleton and pores, is the Darcy velocity, C f and C s are the calcium ion concentration in the pores and the calcium ion concentration on the solid-liquid interface, respectively. e is the diffusion coefficient tensor, k c is the local mass exchange coefficient, a v is the surface area per unit volume for chemical reactions, ε V is the body strain, g is the gravitational acceleration, b is the Biot coefficient, K s is the bulk modulus of matrix particles, μ is the viscosity, k is the permeability, B is the fluid volume coefficient, φ and φ0 are the current porosity and initial porosity, t is the simulation time, q f is the mass source / sink term.

4. The method for numerical simulation of karst reservoir evolution according to claim 3, characterized in that: The finite volume method is used to solve the fluid pressure and solute concentration, and the finite element method is used to solve the stress, strain and displacement.

5. The karst reservoir evolution numerical simulation method according to claim 3, characterized in that: The numerical discretization in step 3.1 yields equation (19): Where: R u , R p and R C are the residuals of the mechanical equilibrium equation, the fluid mass conservation equation and the solute convection-diffusion equation with chemical reaction, K is the stiffness matrix, Δt is the time step, L and M are the poroelastic coupling matrices, The compression matrix related to fluid pressure, is the compression matrix related to fluid stress, F is the conductivity matrix related to fluid pressure, Q is the conductivity matrix related to fluid concentration, superscripts u, p, C, n and k represent displacement, fluid pressure, concentration, time step and iteration step respectively, subscripts i / j represent units, subscripts a / b represent nodes, Δu, Δp and ΔC represent the changes of displacement, fluid pressure and calcium ion concentration in Newton iteration step respectively, C n represents the calcium ion concentration at the initial moment of time step n.

6. The karst reservoir evolution numerical simulation method according to claim 5, characterized in that: In step 3.3, the numerical discretization equation (19) is solved using the Newton-Raphson method to obtain equation (29): Where: X is the variable system to be solved, v is the variable increment symbol, R is the residual within a certain Newton iteration step, P and S represent the main variable and secondary variable respectively, n and k represent the time step and iteration step respectively.

7. The method for numerical simulation of karst reservoir evolution according to claim 6, characterized in that: The Newton-Raphson method is used to solve The following steps are involved: Firstly, the Newton-Raphson method is used to uniformly solve equations (1), (2) and (3) to obtain fluid pressure, displacement, stress, strain and solute concentration; Then, the porosity is calculated according to formula (4) and substituted into the Carman-Kozeny equation to solve for the pore throat radius, permeability, and surface area per unit volume for chemical reaction; Finally, the diffusion coefficient, flow velocity and density are obtained through the relationship between pore structure and physicochemical parameters.

Images