Discrete fracture model-based multiphase multicomponent seepage numerical simulation method

Abstract

The invention provides a multiphase multicomponent seepage numerical simulation method based on a discrete fracture model. The method comprises the following steps: step 1, establishing a multi-phase multi-component seepage mathematical model, and solving; step 2, mesh generation of the discrete fracture model; step 3, discretization of a mass conservation equation; step 4, calculating the conductivity between the matrix and the crack; and 5, establishing a nonlinear equation set, and completing solution by adopting an automatic differential technology and a Newton iteration method. According to the method, a complex fracture network is represented through a discrete fracture model, reservoir stress sensitivity and nonlinear seepage characteristics are considered, adsorption and diffusion effects of a surfactant and salt and mechanisms such as capillary pressure, relative permeability curve and osmotic pressure change under the action of the surfactant and the salt are coupled, and a multi-phase multi-component seepage mathematical model is established. According to the method, fracture explicit representation is achieved through a self-adaptive mesh generation method, a finite volume method is adopted for solving, the model is established, and therefore ultra-low permeability-tight oil reservoir development simulation under the influence of a complex fracture network is achieved.

Description

Technical Field

[0001] The invention belongs to the technical field of oil extraction, and in particular relates to a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model. Background Art

[0002] Compared with ordinary low-permeability reservoirs, ultra-low permeability-tight reservoirs have worse reservoir properties, with pore sizes mainly at the submicron and nanometer levels, weak connectivity, low proportion of movable fluid, and low water drive efficiency. At present, ultra-low permeability-tight reservoirs at home and abroad are mostly developed by horizontal wells, and high initial production capacity is obtained through large-scale multi-stage fracturing transformation, and water injection displacement or throughput are adopted to further improve crude oil recovery. In the development of ultra-low permeability-tight reservoirs, not only a complex fracture network system is involved, but also a multi-phase and multi-component seepage coupling process between oil / water (chemical agent) / rock. Therefore, accurate, comprehensive and efficient characterization of the complex fracture network system and multi-phase and multi-component seepage process of ultra-low permeability-tight reservoirs has become an important research content of numerical simulation of such reservoirs.

[0003] Currently, there is no complete numerical simulation method to characterize the complex fracture network system and multi-phase and multi-component seepage process in ultra-low permeability and tight oil reservoirs. Most mainstream commercial reservoir numerical simulation software adopts the conventional orthogonal grid continuous medium model, which is difficult to accurately characterize the complex oblique fractures.

[0004] Although some software have introduced more accurate discrete fracture models for complex fracture networks in recent years, there is still a lack of accurate characterization of reservoir nonlinear seepage characteristics, osmotic pressure and other phenomena. In addition, during the water injection process, when chemicals are added to the injected water or low-mineralization water drive / throws-up and discharge is used, it is also necessary to fully characterize the resulting changes in oil-water interfacial tension, osmotic pressure, rock wettability and other properties that lead to differences in water injection displacement / imbibition efficiency. However, at present, scholars focus more on the numerical simulation characterization of some seepage mechanisms, lacking a systematic and comprehensive characterization of the multiphase and multi-component seepage of the process. The development model of large-scale fracturing and water injection / chemical agent enhancement for ultra-low permeability-tight oil reservoirs brings new challenges to numerical simulation technology.

[0005] In view of the technical problems existing in the above-mentioned prior art, the present invention characterizes complex fracture networks through discrete fracture models, considers reservoir stress sensitivity and nonlinear seepage characteristics, couples the adsorption and diffusion effects of surfactants and salts, and the capillary pressure, phase permeability curve and osmotic pressure changes under the action of the two, and establishes a multiphase and multicomponent seepage mathematical model. The present invention realizes explicit fracture characterization through an adaptive grid partitioning method, adopts a finite volume method for solving, establishes a multiphase and multicomponent seepage model based on a discrete fracture model, and realizes the simulation of ultra-low permeability-tight oil reservoir development under the influence of complex fracture networks.

[0006] The present invention focuses on the systematic and comprehensive characterization and numerical solution of the complex fracture network system and multi-phase and multi-component seepage process in ultra-low permeability and tight oil reservoirs. Summary of the invention

[0007] The purpose of the present invention is to provide a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model.

[0008] The present invention is achieved through the following technical solutions:

[0009] The present invention relates to a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model, comprising the following steps:

[0010] Step 1, establishing a multi-phase and multi-component seepage mathematical model: including the following steps:

[0011] (1.1) Establish the four-component mass conservation equation;

[0012] (1.2) Non-Darcy effects of matrix and cracks;

[0013] (1.3) Adsorption effect of surfactants and salts;

[0014] (1.4) Diffusion effect of surfactants and salts;

[0015] (1.5) Changes in capillary pressure under the action of surfactants and salts;

[0016] (1.6) Changes in relative permeability curves under the action of surfactants and salts;

[0017] (1.7) Osmotic pressure caused by changes in mineralization;

[0018] Then the established multi-phase and multi-component seepage mathematical model based on the discrete fracture model is solved;

[0019] Step 2: Discrete crack model mesh generation:

[0020] Step 3, discretization of mass conservation equation;

[0021] Step 4: Calculation of conductivity between “matrix and crack”;

[0022] Step 5, establish a nonlinear equation system and use automatic differentiation technology and Newton iteration method to complete the solution;

[0023] Step 6, forming a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model.

[0024] Preferably, in step (1.1), the specific steps of establishing the four-component mass conservation equation are:

[0025] Based on the black oil model, a mathematical model of seepage of four components in oil and water phases is established; the four components refer to: oil, water, surfactant and salt; among them, the water component, salt component and surfactant component only exist in the water phase, and the oil component only exists in the oil phase; considering the adsorption and diffusion effects of surfactants and salt, the mass conservation equations of the four components are established respectively, such as formulas (1)-(4):

[0026]

[0027]

[0028] Among them, the seepage velocities of oil phase and water phase are:

[0029]

[0030] Where φ is the porosity, a decimal; s w ,s o are the saturation of water phase and oil phase, respectively, in decimal; B w , B o are the volume coefficients of water phase and oil phase, respectively, m 3 / m 3 ;v w , v o are the water phase and oil phase seepage velocities at the surface ds, m / s; n is the normal vector at the surface ds; q wS ,q oS are the water production flow rate and oil production flow rate of the well under surface conditions, m 3 / s;c s , c b are the mass concentrations of surfactant and salt, kg / m 3 ρ R is the rock density, kg / m 3 ; are the adsorption concentrations of surfactant and salt, kg / kg, respectively; are the diffusion mass fluxes of surfactant and salt, kg / (m 2 ·s); k is the absolute permeability, m 2 ;k ro and k rw are the relative permeabilities of oil phase and water phase, respectively, dimensionless; μ o and μ w are respectively the viscosity of oil phase and water phase, Pa.s; p o and p w are respectively the oil phase and water phase pressure, Pa; p π is the osmotic pressure, Pa.

[0031] The above formulas (1) to (4) respectively give the mass conservation equations of oil, water, surfactant and salt components. Among them, formulas (3) and (4) take into account the adsorption and diffusion effects of surfactants and salts. On the basis of establishing the mass conservation equations, the seepage characteristics of each component are characterized.

[0032] Preferably, in step (1.2), the non-Darcy effect of the matrix and the crack is specifically:

[0033] There is an obvious non-Darcy effect in the matrix and fracture seepage of ultra-low permeability-tight oil reservoirs. The matrix non-Darcy effect includes stress sensitivity effect and nonlinear flow effect; for stress sensitivity effect, the matrix permeability increases exponentially with increasing pressure, see formula (7):

[0034]

[0035] In formula (7), k m0 is the matrix at the initial pressure p m0 Initial permeability, m 2 ; C m is the matrix permeability stress sensitivity coefficient, dimensionless;

[0036] For nonlinear flow effects, the matrix permeability and pressure gradient have a nonlinear relationship, as shown in formulas (8) and (9):

[0037]

[0038] In formulas (8) and (9), β is the nonlinear coefficient; λ is the starting pressure gradient, Pa / m; α n is the nonlinear coefficient, Pa / m; is the pressure gradient, Pa / m;

[0039] Considering stress sensitivity and nonlinear flow, the matrix permeability is given by formula (10):

[0040]

[0041] The non-Darcy effect of fractures only considers the stress sensitivity effect and does not consider the nonlinear flow effect. For the stress sensitivity effect, the fracture permeability increases exponentially with the increase of pressure, as shown in formula (11):

[0042]

[0043] In formula (11), k f0 is the initial pressure of the crack f0 Initial permeability, m 2 ;c f is the stress sensitivity coefficient of fracture permeability, dimensionless.

[0044] Preferably, in step (1.3), the adsorption effect of the surfactant and the salt is specifically:

[0045] The diffusion effect of surfactants and salts is reflected in the adsorption term in the mass conservation equation, specifically the second term in equation (3) and equation (4), as shown in equations (12) and (13):

[0046]

[0047] In formulas (12) and (13), and are respectively the surfactant and salt adsorption mass per unit volume of rock, kg / m 3 ;

[0048] The adsorption concentration of surfactant and salt is related to their aqueous phase concentration. Before simulation calculation, they need to be given in the form of data tables. and relationship; among them, It refers to the adsorption concentration of surfactant and its concentration in water phase; It refers to the adsorbed concentration of the salt and its concentration in the aqueous phase.

[0049] Preferably, in step (1.4), the diffusion effect of the surfactant and the salt is specifically:

[0050] The diffusion effect of surfactants and salts is reflected in the diffusion term in the mass conservation equation, and the relationship between the diffusion mass flux and concentration is established using Fick's diffusion law, as shown in formulas (14) and (15):

[0051]

[0052] In formulas (14) and (15), D s and D b are the surfactant and salt diffusion coefficients, respectively, m 2 / s; It refers to the concentration gradient of surfactant; Refers to the concentration gradient of salt; before simulation calculation, the diffusion coefficients of surfactant and salt in the matrix and cracks need to be defined respectively.

[0053] Preferably, in step (1.5), the capillary pressure change under the action of the surfactant and the salt is specifically:

[0054] The capillary pressure calculation formula is shown in formula (16):

[0055]

[0056] In formula (16), p cis the capillary pressure, Pa; σ is the interfacial tension, N / m; θ is the wetting angle, °; r is the capillary radius, m;

[0057] In a specific core or reservoir, the capillary pressure is affected by both the oil-water interfacial tension and the wetting angle. Surfactants can change the oil-water interfacial tension and rock wettability, and changes in salt concentration can also cause changes in rock wettability. At this time, the capillary pressure calculation is shown in formula (17):

[0058]

[0059] In formula (17), σ0 and θ0 are the interfacial tension and wetting angle when the surfactant concentration is 0 and the salt concentration is the initial concentration, respectively. c0 is the capillary pressure when the surfactant concentration is 0 and the salt concentration is the initial concentration.

[0060] Preferably, in step (1.6), the change of the phase permeability curve under the action of the surfactant and the salt is specifically:

[0061] Changes in surfactant and salt concentrations can alter the oil-water interfacial tension and rock wettability, thereby affecting the oil-water relative permeability curve.

[0062] Under the influence of interfacial tension, the relative permeability of water or oil phase is:

[0063]

[0064] Among them, k ra (s w ,N c ) refers to the relative permeability of the water phase or oil phase under the influence of the final interfacial tension; m is the interpolation factor; It refers to the relative permeability of water or oil phase under the influence of interfacial tension; It refers to the original relative permeability of water phase or oil phase under the influence of interfacial tension; k ra Refers to the relative permeability of water phase or oil phase; w Refers to the water phase permeability; N c It refers to the number of capillaries;

[0065] For the effect of wettability, define the relative permeability curve under water-wet conditions. and The phase permeability interpolation method is also used to characterize the effect of wettability change on the phase permeability curve, namely:

[0066]

[0067] In formula (19), k ra (s w ,c a) refers to the relative permeability of the water phase or oil phase under the influence of rock wettability; Fs refers to the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration; c a refers to the water-wet condition, It refers to the original relative permeability of the water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state; In the wet state (F s =0) phase permeability curve; F s is the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration and is given in a data table;

[0068] Taking into account the influence of interfacial tension and wettability changes on the oil-water phase permeability curve, and combining equations (18) and (19), the relative permeability of the water phase or oil phase is obtained, as shown in equation (20):

[0069]

[0070] In formula (20), k ra refers to the relative permeability of the water phase or oil phase; m refers to the interpolation factor; It refers to the relative permeability of water phase or oil phase under the influence of interfacial tension; F s (c a ) is a function of surfactant or salt adsorption concentration; It refers to the original relative permeability of water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state.

[0071] Preferably, in step (1.7), the osmotic pressure effect caused by the change in mineralization is specifically:

[0072] Due to the attraction of solute particles to water, there will be a certain osmotic pressure between aqueous solutions of different concentrations on both sides of the semi-permeable membrane. In actual reservoirs, the clay mineral diffuse double layer forms a semi-permeable membrane under the action of surface adsorption and diffusion movement, which allows water molecules to pass through and limits the permeation of salt ions;

[0073] The osmotic pressure generated by the difference in ion mineralization is calculated as shown in formula (21):

[0074]

[0075] In formula (21), γ m is the membrane efficiency, a decimal; T is the absolute temperature, K; V is the molar volume of water, which is 1.8×10 -5 m 3 / mol; A sh and A f are the activities of solution and sandstone pore fluid, mol / m 3 .

[0076] Preferably, in step 2, the specific process of meshing the discrete fracture model is as follows: using the idea of ​​dimensionality reduction, taking the fracture as the edge of the triangular mesh; introducing the fracture as a constraint condition, keeping the position of the fracture point unchanged during the meshing process; although the fracture width is not considered during the meshing process, the influence of the fracture width on the seepage process is considered in the actual simulation calculation. After the meshing is completed, the fracture line coordinates are processed into a rectangular mesh with a width, and the mesh connection relationship is reestablished.

[0077] To illustrate the quality of mesh generation, the following formula is usually used for calculation in finite element analysis:

[0078]

[0079] Among them, a, b and c are the side lengths of the triangular meshes respectively. When q>0.5 of all meshes, the generated mesh quality is considered to be good.

[0080] Therefore, the random fracture network generation program is used to change the fracture parameters (number of fractures, average fracture length and fracture direction) to generate the corresponding random fracture network. The meshing algorithm and process mentioned above are used for meshing, and the quality q values ​​of all meshing are calculated to illustrate the adaptability of the meshing algorithm.

[0081] Preferably, in step 3, the specific process of discretizing the mass conservation equation is:

[0082] The mass conservation equation includes: accumulation term, adsorption term, diffusion term, flow term and source-sink term; different discretization methods are used in the discretization process of each term;

[0083] Among them, the accumulation term and adsorption term are discretized using the forward difference format; the flow term is discretized based on the finite volume standard discretization format and non-Darcy seepage law; the diffusion term is discretized based on the finite volume standard discretization format and Fick's law; the source and sink terms are discretized, and the non-Darcy seepage law needs to be considered; the difference equation group of the mass conservation equation of each component is obtained:

[0084]

[0085]

[0086] In the above equations, the subscript i represents the grid block i, is the mesh pore volume, m 3 ; B r is the rock compression coefficient, Pa -1 ; The superscript n represents a known time step, and the superscript n+1 represents a time step to be determined; η cis the set of adjacent grid blocks of grid block i; p is the pressure at the center of the grid block, Pa; k r is the relative permeability on the grid surface ij connecting grid blocks i and j, dimensionless; μ is the viscosity on the grid surface ij connecting grid blocks i and j, Pa.s; q coSi and q cwSi are the oil and water phase flow rates flowing into or out of the well in grid block i, respectively, m 3 / s, WI is the well index of the node unit where the well is located, m 2 ·m;p bh is the bottom hole pressure, Pa.

[0087] Preferably, in step 4, the conductivity between "matrix and crack" is calculated as follows:

[0088] The non-Darcy effect of matrix and fracture in the calculation is reflected in the conductivity T nD As shown in formula (27):

[0089] T nD =γT (27)

[0090] In formula (27), T nD is the non-Darcy conductivity, m 2 ·m; γ is the non-Darcy coefficient, dimensionless; T is the Darcy conductivity, m 2 ·m; Darcy conductivity is calculated as follows:

[0091] ① Matrix-matrix (MMC)

[0092]

[0093] ②Matrix-crack (MFC) (1 is matrix grid, 2 is crack grid)

[0094]

[0095] ③Fracture-fracture (FFC) (calculated using star-delta conversion)

[0096]

[0097] Where A represents the area of ​​the intersection of two grids, m 2 ; D is the distance from the grid to the intersection surface, m; f is the unit direction vector of D; n is the normal vector of the intersection surface.

[0098] In actual calculation, the semi-conductivity t of the matrix and fracture is calculated first, and then the total conductivity is calculated by the harmonic mean method. The calculation of the matrix non-Darcy semi-conductivity is shown in formula (31):

[0099]

[0100] The calculation of the non-Darcy semi-conductivity of the crack is shown in formula (32):

[0101]

[0102] The stress sensitivity effect is related to pressure and can be directly brought into the calculation. However, the nonlinear flow effect involves the pressure gradient of the grid surface, and the calculation is relatively complicated. The present invention adopts the "half-surface upstream weight method + constructing nonlinear flow operator" method for calculation.

[0103] The discretized mass conservation equations are combined to establish a nonlinear equation group; then, the nonlinear equation group is vectorized and fully implicitly processed; due to the fully implicit processing of the equation group, the equation group has extremely strong nonlinearity, so the Newton iteration method is used to solve it; finally, automatic differentiation technology and Newton iteration method are used to complete the solution, forming a multiphase and multicomponent seepage numerical simulation method based on discrete fracture model (MPMC-DFM); the hardware platform for this numerical simulation method is computer.

[0104] The present invention has the following advantages:

[0105] (1) Based on the seepage and development characteristics of ultra-low permeability-tight oil reservoirs, the present invention characterizes complex fracture networks through discrete fracture models, considers reservoir stress sensitivity and nonlinear seepage characteristics, couples the adsorption and diffusion effects of surfactants and salts, as well as the capillary pressure, phase permeability curve and osmotic pressure changes under the action of the two, and establishes a multi-phase and multi-component seepage numerical simulation method under a complex fracture network system, providing effective technical support for the development of ultra-low permeability-tight oil reservoirs.

[0106] (2) The present invention uses this numerical simulation method to simulate the matrix and fracture stress sensitivity, nonlinear seepage, and the effect of surfactants and salt on the fracturing energy enhancement and imbibition of ultra-low permeability-tight oil reservoirs. The actual simulation results show that the development of ultra-low permeability-tight oil reservoirs should take into account the nonlinear seepage characteristics of the reservoir to accurately evaluate the development scope and oil well production capacity; the energy enhancement and imbibition effect can be improved by appropriately reducing the oil-water interfacial tension by adding surfactants; the osmotic pressure caused by low salinity can improve the energy enhancement and imbibition effect to a certain extent, but the oil enhancement effect is limited. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] Figure 1 It is a technical roadmap of the multi-phase and multi-component seepage numerical simulation method based on discrete fracture model involved in the present invention;

[0108] Figure 2 It is a flow chart of the triangular mesh generation algorithm of crack constraint of the present invention;

[0109] Figure 3 is the histogram of mesh generation and mesh quality distribution under different crack parameters;

[0110] Figure 4 It is a schematic diagram of conductivity calculation;

[0111] Figure 5 It is a model diagram of a multi-stage fracturing horizontal well;

[0112] Figure 6 is a comparison diagram of oil-water relative permeability curves in the embodiment;

[0113] Figure 7 is a graph showing the relationship between daily oil production and time under different matrix stress sensitivity coefficients in the embodiment;

[0114] Figure 8 is a graph showing the relationship between cumulative oil production and time under different matrix stress sensitivity coefficients in the embodiment;

[0115] Fig. 9 is a pressure field distribution diagram under different matrix stress sensitivity coefficients in the embodiment;

[0116] Fig.10 is a graph showing the relationship between daily oil production and time under different fracture stress sensitivity coefficients in the embodiment;

[0117] Fig.11 is a graph showing the relationship between cumulative oil production and time under different fracture stress sensitivity coefficients in the embodiment;

[0118] Fig.12 is a pressure field distribution diagram under different crack stress sensitivity coefficients in the embodiment;

[0119] Fig.13 is a graph showing the relationship between daily oil production and time under different matrix starting pressure gradients in the embodiment;

[0120] Fig.14 is a graph showing the relationship between cumulative oil production and time under different matrix starting pressure gradients in the embodiment;

[0121] Fig.15 is a pressure field distribution diagram under different matrix starting pressure gradients in the embodiment;

[0122] Fig.16 is a graph showing the relationship between daily oil production and time under different nonlinear coefficients in the embodiment;

[0123] Fig.17 is a graph showing the relationship between cumulative oil production and time under different nonlinear coefficients in the embodiment;

[0124] Fig.18 3 is a pressure field distribution diagram under different nonlinear coefficients in the embodiment. DETAILED DESCRIPTION

[0125] The present invention is described in detail below in conjunction with specific embodiments. It should be noted that the following embodiments are only further descriptions of the present invention, but the protection scope of the present invention is not limited to the following embodiments.

[0126] Example 1

[0127] This embodiment relates to a multi-phase multi-component seepage numerical simulation method based on a discrete fracture model. Figure 1 As shown, the following steps are included:

[0128] Step 1, establishing a multi-phase and multi-component seepage mathematical model: including the following steps:

[0129] (1.1) Establish the four-component mass conservation equation;

[0130] (1.2) Non-Darcy effects of matrix and cracks;

[0131] (1.3) Adsorption effect of surfactants and salts;

[0132] (1.4) Diffusion effect of surfactants and salts;

[0133] (1.5) Changes in capillary pressure under the action of surfactants and salts;

[0134] (1.6) Changes in relative permeability curves under the action of surfactants and salts;

[0135] (1.7) Osmotic pressure caused by changes in mineralization;

[0136] Then the established multi-phase and multi-component seepage mathematical model based on the discrete fracture model is solved;

[0137] Step 2: Discrete crack model mesh generation:

[0138] Step 3, discretization of mass conservation equation;

[0139] Step 4: Calculation of conductivity between “matrix and crack”;

[0140] Step 5, establish a nonlinear equation system and use automatic differentiation technology and Newton iteration method to complete the solution;

[0141] Step 6, forming a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model.

[0142] Preferably, in step (1.1), the specific steps of establishing the four-component mass conservation equation are:

[0143] Based on the black oil model, a mathematical model of seepage of four components in oil and water phases is established; the four components refer to: oil, water, surfactant and salt; among them, the water component, salt component and surfactant component only exist in the water phase, and the oil component only exists in the oil phase; considering the adsorption and diffusion effects of surfactants and salt, the mass conservation equations of the four components are established respectively, such as formulas (1)-(4):

[0144]

[0145] Among them, the seepage velocities of oil phase and water phase are:

[0146]

[0147] Where φ is the porosity, a decimal; s w ,s o are the saturation of water phase and oil phase, respectively, in decimal; B w , B o are the volume coefficients of water phase and oil phase, respectively, m 3 / m 3 ;v w , v o are the water phase and oil phase seepage velocities at the surface ds, m / s; n is the normal vector at the surface ds; q wS ,q oS are the water production flow rate and oil production flow rate of the well under surface conditions, m 3 / s;c s , c b are the mass concentrations of surfactant and salt, kg / m 3 ρ R is the rock density, kg / m 3 ; are the adsorption concentrations of surfactant and salt, kg / kg, respectively; are the diffusion mass fluxes of surfactant and salt, kg / (m 2 ·s); k is the absolute permeability, m 2 ;k ro and k rw are the relative permeabilities of oil phase and water phase, respectively, dimensionless; μ o and μ w are respectively the viscosity of oil phase and water phase, Pa.s; p o and p w are the oil phase and water phase pressures, Pa respectively; pπ is the osmotic pressure, Pa.

[0148] The above formulas (1) to (4) respectively give the mass conservation equations of oil, water, surfactant and salt components. Among them, formulas (3) and (4) take into account the adsorption and diffusion effects of surfactants and salts. On the basis of establishing the mass conservation equations, the seepage characteristics of each component are characterized.

[0149] Preferably, in step (1.2), the non-Darcy effect of the matrix and the crack is specifically:

[0150] There is an obvious non-Darcy effect in the matrix and fracture seepage of ultra-low permeability-tight oil reservoirs. The matrix non-Darcy effect includes stress sensitivity effect and nonlinear flow effect; for stress sensitivity effect, the matrix permeability increases exponentially with increasing pressure, see formula (7):

[0151]

[0152] In formula (7), k m0 is the matrix at the initial pressure p m0 Initial permeability, m 2 ; C m is the matrix permeability stress sensitivity coefficient, dimensionless;

[0153] For nonlinear flow effects, the matrix permeability and pressure gradient have a nonlinear relationship, as shown in formulas (8) and (9):

[0154]

[0155]

[0156] In formulas (8) and (9), β is the nonlinear coefficient; λ is the starting pressure gradient, Pa / m; α n is the nonlinear coefficient, Pa / m; is the pressure gradient, Pa / m;

[0157] Considering stress sensitivity and nonlinear flow, the matrix permeability is given by formula (10):

[0158]

[0159] The non-Darcy effect of fractures only considers the stress sensitivity effect and does not consider the nonlinear flow effect. For the stress sensitivity effect, the fracture permeability increases exponentially with the increase of pressure, as shown in formula (11):

[0160]

[0161] In formula (11), k f0 is the initial pressure of the crack f0 Initial permeability, m 2 ;c fis the stress sensitivity coefficient of fracture permeability, dimensionless.

[0162] Preferably, in step (1.3), the adsorption effect of the surfactant and the salt is specifically:

[0163] The diffusion effect of surfactants and salts is reflected in the adsorption term in the mass conservation equation, specifically the second term in equation (3) and equation (4), as shown in equations (12) and (13):

[0164]

[0165] In formulas (12) and (13), and are respectively the surfactant and salt adsorption mass per unit volume of rock, kg / m 3 ;

[0166] The adsorption concentration of surfactant and salt is related to their aqueous phase concentration. Before simulation calculation, they need to be given in the form of data tables. and relationship; among them, It refers to the adsorption concentration of surfactant and its concentration in water phase; It refers to the adsorbed concentration of the salt and its concentration in the aqueous phase.

[0167] Preferably, in step (1.4), the diffusion effect of the surfactant and the salt is specifically:

[0168] The diffusion effect of surfactants and salts is reflected in the diffusion term in the mass conservation equation, and the relationship between the diffusion mass flux and concentration is established using Fick's diffusion law, as shown in formulas (14) and (15):

[0169]

[0170] In formulas (14) and (15), D s and D b are the surfactant and salt diffusion coefficients, respectively, m 2 / s; It refers to the concentration gradient of surfactant; Refers to the concentration gradient of salt; before simulation calculation, the diffusion coefficients of surfactant and salt in the matrix and cracks need to be defined respectively.

[0171] Preferably, in step (1.5), the capillary pressure change under the action of the surfactant and the salt is specifically:

[0172] The capillary pressure calculation formula is shown in formula (16):

[0173]

[0174] In formula (16), p c is the capillary pressure, Pa; σ is the interfacial tension, N / m; θ is the wetting angle, °; r is the capillary radius, m;

[0175] In a specific core or reservoir, the capillary pressure is affected by both the oil-water interfacial tension and the wetting angle. Surfactants can change the oil-water interfacial tension and rock wettability, and changes in salt concentration can also cause changes in rock wettability. At this time, the capillary pressure calculation is shown in formula (17):

[0176]

[0177] In formula (17), σ0 and θ0 are the interfacial tension and wetting angle when the surfactant concentration is 0 and the salt concentration is the initial concentration, respectively. c0 is the capillary pressure when the surfactant concentration is 0 and the salt concentration is the initial concentration.

[0178] Preferably, in step (1.6), the change of the phase permeability curve under the action of the surfactant and the salt is specifically:

[0179] Changes in surfactant and salt concentrations can alter the oil-water interfacial tension and rock wettability, thereby affecting the oil-water relative permeability curve.

[0180] Under the influence of interfacial tension, the relative permeability of water or oil phase is:

[0181]

[0182] In formula (18), k ra (s w ,N c ) refers to the relative permeability of the water phase or oil phase under the influence of the final interfacial tension; m is the interpolation factor; It refers to the relative permeability of water or oil phase under the influence of interfacial tension; It refers to the original relative permeability of water phase or oil phase under the influence of interfacial tension; k ra Refers to the relative permeability of water phase or oil phase; w Refers to the water phase permeability; N c It refers to the number of capillaries;

[0183] For the effect of wettability, define the relative permeability curve under water-wet conditions. and The phase permeability interpolation method is also used to characterize the effect of wettability change on the phase permeability curve, namely:

[0184]

[0185] In formula (19), k ra (s w ,c a) refers to the relative permeability of the water phase or oil phase under the influence of rock wettability; F s refers to the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration; c a refers to the water-wet condition, It refers to the original relative permeability of the water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state; In the wet state (F s =0) phase permeability curve; F s is the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration and is given in a data table;

[0186] Taking into account the influence of interfacial tension and wettability changes on the oil-water phase permeability curve, and combining equations (18) and (19), the relative permeability of the water phase or oil phase is obtained, as shown in equation (20):

[0187]

[0188] In formula (20), k ra refers to the relative permeability of the water phase or oil phase; m refers to the interpolation factor; It refers to the relative permeability of water phase or oil phase under the influence of interfacial tension; F s (c a ) is a function of surfactant or salt adsorption concentration; It refers to the original relative permeability of water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state.

[0189] Preferably, in step (1.7), the osmotic pressure effect caused by the change in mineralization is specifically:

[0190] Due to the attraction of solute particles to water, there will be a certain osmotic pressure between aqueous solutions of different concentrations on both sides of the semi-permeable membrane. In actual reservoirs, the clay mineral diffuse double layer forms a semi-permeable membrane under the action of surface adsorption and diffusion movement, which allows water molecules to pass through and limits the permeation of salt ions;

[0191] The osmotic pressure generated by the difference in ion mineralization is calculated as shown in formula (21):

[0192]

[0193] In formula (21), γ m is the membrane efficiency, a decimal; T is the absolute temperature, K; V is the molar volume of water, which is 1.8×10 -5 m 3 / mol; A sh and A f are the activities of solution and sandstone pore fluid, mol / m3 .

[0194] Preferably, in step 2, the specific process of meshing the discrete fracture model is as follows: adopting the idea of ​​dimensionality reduction, taking fractures as the edges of triangular meshes; introducing fractures as constraints to keep the position of fracture points unchanged during meshing; although the fracture width is not considered during meshing, the influence of fracture width on the seepage process is considered in actual simulation calculations. When meshing is completed, the fracture line coordinates are processed into rectangular meshes with widths, and the mesh connection relationship is reestablished. See the meshing process for details. Figure 2 shown.

[0195] To illustrate the quality of mesh generation, the following formula is usually used for calculation in finite element analysis:

[0196]

[0197] Among them, a, b and c are the side lengths of the triangular meshes respectively. When q>0.5 of all meshes, the generated mesh quality is considered to be good.

[0198] Therefore, the random fracture network generation program is used to change the fracture parameters (number of fractures, average fracture length and fracture direction) to generate the corresponding random fracture network. The meshing algorithm and process mentioned above are used for meshing, and the quality q value of all meshing is calculated to illustrate the adaptability of the meshing algorithm. Figure 3 It can be seen that when the crack parameters are arbitrarily changed for triangular meshing, almost all meshes have q>0.5, indicating that the meshing quality is good and the algorithm has good adaptability to different fracture networks.

[0199] Preferably, in step 3, the specific process of discretizing the mass conservation equation is:

[0200] The mass conservation equation includes: accumulation term, adsorption term, diffusion term, flow term and source-sink term; different discretization methods are used in the discretization process of each term;

[0201] Among them, the accumulation term and adsorption term are discretized using the forward difference format; the flow term is discretized based on the finite volume standard discretization format and non-Darcy seepage law; the diffusion term is discretized based on the finite volume standard discretization format and Fick's law; the source and sink terms are discretized, and the non-Darcy seepage law needs to be considered; the difference equation group of the mass conservation equation of each component is obtained:

[0202]

[0203] In the above equations, the subscript i represents the grid block i, is the mesh pore volume, m 3 ; Br is the rock compression coefficient, Pa -1 ; The superscript n represents a known time step, and the superscript n+1 represents a time step to be determined; η c is the set of adjacent grid blocks of grid block i; p is the pressure at the center of the grid block, Pa; k r is the relative permeability on the grid surface ij connecting grid blocks i and j, dimensionless; μ is the viscosity on the grid surface ij connecting grid blocks i and j, Pa.s; q coSi and q cwSi are the oil and water phase flow rates flowing into or out of the well in grid block i, respectively, m 3 / s, WI is the well index of the node unit where the well is located, m 2 ·m;p bh is the bottom hole pressure, Pa.

[0204] Preferably, in step 4, the conductivity between "matrix and crack" is calculated as follows:

[0205] The non-Darcy effect of matrix and fracture in the calculation is reflected in the conductivity T nD As shown in formula (27):

[0206] T nD =γT (27)

[0207] In formula (27), T nD is the non-Darcy conductivity, m 2 ·m; γ is the non-Darcy coefficient, dimensionless; T is the Darcy conductivity, m 2 ·m; Darcy conductivity is calculated as follows:

[0208] ① Matrix-matrix (MMC)

[0209]

[0210] ②Matrix-crack (MFC) (1 is matrix grid, 2 is crack grid)

[0211]

[0212] ③Fracture-fracture (FFC) (calculated using star-delta conversion)

[0213]

[0214] Where A represents the area of ​​the intersection of two grids, m 2 ; D is the distance from the grid to the intersection surface, m; f is the unit direction vector of D; n is the normal vector of the intersection surface. The grid geometric parameters in the above formulas are as follows: Figure 4 .

[0215] In actual calculation, the semi-conductivity t of the matrix and fracture is calculated first, and then the total conductivity is calculated by the harmonic mean method. The calculation of the matrix non-Darcy semi-conductivity is shown in formula (31):

[0216]

[0217] The calculation of the non-Darcy semi-conductivity of the crack is shown in formula (32):

[0218]

[0219] The stress sensitivity effect is related to pressure and can be directly brought into the calculation. However, the nonlinear flow effect involves the pressure gradient of the grid surface, and the calculation is relatively complicated. The present invention adopts the "half-surface upstream weight method + constructing nonlinear flow operator" method for calculation.

[0220] Example 2

[0221] This example studies the influence of various parameters (i.e., stress sensitivity, nonlinear seepage, surfactant, and salt) in the multiphase and multicomponent seepage numerical simulation method based on discrete fracture model (MPMC-DFM) on the development effect of ultra-low permeability-tight oil reservoirs, and establishes a multi-stage fracturing horizontal well model for ultra-low permeability-tight oil reservoirs. Figure 5 As shown in the figure (the vertical direction is the artificial fracture, and the horizontal direction is the Shuping wellbore). The basic parameters of the model are shown in Table 1. The model uses an unstructured triangular grid, and the horizontal well is produced at a constant rate and limited pressure (daily liquid production of 200m 3 / d, bottom hole pressure 15MPa); among them, the oil-water relative permeability curve comparison is shown in Figure 6 shown.

[0222] Table 1

[0223]

[0224] 1. Influence of stress sensitivity

[0225] 1.1 Effect of matrix stress sensitivity

[0226] When the matrix starting pressure gradient and nonlinear coefficient are 0.2MPa / m and 0.05MPa / m respectively, the fracture stress sensitivity coefficient is 0.1MPa -1 , given different matrix stress sensitivity coefficients (0MPa -1 、0.005MPa -1 、0.01MPa -1 and 0.05MPa -1 ), the simulation results of daily oil production and cumulative oil production are shown in Figure 7 and Figure 8As the matrix stress sensitivity coefficient increases, the daily oil production decreases and the cumulative oil production decreases significantly. When the matrix stress sensitivity coefficient is greater than 0.01MPa -1 , daily oil production and cumulative oil production dropped significantly. In the process of depletion development, the formation pressure continued to drop, the reservoir permeability decreased, the oil production decreased, and as the stress sensitivity coefficient increased, the permeability changed more significantly with the formation pressure, and the oil production decline increased, which could be reduced by 36.35%.

[0227] From the pressure field distribution (see Fig. 9 It can also be seen that after depletion development, the pressure at the bottom of the well and near the artificial fractures is the lowest and extends to the deep part of the reservoir. As the matrix stress sensitivity coefficient increases, the pressure field spreads to the deep part of the reservoir decreases, the effective production area decreases, and the corresponding production capacity also decreases. Therefore, when the reservoir matrix permeability has strong stress sensitivity, the impact of the permeability change caused by the decrease in formation pressure on the production well capacity should be taken seriously.

[0228] 1.2 Influence of crack stress sensitivity

[0229] When the matrix stress sensitivity coefficient is 0.01MPa -1 When different crack stress sensitivity coefficients (0MPa -1 、0.05MPa -1 , 0.1MPa -1 and 0.2MPa -1 ), the simulation results of daily oil production and cumulative oil production are shown in Fig.10 and Fig.11 As the fracture stress sensitivity coefficient increases, the daily oil production and cumulative oil production both show a significant decrease. When the fracture stress sensitivity coefficient is greater than 0.1 MPa -1 , the daily oil production and cumulative oil production decreased significantly. Similarly, in the process of depletion development, as the formation pressure continues to decrease, the reservoir fractures gradually close, the width decreases, the conductivity decreases, and as the stress sensitivity of the fractures increases, the conductivity decreases more significantly, and the oil production decreases more, which can be reduced by 54.51%. From the pressure field distribution after depletion development (see Fig.12 It can also be seen that as the stress sensitivity coefficient of the fracture increases, the pressure field spreads to the deep reservoir, the effective production area decreases, and the production capacity also decreases. The fracture permeability changes with pressure, which is not only related to the fracture opening and roughness, but also to the size, concentration and strength of the proppant in the fracture. The fracture opening that matches the reservoir should be designed and the appropriate proppant should be selected to avoid the problem of a significant decrease in production capacity due to a sharp decrease in fracture permeability or even closure as the pressure decreases.

[0230] 2. Influence of nonlinear seepage

[0231] Due to the poor development of matrix physical properties in ultra-low permeability-tight reservoirs, the pore size is mostly submicron and nanometer, and the fluid seepage shows nonlinear characteristics. In order to explore the influence of matrix nonlinear seepage on the development effect of ultra-low permeability-tight reservoirs, the production capacity and reservoir pressure field changes of multi-stage fracturing horizontal wells were investigated under different matrix starting pressure gradient λ and nonlinear coefficient α.

[0232] 2.1 Effect of startup pressure gradient

[0233] When the matrix and crack stress sensitivity coefficients are 0.01MPa -1 and 0.1MPa -1 When the matrix nonlinear coefficient is 0.05MPa / m, different matrix starting pressure gradients (0.06MPa / m, 0.1MPa / m, 0.2MPa / m, 0.4MPa / m and 0.8MPa / m) are given, the simulation results of daily oil production and cumulative oil production are shown in Fig.13 and Fig.14 . Considering nonlinear seepage, the production capacity of the production well will drop significantly, and as the starting pressure gradient increases, the daily oil production and cumulative oil production will show a significant decline, and the cumulative oil production decline can reach 77.80%. Different from linear seepage, the existence of the starting pressure gradient causes the fluid to produce additional seepage resistance, the pressure drop is accelerated, and the production capacity of the production well decreases. From the pressure field distribution after exhaustion development (such as Fig.15 As shown in the figure, it can be seen that without considering nonlinear seepage, the production pressure of the production well spreads faster and the production range is larger; however, as the starting pressure gradient increases, the pressure field spreads to the deep reservoir and the effective production area decreases, and the production capacity also decreases. Therefore, for ultra-low permeability-tight reservoirs, due to the small pore throat size, obvious micro-scale boundary layer effect, and obvious nonlinear seepage characteristics, the impact of nonlinear seepage on production capacity should be considered in production and development.

[0234] 2.2 Influence of nonlinear coefficient

[0235] When the matrix and crack stress sensitivity coefficients are 0.01MPa -1 and 0.1MPa -1 When the matrix starting pressure gradient is 0.2MPa / m, different nonlinear coefficients (0.01MPa / m, 0.05MPa / m, 0.1MPa / m and 0.18MPa / m) are given, the simulation results of daily oil production and cumulative oil production are shown in Fig.16 and Fig.17 . Considering the nonlinear seepage, the production capacity of the production well will decrease significantly, and as the nonlinear coefficient increases, the daily oil production and cumulative oil production will both decrease. The larger the nonlinear coefficient, the longer the nonlinear seepage section. Under the same displacement pressure gradient, the smaller the seepage velocity and the lower the production capacity of the production well. From the pressure field distribution after depletion development (such as Fig.18It can also be seen from the figure that without considering nonlinear seepage, the production pressure of the production well diffuses faster and has a larger range of production; but as the nonlinear coefficient increases, the pressure field spreads to the deep reservoir and the effective production area decreases, and the production capacity also decreases. Compared with the starting pressure gradient, the nonlinear coefficient has a relatively weak effect on the production capacity and pressure field. In production development, it is necessary to test and evaluate the relationship between the seepage velocity and displacement pressure gradient of reservoir rock samples, so as to obtain nonlinear seepage parameters and provide basic parameters for production capacity evaluation.

[0236] The present invention characterizes complex fracture networks through discrete fracture models, considers reservoir stress sensitivity and nonlinear seepage characteristics, couples the adsorption and diffusion effects of surfactants and salts, and the capillary pressure, phase permeability curve and osmotic pressure changes under the action of the two, and establishes a multiphase and multicomponent seepage mathematical model. The fracture explicit characterization is achieved through an adaptive grid partitioning method, and the finite volume method is used for solution to establish a multiphase and multicomponent seepage model based on the discrete fracture model, thereby realizing the simulation of the development of ultra-low permeability-tight oil reservoirs under the influence of complex fracture networks.

[0237] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art may make various modifications or variations within the scope of the claims, which will not affect the essence of the present invention.

Claims

1. A multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model, characterized in that: The following steps are involved: Step 1, establishing a multi-phase and multi-component seepage mathematical model: including the following steps: (1.1) Establish the four-component mass conservation equation; (1.2) Non-Darcy effects of matrix and cracks; (1.3) Adsorption effect of surfactants and salts; (1.4) Diffusion effect of surfactants and salts; (1.5) Changes in capillary pressure under the action of surfactants and salts; (1.6) Changes in relative permeability curves under the action of surfactants and salts; (1.7) Osmotic pressure caused by changes in mineralization; Then the established multi-phase and multi-component seepage mathematical model based on the discrete fracture model is solved; Step 2: Discrete crack model mesh generation: Step 3, discretization of mass conservation equation; Step 4, calculation of conductivity between "matrix-crack"; Step 5, establish a nonlinear equation system and use automatic differentiation technology and Newton iteration method to complete the solution; Step 6, forming a multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model.

2. The multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model according to claim 1, characterized in that: In step (1.1), the specific steps of establishing the four-component mass conservation equation are: Based on the black oil model, a mathematical model of seepage of four components in oil and water phases is established; the four components refer to: oil, water, surfactant and salt; among them, the water component, salt component and surfactant component only exist in the water phase, and the oil component only exists in the oil phase; the mass conservation equations of the four components are established respectively, such as formulas (1)-(4): Among them, the seepage velocities of oil phase and water phase are: Where φ is the porosity, a decimal; s w ,s o are the saturation of water phase and oil phase, respectively, in decimal; B w , B o are the volume coefficients of water phase and oil phase, respectively, m 3 / m 3 ;v w , v o are the water phase and oil phase seepage velocities at the surface ds, m / s; n is the normal vector at the surface ds; q wS ,q oS are the water production flow rate and oil production flow rate of the well under surface conditions, m 3 / s;c s , c b are the mass concentrations of surfactant and salt, kg / m 3 ρ R is the rock density, kg / m 3 ; are the adsorption concentrations of surfactant and salt, kg / kg, respectively; are the diffusion mass fluxes of surfactant and salt, kg / (m 2 ·s); k is the absolute permeability, m 2 ;k ro , k rw are the relative permeabilities of oil phase and water phase, respectively, dimensionless; μ o , μ w are respectively the viscosity of oil phase and water phase, Pa.s; p o , p w are respectively the oil phase and water phase pressure, Pa; p π is the osmotic pressure, Pa.

3. The multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model according to claim 1, characterized in that: In step (1.2), the non-Darcy effects of the matrix and the cracks are specifically: The matrix non-Darcy effect includes stress sensitivity effect and nonlinear flow effect; For the stress sensitivity effect, the matrix permeability increases exponentially with increasing pressure, see formula (7): In formula (7), k m0 is the matrix at the initial pressure p m0 Initial permeability, m 2 ; C m is the matrix permeability stress sensitivity coefficient, dimensionless; For nonlinear flow effects, the matrix permeability and pressure gradient have a nonlinear relationship, as shown in formulas (8) and (9): In formulas (8) and (9), β is the nonlinear coefficient; λ is the starting pressure gradient, Pa / m; α n is the nonlinear coefficient, Pa / m; is the pressure gradient, Pa / m; Considering stress sensitivity and nonlinear flow, the matrix permeability is given by formula (10): For the stress sensitivity effect, the fracture permeability increases exponentially with increasing pressure, as shown in formula (11): In formula (11), k f0 is the initial pressure of the crack f0 Initial permeability, m 2 ; c f is the stress sensitivity coefficient of fracture permeability, dimensionless.

4. The method for numerical simulation of multi-phase and multi-component seepage based on discrete fracture model according to claim 1, characterized in that: In step (1.3), the adsorption effect of the surfactant and the salt is specifically: The diffusion effect of surfactants and salts is reflected in the adsorption term in the mass conservation equation, specifically the second term in equation (3) and equation (4), as shown in equations (12) and (13): In formulas (12) and (13), and are respectively the surfactant and salt adsorption mass per unit volume of rock, kg / m 3 ; The adsorption concentration of surfactant and salt is related to their aqueous phase concentration. Before simulation calculation, they need to be given in the form of data tables. and relationship; among them, It refers to the adsorption concentration of surfactant and its concentration in water phase; It refers to the adsorbed concentration of the salt and its concentration in the aqueous phase.

5. The multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model according to claim 1, characterized in that: In step (1.4), the diffusion effect of the surfactant and the salt is specifically: The diffusion effect of surfactants and salts is reflected in the diffusion term in the mass conservation equation, and the relationship between the diffusion mass flux and concentration is established using Fick's diffusion law, as shown in formulas (14) and (15): In formulas (14) and (15), D s , D b are the surfactant and salt diffusion coefficients, respectively, m 2 / s; It refers to the concentration gradient of surfactant; Refers to the concentration gradient of salt; before simulation calculation, the diffusion coefficients of surfactant and salt in the matrix and cracks need to be defined respectively.

6. The multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model according to claim 1, characterized in that: In step (1.5), the capillary pressure change under the action of the surfactant and the salt is specifically: The capillary pressure calculation formula is shown in formula (16): In formula (16), p c is the capillary pressure, Pa; σ is the interfacial tension, N / m; θ is the wetting angle, °; r is the capillary radius, m; In a specific core or reservoir, the capillary pressure is affected by both the oil-water interfacial tension and the wetting angle. Surfactants can change the oil-water interfacial tension and rock wettability, and changes in salt concentration can also cause changes in rock wettability. At this time, the capillary pressure calculation is shown in formula (17): In formula (17), σ0 and θ0 are the interfacial tension and wetting angle when the surfactant concentration is 0 and the salt concentration is the initial concentration, respectively. c0 is the capillary pressure when the surfactant concentration is 0 and the salt concentration is the initial concentration.

7. The method for numerical simulation of multi-phase and multi-component seepage based on discrete fracture model according to claim 1, characterized in that: In step (1.6), the change of the phase permeability curve under the action of the surfactant and the salt is specifically as follows: Under the influence of interfacial tension, the relative permeability of water or oil phase is: In formula (18), k ra (s w ,N c ) refers to the relative permeability of the water phase or oil phase under the influence of the final interfacial tension; m is the interpolation factor; It refers to the relative permeability of water or oil phase under the influence of interfacial tension; It refers to the original relative permeability of water phase or oil phase under the influence of interfacial tension; k ra Refers to the relative permeability of water phase or oil phase; w Refers to the water phase permeability; N c It refers to the number of capillaries; For the effect of wettability, define the relative permeability curve under water-wet conditions. and The phase permeability interpolation method is also used to characterize the effect of wettability change on the phase permeability curve, namely: In formula (19), It refers to the relative permeability of water phase or oil phase under the influence of rock wettability; F s refers to the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration; c a refers to the water-wet condition, It refers to the original relative permeability of the water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state; In the wet state (F s =0) phase permeability curve; F s is the interpolation coefficient, which is a function of the surfactant or salt adsorption concentration and is given in a data table; Taking into account the influence of interfacial tension and wettability changes on the oil-water phase permeability curve, and combining equations (18) and (19), the relative permeability of the water phase or oil phase is obtained, as shown in equation (20): In formula (20), k ra refers to the relative permeability of the water phase or oil phase; m refers to the interpolation factor; It refers to the relative permeability of water phase or oil phase under the influence of interfacial tension; F s (c a ) is a function of surfactant or salt adsorption concentration; It refers to the original relative permeability of water phase or oil phase under the influence of wettability; It refers to the relative permeability curve under water-wet state.

8. The multi-phase and multi-component seepage numerical simulation method based on a discrete fracture model according to claim 1, characterized in that: In step (1.7), the osmotic pressure effect caused by the change in mineralization is specifically: The osmotic pressure generated by the difference in ion mineralization is calculated as shown in formula (21): In formula (21), γ m is the membrane efficiency, a decimal; T is the absolute temperature, K; V is the molar volume of water, which is 1.8×10 -5 m 3 / mol; A sh and A f are the activities of solution and sandstone pore fluid, mol / m 3 .

9. The method for numerical simulation of multi-phase and multi-component seepage based on discrete fracture model according to claim 1, characterized in that: In step 2, the specific process of meshing the discrete crack model is as follows: adopting the idea of ​​dimensionality reduction, taking the crack as the edge of the triangular mesh; introducing the crack as a constraint condition, keeping the position of the crack point unchanged during the meshing process; after the meshing is completed, according to the coordinates of the crack line, processing it into a rectangular mesh with a width, and re-establishing the mesh connection relationship; To prove the quality of mesh generation, the following formula (22) is usually used in finite element analysis: Among them, a, b and c are the side lengths of the triangular meshes respectively. When q>0.5 of all meshes, the generated mesh quality is considered to be good.

10. The multi-phase and multi-component seepage numerical simulation method based on discrete fracture model according to claim 1, characterized in that: In step 3, the specific process of discretizing the mass conservation equation is: The mass conservation equation includes: accumulation term, adsorption term, diffusion term, flow term and source-sink term; different discretization methods are used in the discretization process of each term; Among them, the accumulation term and adsorption term are discretized using the forward difference format; the flow term is discretized based on the finite volume standard discretization format and non-Darcy seepage law; the diffusion term is discretized based on the finite volume standard discretization format and Fick's law; the source and sink terms are discretized, and the non-Darcy seepage law needs to be considered; the difference equation group of the mass conservation equation of each component is obtained: In the above equations, the subscript i represents the grid block i, is the mesh pore volume, m 3 ; B r is the rock compression coefficient, Pa -1 ; The superscript n represents a known time step, and the superscript n+1 represents a time step to be determined; η c is the set of adjacent grid blocks of grid block i; p is the pressure at the center of the grid block, Pa; k r is the relative permeability on the grid surface ij connecting grid blocks i and j, dimensionless; μ is the viscosity on the grid surface ij connecting grid blocks i and j, Pa.s; q coSi and q cwSi are the oil and water phase flow rates flowing into or out of the well in grid block i, respectively, m 3 / s, WI is the well index of the node unit where the well is located, m 2 ·m;p bh is the bottom hole pressure, Pa.

11. The multi-phase and multi-component seepage numerical simulation method based on discrete fracture model according to claim 1, characterized in that: In step 4, the conductivity between the matrix and the crack is calculated as follows: The non-Darcy effect of matrix and fracture in the calculation is reflected in the conductivity T nD As shown in formula (27): T nD =γT (27) In formula (27), T nD is the non-Darcy conductivity, m 2 ·m; γ is the non-Darcy coefficient, dimensionless; T is the Darcy conductivity, m 2 ·m; Darcy conductivity is calculated as follows: ① Matrix-matrix (MMC) ②Matrix-crack (MFC) (1 is matrix grid, 2 is crack grid) ③Fracture-fracture (FFC) (calculated using star-delta conversion) Where A represents the area of ​​the intersection of two grids, m 2 ; D is the distance from the grid to the intersection surface, m; f is the unit direction vector of D; n is the normal vector of the intersection surface. In actual calculation, the semi-conductivity t of the matrix and fracture is calculated first, and then the total conductivity is calculated by the harmonic mean method. The calculation of the matrix non-Darcy semi-conductivity is shown in formula (31): The calculation of the non-Darcy semi-conductivity of the crack is shown in formula (32): The stress sensitivity effect is related to pressure and can be directly brought into the calculation.

Images