CN117077577A - Rapid simulation and optimization method suitable for low-permeability fractured reservoir - Google Patents

Abstract

The invention discloses a rapid simulation and optimization method suitable for a low-permeability fractured reservoir, which belongs to the field of unconventional reservoir numerical simulation and comprises the following steps: step 1, abstracting an oil reservoir system into a matrix grid and a crack grid, and constructing a dual medium model according to seepage characteristics of the matrix and the crack; step 2, establishing an automatic history fit mathematical model by utilizing an ES-MDA history fit algorithm and combining constraint conditions; and 3, taking the economic net present value as an objective function of production optimization, establishing an oil reservoir injection and production optimization model by restraining the geological reserves of the oil reservoirs, the pressure of each well and the upper and lower limits of the injection and production liquid, and obtaining a production regulation system for maximizing the economic net present value of the oil reservoirs by adopting a differential evolution algorithm. The invention is helpful for better understanding the geological features and fracture properties of the low permeability oil reservoir, and improves the efficiency and reduces the cost, thereby increasing the economic benefit of oil field development.

Description

Rapid simulation and optimization method suitable for low-permeability fractured reservoir

Technical Field

The invention belongs to the field of unconventional oil reservoir numerical simulation, and particularly relates to a rapid simulation and optimization method suitable for a low-permeability fractured oil reservoir.

Background

Along with the gradual exhaustion of conventional oil and gas resources, the development status of low-permeability fractured reservoirs is increasingly improved, and the low-permeability fractured reservoirs are the main reservoir types for increasing the storage and the production in the future, so how to develop the low-permeability fractured reservoirs efficiently is the important issue of coping with energy crisis in China. Unlike conventional sandstone oil reservoirs, the low-permeability fractured oil reservoir has small gap roar and large seepage resistance, and a starting pressure gradient phenomenon exists in the fluid flowing process, so that the traditional Darcy seepage equation is not suitable, and the seepage channels of the oil reservoirs are increased in a fracturing mode, so that the fractures in the oil reservoirs are fully developed, a single medium model is not suitable for simulation, a dual medium model and a discrete fracture model are adopted in the conventional common model, the two modes can well describe a fracture system, but the discrete fracture model needs to provide specific parameters of the fractures, such as ductility, fracture opening, fracture spacing and the like, and the parameters are relatively difficult to obtain, so that the low-permeability fractured oil reservoir is modeled by adopting the dual medium model and subsequent production optimization is carried out.

Disclosure of Invention

In order to solve the problems, the invention provides a rapid simulation and optimization method suitable for a low-permeability fractured reservoir, which is characterized in that firstly, a series of well-to-well connected communicating pipelines are divided according to a well distribution mode of the reservoir, grids are divided on each pipeline, a divided grid matrix is abstracted into a matrix grid and a fractured grid, each grid has own parameters such as permeability, effective thickness, water saturation and the like, so that a double-hole double-permeability model is constructed, and accurate description and depiction of flow numerical simulation of matrix and fluid in the fracture in the low-permeability fractured reservoir are effectively realized. And meanwhile, the ES-MDA algorithm is utilized to fit the production history of the oil field, the modeled model is updated in real time, an optimization theory is combined, the economic net present value or the accumulated oil production is taken as an objective function, the pressure and the geological reserves of each well and the block are taken as constraints, the subsequent injection and production regulation and control are carried out according to the production history of each well, and the optimization result has guiding significance for the increase of the oil field storage, the upper production and the water and oil control.

The technical scheme of the invention is as follows:

a rapid simulation and optimization method suitable for a low-permeability fractured reservoir comprises the following steps:

step 1, abstracting an oil reservoir system into a matrix grid and a crack grid, and constructing a dual medium model according to seepage characteristics of the matrix and the crack;

step 2, establishing an automatic history fit mathematical model by utilizing an ES-MDA history fit algorithm and combining constraint conditions;

and 3, taking the economic net present value as an objective function of production optimization, establishing an oil reservoir injection and production optimization model by restraining the geological reserves of the oil reservoirs, the pressure of each well and the upper and lower limits of the injection and production liquid, and obtaining a production regulation system for maximizing the economic net present value of the oil reservoirs by adopting a differential evolution algorithm.

Further, the specific process of step 1 is as follows:

step 1.1, dispersing a reservoir into a series of nodes according to well point distribution in an oil reservoir system, and connecting all the nodes pairwise to form an initial well pattern model;

step 1.2, screening an initial well pattern model by introducing judging conditions, generating well connecting pipelines connected in pairs, and dividing the volume of a reservoir into each connecting pipeline;

step 1.3, dividing grids on each connecting pipeline into matrix grids and crack grids, wherein all the matrix grids form a matrix system, and all the crack grids form a crack system; the inter-grid fluid flow capacity is determined by the magnitude of the inter-grid conductivity;

Step 1.4, considering non-Darcy seepage in a matrix system and considering stress sensitivity effect of rock in a fracture system;

step 1.5, constructing a low-permeability fractured reservoir dual medium model by combining seepage rules of the matrix and fluid in the fracture; the dual medium model comprises a continuity equation of a matrix system and a continuity equation of a crack system;

and 1.6, performing linear interpolation approximation on the continuity equation, and fully implicitly solving parameters in the crack and matrix system by using a Newton-Raphson method.

Further, in step 1.2, the determination conditions include a distance condition, an angle condition, and a maximum connection point number condition; wherein the threshold value of the distance condition is 3 times of the minimum well distance, and wells exceeding the threshold value are regarded as unconnected; the angle conditions are specifically set as follows: in a triangle formed by any three connecting pipelines, the minimum internal angle cannot be smaller than the set minimum angle; the maximum connection point number condition is specifically set as follows: each well point allows for a number of well points of 4-8 connections.

Further, in step 1.3, the grid division is that 5-10 grids are equally divided in each connecting pipeline, and when the grid volume of the well point is specified, the calculation formula of the default value of the grid volume of the well point is as follows:

(1)；

Wherein,is +.>The mesh volume where the mesh is located; />Is the total connection volume; />The number of the connecting pipelines;the number of grids in each connecting pipeline; />The number of well points;

after well point grid volume is determined, the calculation formula of the grid volume in the connecting pipeline is as follows:

(2)；

wherein,、/>different well points positioned at two ends of the connecting pipeline; />Indicating the%>A grid; />Representation->、/>The connection pipe of the well point is +.>The volume of the mesh; />Is->、/>The volume of connecting pipes between well points; />Andwell points->And (4) well point->The mesh volume where the mesh is located; />Is +.>The number of pipes connected; />Is +.>The number of pipes connected;

matrix grid and fracture grid, fracture grid volumeMatrix mesh volume->The calculation is as follows:

(3)；

(4)；

(5)；

(6)；

wherein,representing the number of crack spacing units in the coarse grid; />Representing the volume fraction of microcracks in the whole unit; />The micro-crack spacing in the oil reservoir; />The opening degree of micro cracks in an oil reservoir system is set;

the calculation formula of the conductivity between grids is as follows:

(7)；

wherein,for grid->And grid->An inter-conductivity; />For grid->And grid->Is a contact area of (a); />For grid->Center to grid->And grid->Is a distance of the contact surface; / >For grid->Center to grid->And grid->Is a distance of the contact surface; />For grid->Pseudo-permeability->And grid->Pseudo-permeability->Average value of (2);

fluid flow exists between the matrix and the fracture, and the fluid channeling is calculated as follows:

(8)；

wherein,the cross flow between the matrix cracks; />Is a form factor; />Is the fluid density; />Permeability for matrix network; />Is the viscosity of the fluid; />Is the matrix system pressure; />Is fracture system pressure; />Is the matrix grid volume.

Further, in step 1.4, the equation of motion of the matrix system after considering the fidaxy seepage is:

(9)；

wherein,is the matrix seepage velocity; />Permeability for matrix network; />Is the viscosity of the fluid; />Is the flow differential pressure; />And->Two different nonlinear seepage coefficients;

after the stress sensitivity effect of the rock is considered, the motion equation of the fracture system is as follows:

(10)；

wherein,is the crack seepage velocity; />Permeability for the fracture grid; />Is a constant; />Is the viscosity of the fluid;is the stress sensitivity coefficient; />The grid pressure is the crack; />Initial pressure for the fracture grid; />For flow pressure difference；

The formula (9) and the formula (10) are expressed in a unified form as the following equation of motion:

(11)；

wherein, Is a certain fluid; />Representing->Is oil phase->Or water phase->；/>Is->Is not less than the percolation rate; />Is->Is to be determined; />Is->Is a fluid viscosity of (a); />Is the flow differential pressure;

the pseudo-permeability of the matrix and fracture, considering non-darcy seepage and stress-sensitive effects, is expressed as follows:

(12)；

(13)；

wherein,and->Respectively representing a matrix and a crack; />For grid->And grid->Inter-consideration of matrix after Fodaxil seepage->Is to be determined; />Respectively is grid->And grid->Stress-sensitive post-crack->The pseudo permeability of the inner; />For grid->And grid->Interstitial matrix->Is a natural gas; />Representing grid->And grid->Crack between->Is a natural gas; />Is the viscosity of the fluid; />Is a matrix inner grid->And grid->A pressure gradient of the two nodes; />And->Two different nonlinear seepage coefficients; />Is the stress sensitivity coefficient; />For grid->And grid->The current pressure of the inter-fracture; />For grid->And grid->Initial pressure of the inter-fissures.

Further, in step 1.5, the continuity equation of the matrix system is:

(18)；

(19)；

wherein,is oil phase density; />Porosity of the matrix; />Saturation of the oil phase in the matrix; />Is a matrix grid volume; />Is a time step number; / >Representing the divergence; />The seepage velocity of the oil phase in the matrix; />The change value of the volume flow of the oil phase fluid in unit volume of unit time in the matrix grid; />Is the density of the water phase; />Saturation of the aqueous phase within the matrix; />The seepage velocity of the water phase in the matrix; />The change value of the volume flow of the aqueous phase fluid in unit volume of unit time in the matrix grid;

the continuity equation of the fracture system is:

(20)；

(21)；

wherein,porosity of the fracture; />Saturation of the oil phase in the fracture; />Is the seepage velocity of the oil phase in the cracks; />The change value of the volume flow of the oil phase fluid in unit volume of unit time in the fracture grid is obtained; />Saturation of the aqueous phase in the fracture; />Is the seepage velocity of the water phase in the crack; />The change value of the volume flow of the aqueous phase fluid in the unit volume of the unit time of the fracture grid is obtained.

Further, the specific process of step 1.6 is as follows:

step 1.6.1, defining fluidity coefficients as follows:

(22)；

wherein,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>Is the fluidity of a certain phase fluid; />Is the relative permeability of a fluid of a certain phase; />Is the viscosity of a certain fluid;

step 1.6.2, performing linear interpolation approximation on the matrix system and the fracture system by adopting the same format, wherein the linear interpolation approximation adopts the following numerical discrete equation:

(23)；

In the method, in the process of the invention,is the density of a fluid in a certain phase; />Is porosity; />Saturation for a fluid of a certain phase; />Is a mesh volume; />The last moment in the simulation iteration process; />The current time in the simulation iteration process is; />Is the time step; />For grid->A set of adjacent cells; subscript->Representing two adjacent meshes->And->An average value of a parameter on the interface; />For grid->And grid->An inter-conductivity; />For a certain fluid>Time grid->Is a pressure of (2);for a certain fluid>Time grid->Is a pressure of (2); />For grid->Fluid in certain phase at->Injection and production rate at moment;

step 1.6.3, writing the numerical discrete equation into the residual form as follows:

(24)；

in the method, in the process of the invention,equation residual error;

and solving the pressure and water saturation parameters in the fracture and matrix system by adopting a Newton-Paphson iteration solving method, wherein the solving process is as follows: starting from the initial point, a new point approaching to the zero point of the function is obtained through continuous iterative computation, the new point is obtained through Taylor expansion computation at the current point, and then the new point is used as the initial point of the next iteration to continue iterative computation until a solution meeting the requirement is found.

Further, the specific process of step 2 is as follows:

Step 2.1, constructing an objective function of an automatic history fit mathematical model, as follows:

(25)；

in the method, in the process of the invention,is the objective function value; />Representing a model calculation result; />Representing actual observed data;a covariance matrix representing an observation error;

step 2.2, in the history fitting process of the automatic history fitting mathematical model, only the volume of each pipeline, the conductivity between grids and the parameters of the power law phase permeability curve need to be adjusted;

the parameters of the power law phase permeation curve are adjusted by the following power law phase permeation model:

(26)；

wherein,

(27)；

in the method, in the process of the invention,and->Respectively representing the relative permeability of the oil phase and the water phase; />Representing dimensionless saturation; />Anddetermining the bending degree of the relative permeability curves of the water phase and the oil phase respectively; />Representing irreducible water saturation; />Representing residual oil saturation; />Representing the relative permeability value of the aqueous phase at a water saturation of irreducible water saturation; />Is the water saturation value;

step 2.3, parameters for history fitting adjustmentThe following are provided:

(28)；

wherein,for grid->And grid->The inter-conductivity comprises a matrix and a matrix grid, a crack and a crack grid and a matrix and a crack grid; />Is->、/>Matrix-conducting volume of the communication conduit of the well point, < ->Is->、/>Fracture conduction volume of the communication tubing at the well point; / >、/>、/>Is->、/>Different parameters of the communication pipeline phase permeation curve of the well point; />And->Two different nonlinear seepage coefficients;

when all areHistory fit parameters when sharing a set of phase-to-phase data between pairs of wellsThe expression is as follows:

(29)；

wherein,、/>、/>three parameters of the phase permeability curve shared by all well pairs;

step 2.4, in the history fitting process, adding constraint to model parameters; constraints of the automatic history fit mathematical model are as follows:

(30)；

wherein,for grid->And grid->Is defined by the volume of matrix; />Is the total volume of the matrix system; />For grid->And grid->Is defined by the fracture volume of (2); />Is the total volume of the fracture system; />Is the total well count.

Further, the specific process of step 3 is as follows:

and 3.1, adopting an economic net present value as an objective function of production optimization, wherein the mathematical expression is as follows:

(31)；

wherein,is an economic net present value; />Is a control variable; />For time step number, ++>Representing a total number of time steps; />Indicating production well number, +.>Representing the number of production wells; />Indicating the serial number of the water injection well>Representing the number of water injection wells; />Price per ton of crude oil; />For producing wells->At->Oil extraction of the time step; />Cost per ton of sewage to be treated; />Representing a production well- >At->The water yield of the time step; />Representing water injection costs per ton of water; />Indicating water injection well->At->The water injection amount of the time step;

and 3.2, in the production optimization process, controlling the oil reservoir parameters by adopting boundary constraint conditions, wherein the constraint conditions of the oil reservoir parameters are as follows:

(33)；

wherein,is a water injection well->At->The water injection amount of the time step; />For producing wells->At->The liquid extraction amount of the time step;is a filling and production well->At->Time step injection of liquid production volume->Lower limit of (2); />Is a filling and production well->At->Time step injection of liquid production volume->Upper limit of (2); />Is the total geological reserve; />And->Injection well and production well respectively>Control pressure->Lower and upper limits of (2);

and 3.3, obtaining a production regulation system for maximizing the economic net present value of the oil reservoir by adopting a differential evolution algorithm.

Further, the specific process of step 3.3 is as follows: first, a series of control variables are initializedThe series of control variables was mutated in the following manner, and the mutated control variable was designated +.>；

(34)；

In the method, in the process of the invention,represents->Iterating for the second time; />Representing a mutation operator; />Representing three different individual serial numbers;indicate->Sequence number after repeated iterative mutation is->Is a control variable of (2); />Indicate->The sequence number of the iteration isIs used for controlling the initial control variable of the (a); / >Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a); />Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a);

next, the crossover is performed, and the variables after crossover are recorded asFor each->The method comprises the following steps:

(35)；

wherein,a random number representing 0 to 1;/>representing a crossover operator; />Indicate->Serial number->A partial fragment of the control variable of (a); />Indicate->Sequence number->A partial fragment of the control variable of (a); />Indicate->The number of times of iteration is->Partial fragments of the initial control variable of (a);

variable after variation crossoverFor->Each of the sub-variables in (1) is subjected to an objective function operation and then +.>Middle->Sub-variables->And->First->Sub-variables->The function values of (2) are compared if->Function value ratio->Is larger than->As->Iterative +.>。

The invention has the beneficial technical effects that: by simplifying the low-permeability fractured reservoir and constructing a dual medium model of a communication pipeline, the reservoir can be rapidly modeled on the premise of ensuring the simulation precision, and the invention also provides automatic history fitting and automatic production optimization model establishment, so that decisions are more efficient, the geological features and fracture properties of the low-permeability reservoir can be better known, the efficiency is improved, the cost is reduced, and the economic benefit of oilfield development is increased.

Drawings

FIG. 1 is a flow chart of a method of the present invention for rapid simulation and optimization of low permeability fractured reservoirs.

FIG. 2 is a graph comparing daily oil production curves of the first iteration in history matching according to the embodiment of the invention.

FIG. 3 is a graph comparing daily oil production curves of the second iteration in the history fit of the embodiment of the present invention.

FIG. 4 is a graph comparing daily oil production curves of the third iteration in the history matching according to the embodiment of the present invention.

FIG. 5 is a graph showing the comparison of the daily oil production curve of the production well 1 after history-fitting with the daily oil production curve of the numerical simulator according to the embodiment of the present invention.

FIG. 6 is a production optimization regime diagram of four production wells according to an embodiment of the present invention.

FIG. 7 is a graph comparing the results of the production optimization of the examples of the present invention.

Detailed Description

The invention is described in further detail below with reference to the attached drawings and detailed description:

the invention provides a rapid simulation and optimization method suitable for a low-permeability fractured reservoir, which mainly comprises the following steps:

dividing an oil reservoir into communicating pipelines connected with wells through a limiting condition of well pattern subdivision based on relative coordinates of well points in an oil reservoir system, uniformly dividing grids in the pipelines, abstracting a divided grid matrix into matrix grids and fracture grids, wherein each grid has own geological parameters, considering non-Darcy seepage of a low-permeability fractured oil reservoir in the matrix system, considering stress sensitivity effect of the rock in the fracture system, and combining a material balance equation, thereby constructing a dual medium model, and carrying out full implicit solution on parameters such as pressure, water saturation and the like of each grid;

The invention utilizes an ES-MDA history fitting algorithm and combines constraint conditions to establish an automatic history fitting mathematical model. And (3) constructing a history fit objective function, and adjusting model parameters to enable a model simulation result to be matched with actual observed oil reservoir production dynamics to obtain an oil reservoir model capable of guiding subsequent production optimization, wherein in the history fit process, the model is ensured to be matched with actual geological knowledge in order to reduce the randomness of the model, the constraint is added to the model parameters, and an automatic history fit mathematical model is established. Solving by using an ES-MDA history fitting algorithm, and updating the model parameter vector until the whole iteration process is completed or the fitting precision meets the engineering requirement;

the invention takes the economic net present value as an objective function of production optimization, establishes an oil reservoir injection and production optimization model by restraining the geological reserves of the oil reservoir, the pressure of each well and the upper and lower limits of the injection and production liquid, and obtains a production regulation system for maximizing the economic net present value of the oil reservoir through selection, intersection and variation by adopting a differential evolution algorithm.

As shown in FIG. 1, the method for rapidly simulating and optimizing the low permeability fractured reservoir specifically comprises the following steps:

Step 1, abstracting an oil reservoir system into a matrix grid and a crack grid, and constructing a dual medium model according to seepage characteristics of the matrix and the crack.

Based on relative well point coordinates in an oil reservoir system, dividing the oil reservoir into communicating pipelines connected with wells through a limiting condition of well pattern subdivision, uniformly dividing grids in the pipelines, abstracting a divided grid matrix into matrix grids and fracture grids, wherein each grid has own geological parameters, considering non-Darcy seepage of a low-permeability fractured oil reservoir in the matrix system, considering stress sensitivity effect of rock in the fracture system, and combining a material balance equation, thereby constructing a dual medium model, and carrying out full implicit solution on parameters such as pressure, water saturation and the like of each grid.

The specific process of the step 1 is as follows:

step 1.1, dispersing a reservoir into a series of nodes according to well point distribution in an oil reservoir system, and connecting all the nodes pairwise to form an initial well pattern model; the initial well pattern model is a well pattern model obtained by connecting all well points in an oil reservoir in pairs;

step 1.2, fluid exchange is not carried out between two wells corresponding to all connection, so that an initial well pattern model needs to be optimally screened, and therefore, the initial well pattern model is screened by introducing judging conditions, a series of connected inter-well connecting pipelines meeting the conditions are generated, the reservoir volume is divided into each connecting pipeline, and the volume of each connecting pipeline is multiplied by the effective thickness by the area of a quadrangle formed by the well points of the connecting pipelines and the centers of gravity of the triangles communicated with the two sides.

Wherein, the introduced judging conditions comprise a distance condition, an angle condition and a maximum connection point number condition;

each well has an own action range, and the water injection well cannot be applied to a production well far away from the production well, so that the distance screening needs to be carried out on an initial well pattern model, and the threshold value of the distance condition is set as follows: typically taking a value of 2 to 3 times the minimum well spacing, wells exceeding this threshold are considered unconnected.

The angle condition is specifically set as follows: in the triangle formed by any three connecting pipelines, the minimum internal angle cannot be smaller than the set minimum angle, otherwise, the water injection well flows to the other well through the middle well instead of directly taking effect according to the fact that the three wells are in the same path.

The maximum connection point number condition is specifically set as follows: each well point allows for a number of well points of 4-8 connections.

Step 1.3, dividing grids on each connecting pipeline into matrix grids and crack grids, wherein all the matrix grids form a matrix system, and all the crack grids form a crack system; each grid has its own parameters of permeability, porosity, volume, pressure, and water saturation, wherein the matrix grid is the primary reservoir unit, the fracture grid is the primary flow unit, the inter-grid fluid flow capacity is determined by the inter-grid conductivity, fluid flows are all between the matrix and the matrix, between the fracture and between the matrix and the fracture, and the fluid flow between the matrix and the fracture depends on the pressures of the two systems of the matrix system and the fracture system, thereby converting the pipeline connection of the initial well pattern into the connection between the grids in the dual medium pattern.

The dividing grids are divided into 5-10 grids in each connecting pipeline, and the grid volume where well points are located can be specified (such as) If a well point grid volume is specified, the grid volume default value is calculated as follows:

(1)；

wherein,is +.>The mesh volume where the mesh is located; />Is the total connection volume; />The number of the connecting pipelines;the number of grids in each connecting pipeline; />Is the number of well points. The meaning of equation (1) is that the volume of the well point grid = total reservoir volume/(grid number + well point number).

After well point grid volume is determined, the calculation formula of the grid volume in the connecting pipeline is as follows:

(2)；

wherein,、/>different well points positioned at two ends of the connecting pipeline; />Indicating the%>A grid; />Representation->、/>The connection pipe of the well point is +.>The volume of the mesh; />Is->、/>The volume of connecting pipes between well points; />Andwell points->And (4) well point->The mesh volume where the mesh is located; />Is +.>The number of pipes connected; />Is +.>Number of pipes connected.

In the matrix grid and the crack grid, the crack grid volumeMatrix mesh volume->The calculation is as follows:

(3)；

(4)；

(5)；

(6)；

wherein,representing the number of crack spacing units in the coarse grid; />Representing the volume fraction of microcracks in the whole unit, 1; / >The distance m between micro cracks in the oil reservoir; />And m is the opening degree of micro cracks in the oil reservoir system.

Each grid has its own permeability, porosity, pressure and water saturation, which are obtained by the difference between the values of the parameters of two well points, such as the porosity of two well pointsAnd->Division between->Grid, then->The porosity of the individual cells is->The remaining parameters are pushed in this way.

The method for calculating the conductivity between grids comprises the following steps:

(7)；

wherein,for grid->And grid->An inter-conductivity; />For grid->And grid->Contact area, m ² ；/>For grid->Center to grid->And grid->Is a distance of the contact surface; />For grid->Center to grid->And gridIs a distance of the contact surface; />For grid->Pseudo-permeability->And grid->Pseudo-permeability->The pseudo-permeability of which is mentioned in step 1.4.

The fluid flow between the matrix and the crack depends on the pressure of the matrix system and the crack system, and the fluid channeling between the matrix and the crack is calculated as follows:

(8)；/>

wherein,the flow rate is kg/s of channeling among the matrix cracks; />Is of a shape factor of 1/m ² ；/>Is fluid density, kg/m ³ ；/>For the permeability of the matrix network, +. >；/>For fluid viscosity->；/>For matrix system pressure>；/>For fracture system pressure +.>；/>For matrix mesh volume, m ³ 。

Step 1.4, considering non-Darcy seepage in a matrix system and considering stress sensitivity effect of rock in a fracture system for a low-permeability fractured reservoir;

the equation of motion of the matrix system, taking into account the fidaxs percolation, is expressed as follows:

(9)；

wherein,for matrix seepage velocity, +.>；/>For the permeability of the matrix network, +.>；/>In terms of the viscosity of the fluid,；/>for flow pressure difference>；/>And->Two different nonlinear percolation coefficients.

After the stress sensitivity effect of the rock is considered, the motion equation of the fracture system is as follows:

(10)；

wherein,for crack seepage velocity, < >>；/>For permeability of the fracture lattice, +.>；/>Is a constant value, and is used for the treatment of the skin,；/>for fluid viscosity->；/>Is stress sensitivity coefficient>；/>For the fracture lattice pressure to be the same,；/>for the initial pressure of the fracture grid +.>；/>For flow pressure difference>。

The formula (9) and the formula (10) are expressed in a unified form as the following equation of motion:

(11)；

wherein,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>Is->Is>；/>Is->Is of pseudo-permeability,/->；/>Is->Fluid viscosity of>；/>For the flow differential pressure to be the same,。

considering non-darcy seepage and stress-sensitive effects, the permeability can be considered as a function of pressure, and then the pseudo-permeability of the matrix and fracture is expressed as follows:

(12)；

(13)；

Wherein,and->Respectively representing a matrix and a crack; />For grid->And grid->Inter-consideration of matrix after Fodaxil seepage->Is of pseudo-permeability,/->；/>Respectively is grid->And grid->Stress-sensitive post-crack->Pseudo permeability of the interior->；/>For grid->And grid->Interstitial matrix->Is>；/>Representing grid->And grid->Crack between->Is>；/>Viscosity of fluid, ++>；/>Is a matrix inner grid->And grid->Pressure gradient of two nodes->；/>And->Two different nonlinear seepage coefficients, the values of which are determined by experiments; />Is stress sensitivity coefficient>；/>For grid->And grid->Current pressure of interstitial cracks>；/>For grid->And grid->Initial stress of interstitial cracks>。

By this process, the equations of motion of the matrix system and fracture system can be made to take a consistent form, which will facilitate further mathematical model construction and programming work.

And 1.5, constructing a low-permeability fractured reservoir dual medium model according to a formula (11) by combining seepage rules of the matrix and fluid in the fracture.

The dual media model constructed should conform to the following assumptions: the seepage process is isothermal, and the seepage behavior in the oil reservoir is not affected by temperature change; the model considers the simultaneous existence of an oil phase and a water phase, and the two phases are not miscible and all follow the same permeation rule; the influence of capillary force can be considered, and the influence of gravity is not considered temporarily; both the fluids and the rock in the formation are compressible.

The crack system and the matrix system are two abstract grid systems which exist in space at the same time and are mutually overlapped, fluid seepage in the matrix system is nonlinear seepage, fluid flow in the crack system is Darcy seepage, and stress sensitivity effect is considered.

The dual medium model comprises a continuity equation of a matrix system and a continuity equation of a crack system;

wherein the continuity equation of the matrix system is as follows:

(14)；

(15)；

wherein,for oil phase density, < >>；/>Porosity of the matrix, 1; />Is the saturation of the oil phase in the matrix, 1; />For matrix grid volume, +.>；/>Time, s; />For the rate of percolation of the oil phase in the matrix, +.>；/>Representing the divergence; />Is a form factor->；/>For pseudo-permeability in the matrix, +.>；/>Viscosity of oil phase>；For the pressure of the oil phase in the matrix>；/>Is the oil phase pressure in the crack, < >>；/>The injection and extraction rate of the oil phase in the matrix is positive, the injection and extraction are negative, and the +.>；/>For the density of the aqueous phase->；/>Saturation of the aqueous phase in the matrix, 1; />For the rate of percolation of the aqueous phase in the matrix, +.>；/>For the aqueous phase viscosity>；/>For the pressure of the aqueous phase in the matrix>；/>Water phase pressure in crack->；/>The injection and extraction rate of the water phase in the matrix is positive, the injection and extraction rate is negative, and the +.>；

The continuity equation of the fracture system is as follows:

(16)；

(17)；

Wherein,porosity of the crack, 1; />Saturation of oil phase in cracks, 1; />Saturation of water phase in cracks, 1; />For the fracture mesh volume>；/>For the rate of penetration of the oil phase in the cracks, < > and>；/>for the rate of penetration of the aqueous phase in the cracks, < > water>；/>The injection and production rate of the oil phase in the crack is positive, the production is negative,；/>the injection and extraction rate of the water phase in the crack is positive, the injection and extraction rate is negative, and the +.>；

Order theRepresenting the internal phase fluid +.>Variable value of volume flow,/, for>Representing the internal phase fluid +.>The variation of volume flow, i.e.，/>Wherein->Is a phase fluid; />Representing->Is oil phase->Or water phase->；/>Is the fluid in the matrix->Injection and production rate of (2) is positive, injection and production is negative,/->；/>Is fluid->Density of->；/>Is the fluid in the matrix->Pressure of->；/>Is the intra-fracture phase fluid->Pressure of->；/>Is the intra-fracture phase fluid->Injection and production rate of (2) is positive, injection and production is negative,/->The method comprises the steps of carrying out a first treatment on the surface of the At this time, the continuity equation of the matrix system may become:

(18)；

(19)；

wherein,is the change value of the volume flow of the oil phase fluid in unit volume of unit time in the matrix grid,；/>is the change value of the volume flow rate of the aqueous phase fluid in unit volume of unit time in the matrix grid, +. >。

The continuity equation of the fracture system can be changed as:

(20)；

(21)；/>

wherein,is the change value of the volume flow of the oil phase fluid in unit volume of unit time in the fracture grid,；/>is the variation value of the volume flow of the aqueous phase fluid in the unit time unit volume of the crack grid, +.>。

And 1.6, performing linear interpolation approximation on the established continuity equation, wherein a time item of the continuity equation adopts backward first-order difference, a space item adopts a central difference format, a numerical discrete result is written into a residual form, and parameters such as pressure and water saturation in a crack and matrix system are completely and implicitly solved by using a Newton-Raphson method (Newton-Lag Pei Sen method). The specific process is as follows:

step 1.6.1, defining fluidity coefficients as follows:

(22)；

wherein,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>For the fluidity of a certain fluid->；/>For the relative permeability of a phase fluid, +.>；/>For the viscosity of a fluid of a certain phase->。

Step 1.6.2, because the seepage characteristics of the matrix system and the fracture system and the cross flow between the matrix system and the fracture system can be comprehensively considered in the conductivity among grids, the matrix grids and the fracture grids have no essential difference except the parameters, so that the two sets of systems can be approximated by linear interpolation in a uniform format, and the complexity of subsequent programming can be greatly reduced.

The linear interpolation approximation uses the following numerical discrete equation:

(23)；

in the method, in the process of the invention,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>For the density of a certain phase of fluid->；/>Is porosity; />For saturation of a fluid of a certain phaseDegree of sum; />For the mesh volume +.>；/>The last moment in the simulation iteration process; />The current time in the simulation iteration process is; />For the time step +.>；/>For grid->In the model, taking matrix grids as an example, adjacent units only have grids in the same pipeline beside the matrix grids and corresponding crack grids; subscript->Representing two adjacent meshes->And->An average value of a parameter on the interface;for grid->And grid->An inter-conductivity; />For a certain fluid>Time grid->Pressure of->；/>For a certain fluid>Time grid->Pressure of->；/>For grid->Fluid in certain phase at->Injection and production rate at moment, positive injection, negative production, and +.>。

Step 1.6.3, writing the numerical discrete equation into the residual form as follows:

(24)；

in the method, in the process of the invention,equation residual error; />，/>For the total number of grids, this time step option is +.>Therefore, the method is a fully implicit solution, and the pressure and the water saturation in the fracture and the matrix can be solved by using a Newton-Paphson iteration solution method.

The Newton-Paphson iteration solving method is that a new point which is closer to a function zero point is obtained through continuous iteration calculation from an initial point, the new point is obtained through Taylor expansion calculation at a current point, and then the new point is used as the initial point of the next iteration to continue iteration calculation until a solution meeting the requirement is found.

And 2, establishing an automatic history fit mathematical model by utilizing an ES-MDA (ensemble smoothing multiple data assimilation) history fit algorithm and combining constraint conditions.

And (3) constructing a history fit objective function, and adjusting model parameters to enable a model simulation result to be matched with actual observed oil reservoir production dynamics to obtain an oil reservoir model capable of guiding subsequent production optimization, wherein in the history fit process, the model is ensured to be matched with actual geological knowledge in order to reduce the randomness of the model, the constraint is added to the model parameters, and an automatic history fit mathematical model is established. And solving by using an ES-MDA history fitting algorithm, and updating the model parameter vector until the whole iteration process is completed or the fitting precision meets the engineering requirement.

The step 2 specifically comprises the following steps:

step 2.1, the construction of the objective function of the automatic history fit mathematical model is essentially to minimize the error between the model calculation result and the actual observed data, and thus can be expressed as:

(25)；

In the method, in the process of the invention,the objective function value is the difference between the calculated result and the observed data; />Representing the model calculation result, dimension is +.>Wherein->Representing the number of observations; />Representing actual observed data, the dimension is；/>Covariance matrix representing observation error, dimension is +.>。

And 2.2, the automatic history fit mathematical model greatly reduces the characteristic parameters of the oil reservoir model through a simplified method, so that in the history fit process, only the volume of each pipeline, the conductivity between grids and the parameters of the power law phase permeability curve need to be adjusted.

The parameters of the power law phase permeation curve are adjusted by the following power law phase permeation model:

(26)；

wherein,

(27)；/>

in the method, in the process of the invention,and->Respectively representing the relative permeability of the oil phase and the water phase; />Representing dimensionless saturation; />Anddetermining the bending degree of the relative permeability curves of the water phase and the oil phase respectively; />Representing irreducible water saturation; />Representing residual oil saturation; />Representing the relative permeability value of the aqueous phase at a water saturation of irreducible water saturation; />Is the water saturation value. The variation of these parameters and values together describe the variation of the permeation characteristics of the multiphase fluid in the porous medium.

Step 2.3, parameters for history fitting adjustment The following are provided:

(28)；

wherein,for grid->And grid->The inter-conductivity comprises a matrix and a matrix grid, a crack and a crack grid and a matrix and a crack grid; />Is->、/>Matrix-conducting volume of the communication conduit of the well point, < ->Is->、/>Fracture conduction volume of the communication tubing at the well point; />、/>、/>Is->、/>Different parameters of the communication pipeline phase permeation curve of the well point; />And->Two different nonlinear seepage coefficients;

when all pairs of wells share a set of phase permeability data, the parameters of the history fitThe expression is as follows:

(29)；

wherein,、/>、/>three parameters of the permeability curve are common among all pairs of wells.

Step 2.4, in the history fitting process, adding constraint to model parameters; constraints of the automatic history fit mathematical model are as follows:

and 2.4, in the history fitting process, in order to reduce the randomness of the model and ensure that the model accords with the actual geological knowledge, the constraint needs to be added to the model parameters, so that an automatic history fitting mathematical model is established.

According to the definition of conductivity and volume, the value of conductivity and volume is non-negative, the volume of each pipeline is smaller than the total conduction volume, the sum of all pipeline volumes is equal to the total conduction volume, the parameters in the phase permeation model are subject to the following constraint, the characteristic parameters of the nonlinear seepage curve are larger than 0, and the comprehensive compression coefficient is generally 0.00001-0.1, so the constraint condition of the automatic history fit mathematical model is as follows:

(30)；

Wherein,for grid->And grid->Is defined by the volume of matrix; />Is the total volume of the matrix system; />For grid->And grid->Is defined by the fracture volume of (2); />Is the total volume of the fracture system; />Is the total well count.

And 3, taking the economic net present value as an objective function of production optimization, establishing an oil reservoir injection and production optimization model by restraining the geological reserves of the oil reservoirs, the pressure of each well and the upper and lower limits of the injection and production liquid, and obtaining a production regulation system for maximizing the economic net present value of the oil reservoirs by adopting a differential evolution algorithm.

The step 3 specifically comprises the following steps:

and 3.1, adopting an economic net present value as an objective function of production optimization, wherein the mathematical expression is as follows:

(31)；

wherein,is an economic net present value; />Is a control variable; />For time step number, ++>Representing a total number of time steps; />Indicating production well number, +.>Representing the number of production wells; />Indicating the serial number of the water injection well>Representing the number of water injection wells; />Price per ton of crude oil; />For producing wells->At->Oil extraction of the time step; />Cost per ton of sewage to be treated; />Representing a production well->At->The water yield of the time step; />Representing water injection costs per ton of water; />Indicating water injection well->At->The water injection amount of the time step;

And 3.2, in the production optimization process, controlling the oil reservoir parameters by adopting boundary constraint conditions, wherein the constraint conditions of the oil reservoir parameters are as follows:

(33)；

wherein,is a water injection well->At->The water injection amount of the time step; />For producing wells->At->The liquid extraction amount of the time step;is a filling and production well->At->Time step injection of liquid production volume->The injection well comprises a water injection well and a production well; />Is a filling and production well->At->Time step injection of liquid production volume->Upper limit of (2); />Is the total geological reserve; />And->Injection well and production well respectively>Control pressure->Lower and upper limits of (2);

and 3.3, obtaining a production regulation system for maximizing the economic net present value of the oil reservoir by adopting a differential evolution algorithm.

The differential evolution algorithm, i.e. first initializing a series of control variablesThe series of control variables are mutated in the following manner,the control variable after mutation is marked +.>；

(34)；

In the method, in the process of the invention,represents->Iterating for the second time; />Representing a mutation operator, wherein any number between 0 and 2 can be adopted; />Representing three different individual serial numbers; />Indicate->Sequence number after repeated iterative mutation is->Is a control variable of (2); />Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a); />Indicate->Sub-stackThe code number of the generation is- >Is used for controlling the initial control variable of the (a); />Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a); that is to say that three individuals are selected from the initial control variables for the above operation, resulting in the sequence number +.>Is controlled by the control variable of (a).

Next, the crossover is performed, and the variables after crossover are recorded asFor each->The method comprises the following steps:

(35)；

wherein,a random number representing 0 to 1; />Representing a crossover operator; />Indicate->Serial number->A partial fragment of the control variable of (a); />Indicate->Sequence number->A partial fragment of the control variable of (a); />Indicate->The number of times of iteration is->Partial fragments of the initial control variable of (a); if the random number is larger than the crossover operator, the crossover variable is mutated +.>If the random number is smaller than the crossover operator, the original +.>An initial population of iterations.

Variable after variation crossoverFor->Each of the sub-variables in (1) is subjected to an objective function operation and then +.>Middle->Sub-variables->And->First->Sub-variables->The function values of (2) are compared if->Function value ratio->Is larger than->As->Iterative +.>。

In order to demonstrate the feasibility and superiority of the invention, the following examples are given.

Firstly, establishing a dual medium model by using a finite volume method oil reservoir numerical simulator, wherein the grid scale of the dual medium model is as follows: 25 x 1, model size: 100m x 10m, and the porosity of the reservoir matrix is 0.2; crack porosity of 0.02, micro-crack opening of 1×10 ^-5 m is that the crack spacing is 5m, the initial water saturation of the oil deposit is 0.2, the oil-water viscosity is respectively 4 mPa s and 1 mPa s, 1 water injection well is arranged in the oil deposit, 4 production wells are arranged in the center of the oil deposit, four production wells are arranged at four corner points of the oil deposit, the oil deposit numerical simulator is adopted for simulating production for 1500 days, and the water injection rate of the water injection well is 0.8m ³ And/d, the liquid production rates of the production wells are all 0.2m ³ And/d, overall block mining balance.

Then, the method is used, the same reservoir parameter modeling is used, the single forward running time is 15s, compared with the time consumption of a limited volume reservoir numerical simulator, the calculation speed is increased by approximately 5 times, the history fitting is carried out by adopting an ESMDA algorithm, the number of groups is 100, the iteration is 5 times, the oil production speed change of a production well 1 in the history fitting is shown in fig. 2, 3 and 4, the fig. 2, 3 and 4 are respectively oil production speed comparison diagrams in the first iteration, the second iteration and the third iteration of the history fitting, the abscissa is time, the ordinate is daily oil production speed, the gray fine curve is a daily oil production curve obtained by each control variable and reflects the daily oil production speed change of each control variable, and the black coarse curve is a daily oil production curve of the reservoir numerical simulator and reflects the daily oil production speed change of the numerical simulator. Fig. 5 shows a comparison curve of the daily oil production curve of the production well 1 obtained after history fitting and the oil production speed calculated by the conventional reservoir numerical simulator, and it can be seen that the daily oil production curve of the production well 1 after history fitting substantially matches the simulated real daily oil production curve of the numerical simulator. Finally, the production optimization is carried out by using the model after history fitting, and the production optimization system of each production well comprises initial regulation and control and first, second and third regulation, and the liquid production speed after regulation is shown in figure 6. In this embodiment, three time steps are optimized, each time step is 100 days, the total time is 300 days, the optimization result comprises block accumulated oil and block water content, specific data are shown in fig. 7, the optimized accumulated oil is obviously improved, and the water content is relatively stable, so that the reliability of the method in low-permeability fractured reservoir simulation and optimization calculation is verified.

It should be understood that the above description is not intended to limit the invention to the particular embodiments disclosed, but to limit the invention to the particular embodiments disclosed, and that the invention is not limited to the particular embodiments disclosed, but is intended to cover modifications, adaptations, additions and alternatives falling within the spirit and scope of the invention.

Claims (10)

1. The rapid simulation and optimization method suitable for the low-permeability fractured reservoir is characterized by comprising the following steps of:

step 1, abstracting an oil reservoir system into a matrix grid and a crack grid, and constructing a dual medium model according to seepage characteristics of the matrix and the crack;

step 2, establishing an automatic history fit mathematical model by utilizing an ES-MDA history fit algorithm and combining constraint conditions;

and 3, taking the economic net present value as an objective function of production optimization, establishing an oil reservoir injection and production optimization model by restraining the geological reserves of the oil reservoirs, the pressure of each well and the upper and lower limits of the injection and production liquid, and obtaining a production regulation system for maximizing the economic net present value of the oil reservoirs by adopting a differential evolution algorithm.

2. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 1, wherein the specific process of step 1 is as follows:

step 1.1, dispersing a reservoir into a series of nodes according to well point distribution in an oil reservoir system, and connecting all the nodes pairwise to form an initial well pattern model;

Step 1.2, screening an initial well pattern model by introducing judging conditions, generating well connecting pipelines connected in pairs, and dividing the volume of a reservoir into each connecting pipeline;

step 1.3, dividing grids on each connecting pipeline into matrix grids and crack grids, wherein all the matrix grids form a matrix system, and all the crack grids form a crack system; the inter-grid fluid flow capacity is determined by the magnitude of the inter-grid conductivity;

step 1.4, considering non-Darcy seepage in a matrix system and considering stress sensitivity effect of rock in a fracture system;

step 1.5, constructing a low-permeability fractured reservoir dual medium model by combining seepage rules of the matrix and fluid in the fracture; the dual medium model comprises a continuity equation of a matrix system and a continuity equation of a crack system;

and 1.6, performing linear interpolation approximation on the continuity equation, and fully implicitly solving parameters in the crack and matrix system by using a Newton-Raphson method.

3. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 2, wherein in step 1.2, the determination conditions include a distance condition, an angle condition and a maximum number of connection points condition; wherein the threshold value of the distance condition is 3 times of the minimum well distance, and wells exceeding the threshold value are regarded as unconnected; the angle conditions are specifically set as follows: in a triangle formed by any three connecting pipelines, the minimum internal angle cannot be smaller than the set minimum angle; the maximum connection point number condition is specifically set as follows: each well point allows for a number of well points of 4-8 connections.

4. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 3, wherein in the step 1.3, the division of grids is divided into 5-10 grids in each connecting pipeline, and when the well point grid volume is specified, the calculation formula of the well point grid volume default value is as follows:

(1)；

wherein,is +.>The mesh volume where the mesh is located; />Is the total connection volume; />The number of the connecting pipelines; />The number of grids in each connecting pipeline; />The number of well points;

after well point grid volume is determined, the calculation formula of the grid volume in the connecting pipeline is as follows:

(2)；

wherein,、/>different well points positioned at two ends of the connecting pipeline; />Indicating the%>A grid; />Representation of、/>The connection pipe of the well point is +.>The volume of the mesh; />Is->、/>The volume of connecting pipes between well points; />And->Well points->And (4) well point->The mesh volume where the mesh is located; />Is +.>The number of pipes connected; />Is +.>The number of pipes connected;

matrix grid and fracture grid, fracture grid volumeMatrix mesh volume->The calculation is as follows:

(3)；

(4)；

(5)；

(6)；

wherein,representing the number of crack spacing units in the coarse grid; />Representing the volume fraction of microcracks in the whole unit; / >The micro-crack spacing in the oil reservoir; />The opening degree of micro cracks in an oil reservoir system is set;

the calculation formula of the conductivity between grids is as follows:

(7)；

wherein,for grid->And grid->An inter-conductivity; />For grid->And grid->Is a contact area of (a); />Is a gridCenter to grid->And grid->Is a distance of the contact surface; />For grid->Center to grid->And grid->Is a distance of the contact surface; />For grid->Pseudo-permeability->And grid->Pseudo-permeability->Average value of (2);

fluid flow exists between the matrix and the fracture, and the fluid channeling is calculated as follows:

(8)；

wherein,the cross flow between the matrix cracks; />Is a form factor; />Is the fluid density; />Permeability for matrix network; />Is the viscosity of the fluid; />Is the matrix system pressure; />Is fracture system pressure; />Is the matrix grid volume.

5. The method according to claim 4, wherein in step 1.4, the equation of motion of the matrix system after considering the feddalix seepage is:

(9)；

wherein,is the matrix seepage velocity; />Permeability for matrix network; />Is the viscosity of the fluid; />Is the flow differential pressure;and->Two different nonlinear seepage coefficients;

After the stress sensitivity effect of the rock is considered, the motion equation of the fracture system is as follows:

(10)；

wherein,is the crack seepage velocity; />Permeability for the fracture grid; />Is a constant; />Is the viscosity of the fluid; />Is the stress sensitivity coefficient; />The grid pressure is the crack; />Initial pressure for the fracture grid; />Is the flow differential pressure;

the formula (9) and the formula (10) are expressed in a unified form as the following equation of motion:

(11)；

wherein,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>Is->Is not less than the percolation rate; />Is->Is to be determined; />Is->Is a fluid viscosity of (a); />Is the flow differential pressure;

the pseudo-permeability of the matrix and fracture, considering non-darcy seepage and stress-sensitive effects, is expressed as follows:

(12)；

(13)；

wherein,and->Respectively representing a matrix and a crack; />For grid->And grid->Inter-consideration of matrix after Fodaxil seepage->Is to be determined; />Respectively is grid->And grid->Stress-sensitive post-crack->The pseudo permeability of the inner;for grid->And grid->Interstitial matrix->Is a natural gas; />Representing grid->And grid->Crack between->Is a natural gas; />Is the viscosity of the fluid; />Is a matrix inner grid->And grid->A pressure gradient of the two nodes; />And->Two different nonlinear seepage coefficients; / >Is the stress sensitivity coefficient; />For grid->And grid->The current pressure of the inter-fracture; />For grid->And grid->Initial pressure of the inter-fissures.

6. The method according to claim 5, wherein in step 1.5, the continuity equation of the matrix system is:

(18)；

(19)；

wherein,is oil phase density; />Porosity of the matrix; />Saturation of the oil phase in the matrix; />Is a matrix grid volume; />Is a time step number; />Representing the divergence; />The seepage velocity of the oil phase in the matrix; />The change value of the volume flow of the oil phase fluid in unit volume of unit time in the matrix grid; />Is the density of the water phase; />Saturation of the aqueous phase within the matrix; />The seepage velocity of the water phase in the matrix; />The change value of the volume flow of the aqueous phase fluid in unit volume of unit time in the matrix grid;

the continuity equation of the fracture system is:

(20)；

(21)；

wherein,porosity of the fracture; />Saturation of the oil phase in the fracture; />Is the seepage velocity of the oil phase in the cracks; />The change value of the volume flow of the oil phase fluid in unit volume of unit time in the fracture grid is obtained; />Saturation of the aqueous phase in the fracture; / >Is the seepage velocity of the water phase in the crack; />The change value of the volume flow of the aqueous phase fluid in the unit volume of the unit time of the fracture grid is obtained.

7. The method for rapid simulation and optimization of a low permeability fractured reservoir according to claim 6, wherein the specific process of step 1.6 is as follows:

step 1.6.1, defining fluidity coefficients as follows:

(22)；

wherein,is a certain fluid; />Representing->Is oil phase->Or water phase->；/>Is the fluidity of a certain phase fluid; />Is the relative permeability of a fluid of a certain phase; />Is the viscosity of a certain fluid;

step 1.6.2, performing linear interpolation approximation on the matrix system and the fracture system by adopting the same format, wherein the linear interpolation approximation adopts the following numerical discrete equation:

(23)；

in the method, in the process of the invention,is the density of a fluid in a certain phase; />Is porosity; />Saturation for a fluid of a certain phase; />Is a mesh volume; />The last moment in the simulation iteration process; />The current time in the simulation iteration process is; />Is the time step; />For grid->A set of adjacent cells; subscript->Representing two adjacent meshes->And->An average value of a parameter on the interface; />For grid->And grid->An inter-conductivity; />For a certain fluid>Time grid->Is a pressure of (2); / >For a certain fluid>Time grid->Is a pressure of (2); />For grid->Fluid in certain phase at->Injection and production rate at moment;

step 1.6.3, writing the numerical discrete equation into the residual form as follows:

(24)；

in the method, in the process of the invention,equation residual error;

and solving the pressure and water saturation parameters in the fracture and matrix system by adopting a Newton-Paphson iteration solving method, wherein the solving process is as follows: starting from the initial point, a new point approaching to the zero point of the function is obtained through continuous iterative computation, the new point is obtained through Taylor expansion computation at the current point, and then the new point is used as the initial point of the next iteration to continue iterative computation until a solution meeting the requirement is found.

8. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 7, wherein the specific process of step 2 is as follows:

step 2.1, constructing an objective function of an automatic history fit mathematical model, as follows:

(25)；

in the method, in the process of the invention,is the objective function value; />Representing a model calculation result; />Representing actual observed data; />A covariance matrix representing an observation error;

step 2.2, in the history fitting process of the automatic history fitting mathematical model, only the volume of each pipeline, the conductivity between grids and the parameters of the power law phase permeability curve need to be adjusted;

The parameters of the power law phase permeation curve are adjusted by the following power law phase permeation model:

(26)；

wherein,

(27)；

in the method, in the process of the invention,and->Respectively representing the relative permeability of the oil phase and the water phase; />Representing dimensionless saturation; />And->Determining the bending degree of the relative permeability curves of the water phase and the oil phase respectively; />Representing irreducible water saturation; />Representing residual oil saturation; />Representing the relative permeability value of the aqueous phase at a water saturation of irreducible water saturation; />Is saturated with waterA degree value;

step 2.3, parameters for history fitting adjustmentThe following are provided:

(28)；

wherein,for grid->And grid->The inter-conductivity comprises a matrix and a matrix grid, a crack and a crack grid and a matrix and a crack grid; />Is->、/>Matrix-conducting volume of the communication conduit of the well point, < ->Is->、/>Fracture conduction volume of the communication tubing at the well point; />、/>、/>Is->、/>Different parameters of the communication pipeline phase permeation curve of the well point; />And->Two different nonlinear seepage coefficients;

when all pairs of wells share a set of phase permeability data, the parameters of the history fitThe expression is as follows:

(29)；

wherein,、/>、/>three parameters of the phase permeability curve shared by all well pairs;

step 2.4, in the history fitting process, adding constraint to model parameters; constraints of the automatic history fit mathematical model are as follows:

(30)；

Wherein,for grid->And grid->Is defined by the volume of matrix; />Is the total volume of the matrix system; />For grid->And grid->Is defined by the fracture volume of (2); />Is the total volume of the fracture system; />Is the total well count.

9. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 8, wherein the specific process of step 3 is as follows:

and 3.1, adopting an economic net present value as an objective function of production optimization, wherein the mathematical expression is as follows:

(31)；

wherein,is an economic net present value; />Is a control variable; />For time step number, ++>Representing a total number of time steps; />Indicating production well number, +.>Representing the number of production wells; />Indicating the serial number of the water injection well>Representing the number of water injection wells; />Price per ton of crude oil; />For producing wells->At->Oil extraction of the time step; />Cost per ton of sewage to be treated; />Representing a production well->At->The water yield of the time step; />Representing water injection costs per ton of water; />Indicating water injection well->At->The water injection amount of the time step;

and 3.2, in the production optimization process, controlling the oil reservoir parameters by adopting boundary constraint conditions, wherein the constraint conditions of the oil reservoir parameters are as follows:

(33)；

wherein,is a water injection well- >At->The water injection amount of the time step; />For producing wells->At->The liquid extraction amount of the time step;is a filling and production well->At->Time step injection of liquid production volume->Lower limit of (2); />Is a filling and production well->At->Time step injection of liquid production volume->Upper limit of (2); />Is the total geological reserve; />And->Injection well and production well respectively>Control pressure->Lower and upper limits of (2);

and 3.3, obtaining a production regulation system for maximizing the economic net present value of the oil reservoir by adopting a differential evolution algorithm.

10. The method for rapid simulation and optimization of low permeability fractured reservoirs according to claim 9, wherein the specific process of step 3.3 is: first, a series of control variables are initializedThe series of control variables was mutated in the following manner, and the mutated control variable was designated +.>；

(34)；

In the method, in the process of the invention,represents->Iterating for the second time; />Representing a mutation operator; />Representing three different individual serial numbers;indicate->Sequence number after repeated iterative mutation is->Is a control variable of (2); />Indicate->The sequence number of the iteration isIs used for controlling the initial control variable of the (a); />Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a); />Indicate->The sequence number of the iteration is->Is used for controlling the initial control variable of the (a);

next, the crossover is performed, and the variables after crossover are recorded as For each->The method comprises the following steps:

(35)；

wherein,a random number representing 0 to 1; />Representing a crossover operator; />Indicate->Serial number->A partial fragment of the control variable of (a); />Indicate->Sequence number->A partial fragment of the control variable of (a); />Indicate->The number of times of iteration is->Part of the initial control variable of (a)Dividing into segments;

variable after variation crossoverFor->Each of the sub-variables in (1) is subjected to an objective function operation and then +.>Middle->Sub-variables->And->First->Sub-variables->The function values of (2) are compared if->Function value ratio->Is larger than->As->Iterative +.>。