CN114896914A

CN114896914A - Unconventional reservoir high-precision intelligent fracturing regulation and control method and device

Info

Publication number: CN114896914A
Application number: CN202210661184.0A
Authority: CN
Inventors: 曾凡辉; 张宇; 胡大淦; 郭建春; 刘学伟
Original assignee: Southwest Petroleum University
Current assignee: Southwest Petroleum University
Priority date: 2022-06-13
Filing date: 2022-06-13
Publication date: 2022-08-12

Abstract

The invention discloses a method and a device for high-precision intelligent fracturing regulation and control of unconventional reservoirs, wherein the method comprises the following steps: collecting basic parameters required by calculation; establishing a perforation cluster cracking and expansion fluid-solid coupling model by utilizing a fluid-solid coupling control equation of rock solid deformation, fluid flow in the seam and matrix pressure diffusion; coupling dynamic flow distribution, initial fracture perforation cluster extension and fracture induced stress in the process of multi-fracture competitive fracture initiation and expansion, and establishing a planar multi-cluster fracture initiation-expansion coupling model; establishing a non-planar multi-cluster crack transformation volume model; based on a non-planar crack initiation and propagation mechanism model, adopting DOE design to establish a sample database, utilizing a machine learning model to train, and establishing a multi-input multi-output intelligent agent model for initiating the simulation of the mechanism model; the method comprises the following steps of taking the large reconstruction volume and the low construction pressure of each cluster of cracks matched with reservoir geological characteristics as a multi-objective function; and comprehensively solving the optimal solution of the whole objective function by using a genetic algorithm.

Description

Unconventional reservoir high-precision intelligent fracturing regulation and control method and device

Technical Field

The invention relates to the field of petroleum and natural gas engineering, in particular to a method and a device for high-precision intelligent fracturing regulation and control of an unconventional reservoir.

Background

The permeability of an unconventional reservoir is extremely low, and the formation of a balanced and expanded 'artificially modified gas reservoir' high permeability zone in a 'sweet spot region' through horizontal well staged multi-cluster fracturing is the key for realizing the efficient development of unconventional oil and gas. Due to the fact that stress anisotropy and permeability anisotropy of unconventional oil and gas reservoirs cause different degrees of difficulty in crack initiation and expansion of perforation clusters, uneven expansion of the perforation clusters is caused; in addition, the irregular reservoir multi-crack non-equilibrium expansion regulation and control relates to the multi-target optimization problems of the transformation volume of each cluster, the length of each cluster of cracks, the construction pressure and the like, and the conventional multi-target optimization algorithm for carrying out sensitive factor optimization analysis based on a multi-cluster perforation simultaneous initiation expansion mechanism model has large calculated amount and is difficult to realize global optimization. Under the condition of a heterogeneous reservoir, the flow of each cluster is different in modification volume, and how to select a perforation position and adjust perforation parameters and construction parameters under the condition of a certain construction displacement enables the flow of each cluster to play the maximum role of the cluster, namely the modification volume is maximized, and the non-equilibrium expansion regulation and control research of the fracturing fracture of the unconventional reservoir needs to be carried out by combining a multi-target intelligent algorithm.

Therefore, the research on multi-cluster fracture initiation-expansion coupling and multi-target fracturing parameter optimization is needed, and the uniform transformation of the unconventional reservoir fractured horizontal well is realized.

Disclosure of Invention

In view of this, the invention aims to provide a high-precision intelligent fracturing regulation and control method and device for an unconventional reservoir.

In order to achieve the above technical objects, the present invention provides the following technical solutions.

The method for regulating and controlling the highly-accurate intelligent fracturing of the unconventional reservoir is characterized by comprising the following steps of:

collecting basic parameters required by calculation, wherein the basic parameters required by calculation comprise reservoir parameters, mechanical parameters, fracturing parameters and shaft parameters, the reservoir parameters comprise original stratum pore pressure, reservoir permeability, reservoir thickness, comprehensive compression coefficient, reservoir porosity, ground stress and fracture zone parameters, the mechanical parameters comprise parallel bedding Young modulus, vertical bedding Young modulus, parallel bedding Poisson ratio, vertical bedding Poisson ratio, rock tensile strength, weak surface cohesive force, formation tendency, formation inclination angle and internal friction angle parameters, the fracturing parameters comprise construction displacement, fracturing fluid viscosity, cluster number, perforation depth, perforation radius, perforation length, perforation density and perforation phase angle parameters, and the shaft parameters comprise well bore radius, well bore inclination angle and shaft azimuth angle;

based on a perforation cluster crack initiation and expansion physical model, establishing a perforation cluster crack initiation and expansion fluid-solid coupling model by utilizing a fluid-solid coupling control equation of rock solid deformation, fluid flow in cracks and matrix pressure diffusion and based on a perforation cluster crack initiation criterion and an I, II type composite fracture expansion criterion;

based on a perforation cluster initiation and expansion fluid-solid coupling model, solving a rock deformation-pore pressure diffusion-fluid flow multi-field coupling variable in a seam by adopting a displacement discontinuity method, calculating a perforation cluster crack discrete point source normal displacement discontinuity quantity, a tangential displacement discontinuity quantity, a filtration loss quantity, an in-seam pressure and a perforation crack extension direction angle, and then substituting into an in-seam fluid flow equation to obtain a seam length, a seam width, an in-seam pressure and an expansion track;

on the basis of a planar perforation cluster crack initiation-expansion coupling model and multi-field coupling variable solving, reservoir physical properties, ground stress heterogeneity and perforation hole friction resistance are considered, dynamic flow distribution, initial crack perforation cluster extension and crack induced stress in the process of competing initiation and expansion of multi-cracks are coupled, the planar multi-cluster crack initiation-expansion coupling model is established, and the flow, the pressure and the multi-field coupling variable of each cluster of cracks are obtained by solving through a Newtonian Raphson method;

on the basis of a planar multi-cluster fracture initiation-expansion coupling model and multi-field coupling variable solving, taking non-planar fracture network fluid flow into consideration, dynamically distributing coupling multi-fracture and bifurcation fracture flow, and establishing a non-planar fracture initiation-expansion coupling model based on a hydraulic fracture and natural fracture interaction judgment criterion taking anisotropic induced stress into consideration; realizing non-synchronous fracture initiation-expansion simulation of the non-planar fracture; then, establishing a non-planar multi-cluster fracture reconstruction volume model based on the shear slip and tensile failure mechanical conditions of reservoir rock and a fracture network permeability model according to the change of the reservoir pore pressure field caused by the induction of the non-planar multi-cluster fractures;

based on a non-planar crack initiation and propagation mechanism model, adopting DOE design to establish a sample database, training by using a Gauss-Kriging machine learning model, and establishing a multi-input multi-output intelligent agent model for simulating the mechanism model; on the basis, the reconstruction volume of each cluster of cracks matched with the geological characteristics of the reservoir is large, the construction pressure is low, the crack length is expanded and balanced, and the geological engineering dessert is restrained; and comprehensively solving the optimal solution of the whole objective function by using an improved genetic algorithm, and establishing the unconventional reservoir fracturing multi-fracture unbalanced expansion intelligent control method.

In addition, based on the unconventional high-precision intelligent fracturing regulation and control method of reservoir that this application provided, set up the high-precision intelligent fracturing regulation and control device of unconventional reservoir that corresponds, include:

the basic parameter acquisition module is used for collecting basic parameters required by calculation;

the hydraulic fracture seam length, seam width, seam internal pressure and expansion track acquisition module is used for solving rock deformation-pore pressure diffusion-fluid flow multi-field coupling variables in the seam through a displacement discontinuity method and calculating seam length, seam width, seam internal pressure and expansion track at any time;

the system comprises a planar multi-cluster fracture flow, an internal fracture pressure and multi-field coupling variable acquisition module for each cluster of fractures, wherein the planar multi-cluster fracture flow, the internal fracture pressure and the multi-field coupling variable acquisition module is used for considering reservoir physical properties, ground stress heterogeneity and perforation friction resistance, coupling dynamic flow distribution, initial perforation cluster extension and fracture induced stress in the multi-fracture competition initiation and extension process, establishing a planar multi-cluster fracture initiation-extension coupling model, and solving by adopting a Newtonian-Raphson method to obtain planar multi-cluster fracture flow, internal fracture pressure and multi-field coupling variables of each cluster of fractures;

the non-planar multi-cluster fracture reconstruction volume prediction module is used for coupling multi-fracture and bifurcation fracture flow dynamic distribution and establishing a non-planar fracture initiation-expansion coupling model based on the interactive judgment criterion of hydraulic fractures and natural fractures considering anisotropic induced stress; according to the change of a reservoir pore pressure field caused by the induction of the non-planar multi-cluster fractures, obtaining non-planar multi-cluster fracture modification volume parameters based on reservoir rock shear slip and tensile failure mechanical conditions and a fracture network permeability model;

the intelligent regulation and control optimized perforation and construction parameter acquisition module is used for establishing a multi-input multi-output intelligent agent model to replace a mechanism model by utilizing a Gauss Criging machine learning algorithm based on a crack initiation and propagation mechanism model; and modifying a multi-target function with large volume, balanced expansion of the length of each cluster of cracks and low construction pressure matched with the geological characteristics of the reservoir to select geological engineering dessert perforation and pressure limiting constraint conditions, and comprehensively optimizing the number, position, cluster length, hole density, hole diameter, construction discharge, fracturing fluid viscosity perforation and construction parameters corresponding to the optimal solution of the whole target function by applying an improved genetic algorithm.

The device further comprises:

the flow calculation unit is used for establishing the relation between the total displacement and the flow of each cluster of all cracks and the flow of the bifurcation cracks;

the pressure calculation unit is used for the relationship between the shaft pressure and the fluid pressure and perforation at the entrance of each cluster of non-planar cracks;

and the overall objective function optimal solution calculation unit is used for calculating a multi-objective function which is matched with the geological characteristics of the reservoir and has large reconstruction volume of each cluster of cracks, balanced expansion of each cluster of cracks and low construction pressure so as to select perforation parameters and construction parameters under geological engineering dessert perforation and pressure limiting constraint conditions.

By the method and the device, high-precision intelligent regulation and control can be performed during fracturing of the unconventional reservoir, and uniform transformation of the unconventional reservoir fractured horizontal well is achieved.

Drawings

FIG. 1 is a schematic diagram of basic parameters of a Changning X well 10 th fracturing segment.

FIG. 2 is a schematic diagram of the effect of different reservoir permeability and porosity parameter combinations on the seam length, the modification volume and the construction pressure.

FIG. 3 is a schematic diagram of 2 clusters of construction parameters optimized by different genetic iteration times and an optimization process under perforation.

FIG. 4 is a schematic diagram of the construction parameters and the optimization process under perforation for 6 clusters of different genetic iteration times.

Detailed Description

The details of the present invention can be more clearly understood in conjunction with the accompanying drawings and the description of the embodiments of the present invention. However, the specific embodiments of the present invention described herein are for the purpose of illustration only and are not to be construed as limiting the invention in any way. Any possible variations based on the present invention may be conceived by the skilled person in the light of the teachings of the present invention, and these should be considered to fall within the scope of the present invention.

The invention provides an intelligent regulation and control method and device for unconventional reservoir fracturing multi-fracture non-equilibrium expansion, wherein the method comprises the following steps:

1. collecting and calculating required basic parameters including reservoir parameters, mechanical parameters, fracturing parameters and shaft parameters, wherein the reservoir parameters include original stratum pore pressure, reservoir permeability, reservoir thickness, comprehensive compression coefficient, reservoir porosity, ground stress and fracture zone parameters, the mechanical parameters include parallel bedding Young modulus, vertical bedding Young modulus, parallel bedding Poisson ratio, vertical bedding Poisson ratio, rock tensile strength, weak surface cohesive force, formation tendency, formation inclination angle and internal friction angle parameters, the fracturing parameters include construction displacement, fracturing fluid viscosity, cluster number, perforation depth, perforation radius, perforation length, perforation density and perforation phase angle parameters, and the shaft parameters include shaft radius, shaft inclination angle and shaft azimuth angle.

2. Based on a perforation cluster fracture initiation and expansion physical model, a perforation cluster fracture initiation and expansion fluid-solid coupling model is established by utilizing a fluid-solid coupling control equation of rock solid deformation, fluid flow in a fracture and matrix pressure diffusion and based on a perforation cluster fracture initiation criterion and an I, II type composite fracture expansion criterion. It should be noted that the model presented in this application is only an example and does not constitute a limitation of the present invention. The person skilled in the art can also establish a multi-fracture propagation model by adopting other methods in the field and perform high-precision intelligent fracture regulation without creative labor.

(1) Deformation of rock solids

Based on a perforation cluster fracture initiation expansion physical model, the relationship among stress, strain and pore pressure of a linear pore elastic medium is given by utilizing a Biot pore elasticity theory:

σ _ij ＝2Gε _ij +2Gvδ _ij ε _kk /(1-2v)-αδ _ij p；(i,j＝1,2,3) (1)

in the formula: epsilon _ij Is a strain tensor, dimensionless; sigma _ij Is the stress tensor, MPa; sigma _kk Is the sum of the total volume stress, MPa; epsilon _kk Is the sum of the total volume strains, dimensionless; p is pore pressure, MPa; g is shear modulus, MPa; nu is Poisson's ratio and is dimensionless; delta _ij Is a Kronecker-delta function and has no dimension; α is the Biot coefficient, dimensionless.

The combined stress and displacement static equilibrium equation is:

ε _ij ＝0.5(u _i,j +u _j,i ) (2)

in the formula: u. of _i,j 、u _j,i Is the displacement component, m.

Combining the static equilibrium equation (2) with the formula (1), neglecting physical strength, and obtaining a Navier equation of pore elasticity:

Gu _i,kk +Gu _k,ki /1-2v-αp _,i ＝0 (3)

the only unknown quantity in equations (1) and (3) is the displacement u, which will be solved by coupling the displacement discontinuity quantities.

(2) Fluid flow in the slot

Assuming that the flow of the fracturing fluid in the crack is Newtonian fluid and laminar flow, the fluid flux in the crack can be calculated according to the Poiseue's law. Considering that the fracturing fluid of an injection source is not compressible, based on the fluid substance balance principle in the fracture, the flow of the fracturing fluid in any hydraulic fracture infinitesimal in the length direction should be equal to the sum of the volume change of the fracture infinitesimal and the filtration loss of the fracturing fluid:

in the formula: q _k Flow of the fracturing fluid of the kth hydraulic fracture, m ³ /min；c _L Is the fluid loss coefficient, m/min ^0.5 (ii) a t is the operation injection time, min; q. q.s _L The fluid loss rate of the fracturing fluid at the hydraulic fracture (x, y) at the time t, m ² /min。

The first item on the left side in the formula (4) represents the flow change of inflow and outflow fracturing fluids in any hydraulic fracture infinitesimal; the first term on the right represents the injection flow strength; the second term on the right represents the fracturing fluid filtrate loss; the third term on the right represents the amount of fracture infinitesimal volume change.

(3) Pressure diffusion of matrix

According to the principle of conservation of mass, the net fluid flow from the point source of the substrate given by the fluid mass balance equation is equal to the sum of the fluid mass increase in the pore space and the injected/discharged fluid, neglecting

And assuming that the pore volume changes very little, there are:

in the formula: k is the matrix permeability, μm ² ；μ _g Is the viscosity of the matrix fluid, mPa · s; q _s For point source injection or production flow, m ³ /min；c _f Is the reservoir compressibility factor, MPa ^-1 (ii) a Alpha is the hole elasticity coefficient and is dimensionless; phi is porosity, dimensionless.

(4) Fracture initiation propagation criterion

When the bottom hole pressure reaches the fracture initiation pressure of the k perforation cluster fracture, namely the fracture initiation condition is met, the fracture propagation can be regarded as the type I and type II composite fracture problem. Discontinuity amount of displacement near the tip of the crackCan calculate K _I 、K _II The stress intensity factor, and then the propagation direction of the crack can be calculated. Based on the fracture mechanics theory, when the condition that the equivalent stress intensity factor is equal to the critical stress intensity factor is met, the crack is expanded, namely:

in the formula:

K _I is a type I strength factor, MPa.m ^1/2 ；K _II Is a type II intensity factor, MPa.m ^1/2 (ii) a a is half of the extension step, m; e is Young's modulus, MPa; v is Poisson's ratio, dimensionless; d _n 、D _s The corresponding displacement discontinuity, m, for the crack propagation tip.

And (3) solving a first derivative on the left side of the formula (6) to make the value of the first derivative be 0, and simultaneously ensuring that the second derivative is less than 0 to obtain a perforation crack extension direction angle omega as follows:

3. based on a perforation cluster initiation and expansion fluid-solid coupling model, a displacement discontinuity method is adopted to solve rock deformation, pore pressure diffusion and multi-field coupling variables of fluid flow in the seams, discrete point source normal displacement discontinuity, tangential displacement discontinuity, filtration loss, pressure in the seams and perforation crack extension direction angles of the perforation cluster cracks are calculated, and then the discrete point source normal displacement discontinuity, the tangential displacement discontinuity, the filtration loss, the pressure in the seams and the perforation crack extension direction angles are substituted into a fluid flow equation in the seams to obtain the length, the width, the pressure in the seams and the expansion track.

The fracture initiation and propagation coupling is to synthesize fluid-solid coupling control equations of rock solid deformation (strain and stress), fluid flow in the fracture (fracture internal pressure and fracture width), matrix pressure diffusion (pore pressure) and the like, and perform fracture initiation and propagation under the condition of meeting the fracture propagation criterion. Wherein the common variable acting as a bridging link is the amount of undisplaced discontinuity. Stress, strain, internal pressure of jointForce, gap width, pore pressure and displacement discontinuity D _n 、D _s Are directly related.

(1) Solid deformation equation stress and displacement solution

For the plane strain condition, a single fracture section of length 2a begins injecting fracturing fluid at a constant flow rate at time t-0. The inner boundary and the outer boundary condition are determined according to the initial condition and the formulas (3) and (4) and are defined as the formulas (8) to (10).

The initial conditions are given by:

t＝0,(x,y):p＝0,u _x ＝u _y ＝0,σ _xx ＝σ _yy ＝σ _xy ＝0 (8)

inner boundary conditions:

y＝0,|x|≤a:u _x (x,0 ^- )-u _x (x,0 ⁺ )＝D _s ,u _y (x,0 ^- )-u _y (x,0 ⁺ )＝D _n ,Q _s ＝-Q _k (9)

the outer boundary condition is

The basic solution of the stress strain and the displacement of the elastic rock of the hole is obtained by adopting a displacement discontinuity method and combining fluid constant flow injection/output, boundary conditions and a displacement discontinuity solution as follows.

Displacement induced by continuous unit fluid source along straight fracture section

And stress

Solve as follows, let:

r ² ＝(x-x′) ² +y ² (11)

where x' varies from-a to + a. Defining:

order:

the fluid source induced displacement is:

the fluid source induced stress is:

wherein the displacement is caused by discrete amounts of continuous cell normal displacement along a linear fracture segment

And stress

The solution is as follows:

normal displacement discontinuity induced displacement:

normal displacement discontinuity induced stress:

displacement caused by discontinuity of continuous unit shear displacement along a linear fracture

And stress

The solution is as follows:

discontinuous amount induced displacement of tangential displacement:

tangential displacement discontinuity induced stress:

in the formula: d _n And D _s As sources of discontinuity in normal and tangential displacements, m; q. q.s _L The fluid loss velocity in the fracture (interfacial flow velocity between fracture and matrix), m ³ Min; superscripts dn, ds, and q are for normal displacement discontinuity, tangential displacement discontinuity, and fluid source, respectively;

displacement in x and y directions, m;

is the fluid source stress component, MPa;

x, y direction displacements, m, generated for the displacement discontinuity;

the stress component is the stress component generated by the displacement discontinuity, MPa.

(2) Matrix pressure diffusion solution

By adopting a displacement discontinuity method, injecting a displacement discontinuity quantity induced pore pressure into a fracture section fluid at a constant flow rate is as follows:

p(x,y,t)＝p ^dn (x,y,t)D _n +p ^ds (x,y,t)D _s +p ^q (x,y,t)q _L (29)

in the formula: p is a radical of ^dn 、p ^ds For generating discrete amounts of displacementThe induced pore pressure of (3), MPa; p is a radical of ^q Induced pore pressure, MPa, generated for constant velocity fluid injection/production at the fracture section.

Wherein the fluid source induces a pore pressure p ^q Given by:

normal displacement discontinuity induced pore pressure p ^dn ：

Discontinuous amount of tangential displacement induced pore pressure p ^ds ：

In the formula: v. of _u The Poisson's ratio of the rock under the condition of no water drainage is dimensionless and is generally 0.31.

(3) Fluid flow equation dispersion within a slot

The fluid flow equation in the gap established according to the formula (4) has strong nonlinear characteristics, and is solved after being dispersed by adopting a finite difference method. The finite difference equation of the hydraulic fracture mass balance equation is as follows:

in formula (54), Δ x ═ x _i+1/2 -x _i-1/2 ，Δt＝t ⁿ -t ^n-1 N represents the current time and n-1 represents the previous time, further discretizing equations rearranging terms:

(4) displacement discontinuity method coupling solution

The solution to the displacement discontinuity amount is space and time dependent, with the basic solution based on constant displacement discontinuity and constant interface fluid loss velocity. In numerically integrating the source intensity at each time step, all previous increments of source intensity must be included. The induced stress and pore pressure increase on the ith fracture section were calculated from the source strength increase:

in the formula:

and

for the current time τ _ξ Source intensity increment for jth fracture section;

and

is tau _h The source intensity increment for the jth fracture segment in time, h summed from 1 to ξ -1;

and

for the jth broken cell at time step τ _h The influence coefficient of the conversion from the formula (24) to the formula (32) to the i-th cleavage unit.

The normal induced stress increment is used for establishing an equation set as follows:

the tangential induced stress increment is used for establishing an equation set as follows:

the system of equations is constructed using the induced pore pressure increase as:

the flow equation (34) within the joint seam translates into a displacement discontinuity form:

in the formula: a is half of the extension step, m;

to represent tau _ξ The ith gap width at time is tau _ξ-1 Discontinuity amount of normal displacement

And τ _ξ Discontinuous increment of time normal displacement

Sum, m;

to represent tau _ξ The pressure in the ith point seam is MPa;

to represent tau _ξ Fluid loss velocity at point i in time, m ² /min。

In the formula (36) -in the formula (39) with τ _h The parameters at time are all known quantities at _ξ Ginseng of all agesThe quantities are all unknown quantities, thereby building a system of nonlinear equations:

the above equation is a 4N-dimensional nonlinear system of equations, and the component form of F is:

solving the nonlinear equation set (40) by using a Newton Raphson iterative method, and enabling x ^(k) ＝[x ₁ ^(k) ,x ₂ ^(k) ,…,x _4N ^(k) ] ^T Component f of function F (x) _i (x) (i ═ 1, …,4N) in x ^(k) Using a multivariate function Taylor expansion and taking the linear part thereof, the Newton iteration format for solving the nonlinear equation set (40) is:

x ^(k+1) ＝x ^(k) -F′(x ^(k) ) ^-1 F(x ^(k) )(k＝0,1,…) (42)

thus, the adjacent iteration roots | x according to equation (42) ^(k) -x ^(k-1) If the | norm is less than a certain error, the discontinuous increment delta D of the discrete point source displacement of the perforation cluster crack can be obtained by solving _s 、ΔD _n 、q _L 、p _f And the extension direction angle of the perforation crack is substituted into the flow equation of the fluid in the crack to obtain the length, width, pressure and expansion track of the crack.

4. On the basis of a planar perforation cluster fracture initiation-expansion coupling model and multi-field coupling variable solving, reservoir physical properties, ground stress heterogeneity, perforation hole friction resistance and the like are considered, dynamic flow distribution, initial perforation cluster extension and fracture induced stress in the process of multi-fracture competition initiation and expansion are coupled, the planar multi-cluster fracture initiation-expansion coupling model is established, and the flow, the pressure and the multi-field coupling variable of each cluster of fractures are obtained by adopting a Newtonian Raphson method.

(1) Hydraulic fracture induced stress

The hydraulic fracture propagation stress field distribution is the superposition of multi-fracture induced stress, in-situ stress and anisotropic induced stress; supposing that M perforation clusters are arranged, each perforation cluster is dispersed into N infinitesimal sections, reservoir physical properties, ground stress heterogeneity and the like are considered on the basis of a perforation cluster expansion induced stress field formula (40), and normal stress and shear stress components in (x, y) coordinates obtained by superposing multi-crack induced stresses are as follows:

in the formula: sigma _mx 、σ _my 、τ _mxy Generating a normal stress component and a shearing stress component, MPa, for the perforation cluster fracture under x and y coordinates; sigma _H 、σ _h 、σ _v Maximum horizontal principal stress, minimum horizontal principal stress and vertical stress, MPa; beta and psi are the well bore azimuth angle and the well bore inclination angle; sigma _xciso 、σ _yciso 、τ _xyciso Anisotropic well circumferential stress in x, y and xy directions, MPa;

and (3) expanding the m-th perforation cluster crack to the i-th section of positive and tangential induced stress (MPa) at the t moment under x and y coordinates.

On the basis of the calculation of the single-cluster crack induced stress (35), the multi-cluster crack induced stress fields are superposed to obtain the multi-cluster crack induced stress field. And if M clusters of cracks are arranged, the induced stress of the mth cluster of j points in the multi-cluster cracks is as follows:

in the formula:

and

the seismic source intensity increment of the jth fracture section at the current time step;

and

summing from 1 to ξ -1 for the previous source intensity increment for the jth break at time step h;

it is the case that equation (35) is switched to the influence coefficient of the jth broken cell on the ith broken cell at time step h in the case of multi-cluster.

(2) Planar multi-cluster fracture flow dynamic allocation

Based on the Kirchoff first law, when the staged multi-cluster fracturing of the horizontal well is carried out, the total displacement of the fracturing pump is Q _t The total flow is divided into clusters, and the displacement of the k perforation cluster of the horizontal well is expanded to the j section is Q _k,j The total displacement of fluid is equal to the sum of the displacements per cluster of all fractures, i.e.:

based on Kirchoff's second law, the fluid pressure balance criteria in the horizontal wellbore are established using the point of reference in the horizontal wellbore (O-target). The wellbore pressure of the O target point is equal to the sum of the fluid pressure at the inlet of each cluster of cracks and the friction resistance of the perforation holes:

p _w ＝p _pfk,j +p _fk,j (j＝1,2,…,n) (46)

in the formula: p is a radical of _w,k The pressure of a wellbore fluid when the horizontal well crack k expands to the jth section is MPa; p is a radical of _pfk,j Friction resistance pressure of a perforation hole is MPa when the horizontal well crack k expands to the jth section; p is a radical of _fk,j For the fracture k of the horizontal well to expand to the jth sectionBody pressure, MPa.

Eyelet friction resistance p _pfk,j The calculation formula of (a) is as follows:

in the formula: rho _s In terms of fracturing fluid density, kg/m ³ ；N _pk The number of the perforation clusters of the k perforation of the horizontal well and the density of the perforation holes G _k Length L of perforation cluster _k Product of (A), N _pk ＝G _k L _k ，m；d _pk The diameter of a perforation cluster hole of a horizontal well k is m; c _d The flow coefficient of the horizontal well perforation is zero dimension.

Fracture fluid pressure when horizontal well fracture k expands to jth section:

p _fk,j ＝p _netk,j +σ _n,k,j (48)

wherein the net pressure in the gap is:

p _netk,j ＝p _fk,j -σ _n,k,j (49)

the circumferential stress of the fracture end when the horizontal well perforation cluster fracture k expands to the jth section is as follows:

in the formula: omega _k,j The angle between the fracture extension direction and the horizontal direction when the fracture k of the perforation cluster expands to the j section is degree.

Simultaneous equations (45) and (48) can be constructed by pressure-flow coupling of a nonlinear system of equations:

F(Q ₁ ,Q ₂ ,…,Q _M ,p _w )＝0 (51)

the above equation is a nonlinear system of equations in dimension M +1, and the component form of F is:

solving the nonlinear equation set (51) by using the Newton Raphson iteration method can be written as

F(x ₁ ,x ₂ ,…,x _M ,x _M+1 )＝0 (53)

Let x ^(k) ＝[x ₁ ^(k) ,x ₂ ^(k) ,…,x _M+1 ^(k) ] ^T Component f of function F (x) _i (x) (i ═ 1, …, M) in x ^(k) Using a multivariate function taylor expansion and taking the linear part of the taylor expansion, the newton iteration format for solving the nonlinear system of equations (51) is:

x ^(k+1) ＝x ^(k) -F′(x ^(k) ) ^-1 F(x ^(k) )(k＝0,1,…) (54)

thus, the adjacent iteration roots | x according to equation (54) ^(k) -x ^(k-1) If the | | norm is less than a certain error, the flow and pressure of each perforation cluster and the multi-field coupling variables of the cracks of each cluster can be solved.

5. On the basis of a planar multi-cluster fracture initiation-expansion coupling model and multi-field coupling variable solving, taking non-planar fracture network fluid flow into consideration, dynamically distributing coupling multi-fracture and bifurcation fracture flow, and establishing a non-planar fracture initiation-expansion coupling model based on a hydraulic fracture and natural fracture interaction judgment criterion taking anisotropic induced stress into consideration; realizing non-synchronous fracture initiation-expansion simulation of the non-planar fracture; and establishing a non-planar multi-cluster fracture reconstruction volume model based on the shear slip and tensile failure mechanical conditions of reservoir rock and a fracture network permeability model according to the change of the reservoir pore pressure field caused by the induction of the non-planar multi-cluster fractures.

(1) Hydraulic fracture and natural fracture interaction

The hydraulic fractures induce stresses in the formation during propagation in the formation, and thus the stress applied to the natural fractures is effectively a superposition of the in situ stress and the induced stresses created by the hydraulic fractures. When the hydraulic fracture is close to the natural fracture, the stress distribution on the natural fracture surface is that the approach angle of the hydraulic fracture and the natural fracture is beta _NF . The normal stress and the shear stress on the natural crack surface are as follows:

in the formula: sigma _xxNF 、σ _yyNF 、τ _xyNF The normal stress and the shear stress components of one point of the natural crack surface are MPa; sigma _H 、σ _h The maximum and minimum horizontal principal stress of the stratum, MPa; k is _I Is a type I stress intensity factor, MPa.m ^0.5 (ii) a r is the distance from the tip of the hydraulic fracture to one point of the natural fracture surface, m; theta is an included angle between the direction of the maximum horizontal principal stress and a connecting line from the tip of the hydraulic fracture to one point of the natural fracture surface of the natural fracture.

Order to

Conversion of the stress in equation (55) to the coordinate β _xNF And beta _yNF Then natural fracture face stress:

when a hydraulic fracture interacts with a natural fracture, the natural fracture may open, shear slip, and cross through, all of which may greatly affect the hydraulic fracture propagation path. Thus, fracture criteria models for opening, shearing, and crossing of the natural fracture are established separately as the hydraulic fracture interacts with the natural fracture.

Natural fracture reopening

When the positive stress at the fracture tip position is less than the fluid pressure in the mth hydraulic fracture, the natural fracture will open:

p _f ≥σ _n,m (57)

the fluid pressure within the fracture is the sum of the fracture tip circumferential stress and the fracture internal net pressure:

p _f ＝σ _n,m +p _net,m (58)

when tensile failure occurs when the natural fracture surface satisfies the condition (57), the opening width of the natural fracture is as follows according to the theory of elastic mechanics:

w _m ＝2(1-v)(p _f -σ _n,m )H _NF,s /E (59)

in the formula: w is a _m M is the opening width of the natural fracture; h _NF,s Is the s-th natural fracture height, m.

② natural crack slippage

The critical conditions for the shear slip of the natural fracture are as follows:

|τ _βNF,m |＞s _O -μ _f (σ _n,m -p _f ) (60)

when the (60) natural fracture slip condition is met, the shear displacement is:

in the formula: s ₀ The cohesive force of the natural crack surface is MPa; u. of _s,k M is the shear displacement; k is a radical of _o Is a Kolosov constant, k _o 3-4 v; g is shear modulus, G ═ E/2/(1+ v), MPa; x is the coordinate of any point on the crack surface, m; l is half the length of the natural crack, m.

Thirdly, hydraulic fracture passes through natural fracture

The hydraulic fracture will pass through the natural fracture conditions.

Critical transit time sigma ₁ To achieve tensile strength sigma _T ：

σ ₁ ＝σ _T (63)

In addition to satisfying equation (63), | τ must also be satisfied _βNF,m |＞s _O -μ _f (σ _n,m -p _f ) A condition that, when both conditions are met, the hydraulic fracture will continue to extend through the natural fracture. Considering the common influence of the anisotropy induced stress and the in-situ stress induced stress, let T ═ sigma ^ sigma _T -(σ _H +σ _xciso -σ _h -σ _yciso ) And/2, substituting the formulas (57) to (62) into the formula (63) to obtain:

critical distance r corresponding to the solution of equation (64) _c ：

The steering angle corresponding to the solution of equation (64) is γ _i ：

In the formula: gamma is the crack steering angle, deg..

(2) Non-planar multi-cluster fracture flow dynamic allocation

The flow dynamic distribution of the multi-cluster non-planar fractures is different from the flow dynamic distribution of the multi-cluster planar fractures in that the judgment criterion of hydraulic fractures penetrating through natural fractures, the flow distribution of the hydraulic fractures penetrating through the bifurcate fractures of the natural fractures, the pressure balance of nodes, the induced stress of activated natural fractures and the filtration loss are increased. Therefore, the factors are considered on the basis of the dynamic distribution of the multi-cluster plane fracture flow.

Flow balance

Based on the Kirchoff first law, when the staged multi-cluster fracturing of the horizontal well is carried out, the total displacement of the fracturing pump is Q _t The total flow is divided into clusters, the hydraulic fracture of the kth perforation cluster of the horizontal well is expanded to the jth section, and the displacement is Q _k，j And the total discharge capacity of the fluid is equal to the sum of the flow rates of all the fractures in each cluster, and the hydraulic fractures are penetrated by natural fractures to form bifurcated fractures, wherein the flow rates are as follows:

in the formula: q. q of _Si Flow rate of the hydraulic fracture to the previous step size fracture from the ith node of n nodes formed by the hydraulic fracture through the S natural fracture，m ³ /min；T _ki Is a transmissible geometric factor without dimension; lambda is a fluid flow resistance coefficient, and is dimensionless; p is a radical of _fk 、p _fi The pressure of fluid in the current i-node seam of the kth cluster in the previous step length is MPa.

For a split with n splits, the transitive geometry factor is:

in the formula: alpha is alpha _k Coefficient of flow Capacity of the kth bifurcation fracture, m ³ . Wherein:

in the formula: w is a _NF The width of the natural crack is m; a. the _i Is the i-th bifurcation crack area, m ² ；D _i Is the distance, m, of the ith bifurcation fracture from the central node.

Pressure balance

Based on Kirchoff's second law, the fluid pressure balance criteria in the horizontal wellbore are established using the horizontal wellbore (O target) as a reference point. The wellbore pressure at the O-target is equal to the sum of the fluid pressure at the entrance of each cluster of non-planar fractures and the friction at the perforation:

p _w ＝p _pfk,j +p _fk,j (j＝1,2,…,N) (70)

fracture fluid pressure when horizontal well fracture k expands to jth section:

p _fk,j ＝p _netk,j +σ _nk,j (71)

wherein the net pressure in the gap is:

p _netk,j ＝p _fk,j -σ _nk,j (72)

in the formula: sigma _mxNF 、σ _myNF 、τ _mxyNF Generating normal stress and shear stress components, MPa, for the non-planar fracture under x and y coordinates; beta is a beta _k,j The included angle between the extension direction of the cracks and the horizontal direction when the cracks of the k perforation cluster extend to the j section is degree.

The multi-cluster non-planar fracture propagation stress is the superposition of multi-fracture induced stress, in-situ stress, anisotropic induced stress and activated natural fracture induced stress, and the natural fracture induced stress is superposed on the basis of a multi-cluster propagation stress field formula (44) to obtain the non-planar multi-cluster fracture normal stress and the shear stress components under x and y coordinates as follows:

in the formula: sigma _H 、σ _h 、σ _v Maximum horizontal principal stress, minimum horizontal principal stress and vertical stress, MPa; beta and psi are the well bore azimuth angle and the well bore inclination angle; sigma _xciso 、σ _yciso 、τ _xyciso The calculated anisotropic well circumferential stress in the x, y and xy directions is MPa;

the non-planar fracture of the mth perforation cluster is expanded to the ith section of positive and tangential induced stress, MPa, at the time t under x and y coordinates; theta _NF,s Is the azimuth of the s-th natural fracture.

On the basis of the calculation of the induced stress of the multiple clusters of fractures, each activated natural fracture is regarded as the induced stress generated discontinuously by the displacement of a single cluster of fractures; and superposing the multiple clusters of crack induced stress fields to activate the natural crack induced stress to obtain the multiple clusters of non-planar crack induced stress fields. S natural cracks are activated, and the m & ltth & gt point i & lt/th & gt induced stress in the non-planar natural cracks is as follows:

the equations (45), (68) and (70) may form a nonlinear system of equations:

F(Q ₁ ,Q ₂ ,…,Q _M ,p _w ,q ₁₁ ,q ₁₂ ,…,q _Sn ,Q _1(j-1) ,…,Q _S(j-1) ,p _f11 ,p _f12 ,…,p _fSn )＝0 (76)

the above formula is a (M +1) + S x (2n +1) -dimensional nonlinear system of equations, and the component form of F is:

in the formula:

and

and

is the first step of time hPrevious source intensity increments for j breakouts, summed from 1 to ξ -1;

the influence coefficient of the jth fracture unit on the ith fracture unit in the time step h is shown; tau is _NFξ Current time to traverse the natural fracture, s; tau is _NFh To pass through natural fractures _NFh Time, s.

The non-linear equation set (76) is solved by using the Newton Raphson iteration method, and can be written as

F(x ₁ ,x ₂ ,…,x _M+1+Sn+Sn+S )＝0 (78)

Let x ^(k) ＝[x ₁ ^(k) ,x ₂ ^(k) ,…,x _M+1+2Sn+S ^(k) ] ^T Component f of function F (x) _i (x) (i ═ 1, …, M +1+2Sn + S) in x ^(k) Using a multivariate function Taylor expansion and taking the linear part thereof, the Newton iteration format for solving the nonlinear equation set (78) is:

x ^(k+1) ＝x ^(k) -F′(x ^(k) ) ^-1 F(x ^(k) )(k＝0,1,…) (79)

thus, the adjacent iteration roots | | | x according to equation (79) ^(k) -x ^(k-1) If the | | norm is less than a certain error, the flow and pressure of each perforation cluster and the flow and pressure of the bifurcation fracture can be solved.

(3) Non-planar multi-cluster fracture reformation volume prediction

Firstly, dispersing the formed ith hydraulic fracture into a section u (i) at a spacing a, and then dispersing the natural fracture network into a section v (k) at a spacing b; each section is regarded as the pressure change generated by point source injection to the plane, and based on the pressure drop superposition principle, the pressure generated by the non-planar fracture fluid to the plane is as follows:

in the formula: p (x, y) is the pore pressure at the plane (x, y) after the formation of a non-planar fracture pressure response；p _fi,j (x, y) calculating the pressure in the jth section of the ith perforation cluster fracture according to a formula (79); p is a radical of _fi,j And (x, y) calculating the fracture pressure of the jth section of the ith natural fracture at the plane (x, y) according to the formula (79).

The pore pressure after the pressure of the planar multi-cluster non-planar fracture is responded is obtained according to a formula (80), and then the non-planar fracture reconstruction volume simulation can be carried out by utilizing the multi-cluster fracture reconstruction volume calculation process to superpose and activate the natural fracture reconstruction volume in combination with the opening shear failure criterion of the planar weak point.

6. Based on a non-planar crack initiation and propagation mechanism model, adopting DOE design to establish a sample database, training by using a Gauss-Kriging machine learning model, and establishing a multi-input multi-output intelligent agent model for simulating the mechanism model; on the basis, the restriction of large reconstruction volume of each cluster of cracks matched with the geological characteristics of the reservoir, low construction pressure, balanced crack length expansion, geological engineering desserts and the like is taken as a multi-target function; and comprehensively solving the optimal solution of the whole objective function by using an improved genetic algorithm, and establishing the unconventional reservoir fracturing multi-fracture unbalanced expansion intelligent control method.

(1) Multi-input multi-output intelligent agent model

The research of the intelligent agent model is developed, and firstly, the design of an experimental scheme is required to be carried out, and a sample large database is constructed. Carrying out test design and database establishment according to parameter ranges such as vertical crustal stress, minimum horizontal main stress, maximum horizontal main stress, Young modulus, Poisson's ratio, permeability, porosity, discharge capacity, fracturing fluid viscosity, perforation diameter, perforation density, perforation cluster length, perforation position and the like of a horizontal section of a fracturing well; for vertical stress x ₁ Minimum horizontal principal stress x ₂ Maximum horizontal principal stress x ₃ Young's modulus x ₄ Poisson's ratio x ₅ Permeability x ₆ Porosity x ₇ Discharge x ₈ Viscosity x of fracturing fluid ₉ Diameter x of the hole ₁₀ Cell density x ₁₁ Length x of perforation cluster ₁₂ Perforation cluster position x ₁₃ Equal independent variables are used as input parameters, DOE experimental Design (Design of experiment) is adopted, and the output parameters are based on the non-plane established in the frontCrack propagation mechanism model calculation of crack length y of each design point under DOE (design of experiment) design parameters ₁ Improved volume y ₂ Pressure y ₃ And the like. According to design point S (S) ₁ ,s ₂ ,...,s _m ) And a training sample big database X is constructed, so that the precision of the training model adopting the sample set is fully ensured.

The input parameter training sample big database X is:

in the formula: x is an input parameter training sample big database matrix; x is the number of _j (i) J is more than or equal to 1 and less than or equal to n, and i is more than or equal to 1 and less than or equal to m, which are the ith design point of the jth component representing x.

The big database Y of the output parameter training sample is as follows:

in the formula: y is a large database matrix of the output parameter training sample; y is _q (i) To represent the target value of the ith design point of the qth component of y, 1. ltoreq. q, 1. ltoreq. i. ltoreq.m.

Aiming at the research of the unbalanced expansion regulation and control of the multi-cluster perforation cracks, the method is based on multiple input parameters such as geology, construction parameters and the like, takes multiple output parameters such as large reconstruction volume, low construction pressure and the like of each cluster of cracks matched with reservoir geological characteristics as a multi-target function, and solves the problems of balanced crack length expansion and balanced geological engineering desserts as constraints. Approximate processing for a multiple-input multiple-output Problem (Multiinput-Multioutput Problem) belongs to the category of multiple-input multiple-output intelligent agents; a multi-input multi-output based experimental scheme is subjected to Kriging machine learning training and MSE (minimum mean square error) constraint by adopting 0-order, 1-order and 2-order regression models and Gauss correlation models, and a multi-input multi-output intelligent agent model is established.

Approximate agent model for crack initiation and propagation of anisotropic unconventional reservoir non-planar fractures is established by adopting Kriging intelligent algorithm. Given m design points S ═ S ₁ ,...,s _m ] ^T (s _i ∈R ⁿ ) And its response Y ═ Y ₁ ,...,y _m ] ^T (y _i ∈R ^q ) Combinations of (a) and (b). Inputting the n dimensions on the assumption that the data are subject to the standardization condition

The deterministic output y (x) e R ^q Expressed as an approximate combination of a regression model F and a random function z (stochastic process):

wherein F is a function using p preference functions

The regression model of the linear combination of (1):

wherein the coefficient { beta _k,l Is the regression parameter.

Assuming that the mean of the random process z is 0, the covariance between z (w) and z (x) is:

wherein

Is the process variance of the system response to the ith component,

representing the correlation model with the parameter theta.

Equation (85) is approximately equivalent to the regression model for the objective function plus a preferred stochastic process z bias.

And the true value of y can be written as:

y _l (x)＝F(β _:,l ,x)+α(β _:,l ,x) (86)

where alpha represents the approximation error.

Thus, the intelligent agent process is to approximate the model in the target region (e.g., x ∈ D) by choosing β appropriately

And true value y _l (x) The error of (2) is minimized.

Define R as the random process correlation coefficient between z at different design points:

at the untested point x, let:

in the formula: r (x) is a correlation coefficient between z representing each design point and x.

Taking linear correlation prediction as an example for introduction, the linear predictor is:

wherein c ═ c (x) e R ^m . The error is as follows:

wherein Z is [ Z ] ₁ ,...,z _m ] ^T The error at each design point is shown. To ensure that the predictor is unbiased, it is necessary to make F ^T c-F (x) 0, or F ^T c(x)＝f(x)。

Under this condition, the Mean Square Error (MSE) of the predictor (89) is:

at F ^T c-F (x) 0 or F ^T c (x) f (x) under the constraint of _c Error for independent variables

Minimizing, constructing Lagrange function, and obtaining according to the first-order necessary condition of optimality:

its generalized least squares solution (relative to R) is:

β ^* ＝(F ^T R ^-1 F) ^-1 F ^T R ^-1 Y (93)

substituting equation (93) into equation (92) yields the expression for the predictor:

for multiple response cases (i.e. q)>1) Equation (140) holds for each column in Y. Therefore, equation (94) can be obtained by calculating β from equation (93) ^* ∈R ^p×q And from the residual R gamma ^* ＝Y-Fβ ^* Calculating gamma ^* ∈R ^m×q 。

While for a fixed set of design points, the matrix beta ^* And gamma ^* Is also fixed. For each new x, only the vector f (x) e R needs to be computed ^p And R (x) ε R ^m And two simple products are added.

The following expression for MSE of the crikin predictor:

wherein u ═ F ^T R ^-1 r-f and σ ² Obtained by generalized least squares fitting. Equation (95) can be generalized directly to the multi-response case: for the ith response, σ may be replaced by σ _l (process variance of the ith response function).

When n is 1, the model is a 1 st order regression model. Let x _j The j-th component representing x, constant term, p ═ 1:

f ₁ (x)＝1 (96)

when the first term p is n + 1:

f ₁ (x)＝1，f ₂ (x)＝x ₁ ，...，f _n+1 (x)＝x _n (97)

the second order term is the sum of the first order,

f ₁ (x)＝1 (98)

f ₂ (x)＝x ₁ ，...，f _n+1 (x)＝x _n (99)

the corresponding jacobian determinant is (subscript n × q denotes the size of the matrix, O denotes the all 0 matrix):

constant term: j. the design is a square _f ＝[O _n×1 ]The first item: j. the design is a square _f ＝[O _n×1 I _n×n ]The second order term: j. the design is a square _f ＝[O _n×1 I _n×n H]

When n is>The model is an n-order regression model at 1. For H ∈ R ^n×(p-n-1) And n is 2, 3, the H matrix is:

when n is 2:

when n is 3:

the correlation model is a product of static first-order correlation models, taking into account a correlation model of the form:

the correlation function is a non-linear relationship, such as Gauss, which is approximate to a parabolic relationship, wherein the optimal coefficient theta of the Gauss correlation function ^* This can be obtained by solving the following equation:

where | R | is the determinant of R; theta ^* Is the optimal coefficient.

(2) High-dimensional multi-objective function and constraint condition

Aiming at the research of the unbalanced expansion regulation and control of the multi-cluster perforation cracks, the method aims at solving the problem of multi-target optimization of constraints such as large modification volume of each cluster of cracks, low construction pressure, balanced crack length expansion, geological engineering desserts and the like which are matched with the geological characteristics of a reservoir. Therefore, the unconventional reservoir multi-cluster perforation crack parameter multi-objective optimization function established by the method is as follows:

min{H ₁ [x(ξ)],H ₂ [x(ξ)],H ₃ [x(ξ)]，,H ₄ [x(ξ)]} (105)

the multi-objective optimization function in the formula (105) represents the optimized discharge capacity, the viscosity of the fracturing fluid, the diameter of the perforation, the density of the perforation and the length of the perforation cluster under the condition of meeting the multi-objective optimization and constraint conditions, and has the characteristics of high-dimensional multi-input parameters, multi-output parameters, high nonlinearity and irreducibility; the first term indicates that the engineered dessert proximity for the selected perforation location is small; the second term represents the balanced expansion of the crack length of each cluster; the third item represents that each cluster is large in modification volume; the fourth term indicates that the construction pressure is small. Wherein:

in the formula: xi is an independent variable which represents formation parameters (a minimum level principal stress range, a maximum level principal stress, a Young modulus, a Poisson ratio, a permeability, a porosity and the like) corresponding to the optimized positions of each cluster of perforation and is constructed in the optimized discharge capacity, the viscosity of fracturing fluid, the diameter of each perforation, the density of each perforation and the length of each perforation cluster; d [ x (xi)]Calculating the fracturing closeness of geological engineering at the perforation position; std { } is a variance function representing the length of all the cluster seams, and the smaller the variance is, the more balanced the seam length expansion is;

to represent the independent variable and the slot length y for the database x ₁ Improved volume y ₂ Construction pressure y ₃ The dependent variable adopts a regression weight coefficient established by a Gauss-Kriging agent model;

to represent the independent variable and the slot length y for the database x ₁ Improved volume y ₂ Construction pressure y ₃ The dependent variable adopts a correlation weight coefficient established by a Gauss-Kriging agent model; sigma _n [x(ξ)]The circumferential stress increment generated by each cluster of induced stresses is calculated by using formula (73).

Constraint conditions are as follows:

the first term in the formula (107) represents that the sum of the flow rates of all clusters of cracks is equal to the total construction discharge capacity Q _t The constraint of (2); the second term represents the constraint that the bottom hole pressures of the clusters are equal; the third term represents a constraint condition that the cluster spacing is greater than the minimum value of the cluster spacing; the fourth item represents the formation parameters, construction parameters and perforation parameters corresponding to each clusterAnd a constraint between the minimum and maximum values of the fractured well.

(3) High-dimensional multi-objective parallel evolutionary optimization method

In order to solve the nonlinear multi-objective functions and constraint functions, a sequential linear programming, a sequential quadratic programming, a sequence unconstrained minimization method, a genetic algorithm and an improved genetic algorithm are developed based on an optimization mathematical principle. The genetic algorithm generally has a slow convergence speed, but can successfully find a reliable optimal solution of the high variable noise problem. However, multiple objectives are inherently conflicting with each other, and there is rarely a single optimal solution for multiple objectives, but rather a set of solutions known as Pareto optimal sets. The high-dimensional multi-objective parallel evolution optimization method is naturally applicable to the problems because a group of solutions can be evolved to approach the whole Pareto optimal set in one operation.

Therefore, the unconventional reservoir multi-cluster perforation crack parameter multi-objective optimization problem is decomposed into a plurality of sub-optimization problems, each sub-optimization problem comprises a few objective functions of the original optimization problem and also comprises an objective function formed by aggregating other objective functions, and the problem solving difficulty is reduced; secondly, the improved multi-target genetic optimization algorithm is improved in fitness, individual intensity calculation and environment selection, the problems that a boundary solution in the genetic algorithm is shifted out of a population and the convergence speed is low are solved, and the parallel evolutionary optimization algorithm developed aiming at the high-dimensional multi-target function comprises the following steps.

Assuming that N is the size of the population,

for non-dominated solution sets, t _max Is the maximum number of iterations.

The method comprises the following steps: firstly, parallel evolutionary algorithm coding is carried out aiming at multi-target problems, and an initial population P is set ₀ Size, non-dominated solution set

The time t is 0. Secondly, the fitness is calculated by combining the strength value, the original fitness and the neighborhood density parameter to obtain a dominant solution set and a non-dominant solution setAnd (5) matching the solution set. ③ P _t And

all non-dominant solution sets in (1) are copied to

And selecting the environment. Fourthly, the number of iterations t>t _max When the termination condition is reached, outputting the decision variable set represented by the non-dominant solution to P _t+1 Otherwise, the algorithm continues. Fifthly, to

A competitive selection is made to place in the preferred set. Sixthly, the individuals in the preferred set are crossed and mutated, and the result is stored

Adding one iteration to two until the maximum iteration time t _max 。

The specific evolutionary algorithm coding, fitness calculation and environment variable selection method comprises the following steps:

coding mode of parallel evolutionary algorithm

Using a random sequence b, coded in decimal notation ₁ b ₂ b ₃ …b ₂₄ As NP (i) chromosome (all NP (i) chromosomes constitute Pareto feasible solutions), where b ₁ For fracturing fluid viscosity corresponding to x ₉ (ξ)；b ₂ The number of the perforation clusters corresponds to xi; b _3-12 The discharge capacity of the 1 st-10 th perforation cluster corresponds to x ₈ (ξ)；b _13-22 For the 1 st-10 th perforation cluster position x of a specific fracturing well ₁₃ (xi), the position of each perforation cluster corresponds to the stratum parameter, and the non-perforation clusters when the cluster number is less than 10 clusters do not participate in multi-objective optimization; b _23-32 The diameter of the corresponding perforation hole of each perforation cluster position is determined, and the diameter of the perforation holes of the non-perforation clusters when the number of the clusters is less than 10 clusters does not participate in multi-objective optimization; b _33-42 Corresponding hole density for each perforation cluster position, wherein the non-perforation cluster hole density does not participate in multi-objective optimization when the cluster number is less than 10 clusters; b _43-52 Each perforation cluster position corresponds to the length of the perforation cluster, and the length of the non-perforation cluster is not involved when the cluster number is less than 10 clustersPerforming multi-objective optimization; each random sequence corresponds to an individual in the population.

Fitness calculation

And the fitness distribution adopts an improved fitness function to carry out assignment. Firstly, parent individuals are paired according to the principle of 'house-to-house', namely, the parents are sorted by fitness function (objective function) values, and the parent individuals are paired with small objective function values and large objective function values. And then determining the position of the intersection point by using the chaotic sequence, and finally intersecting the determined intersection item. In order to avoid the situation that the individuals in the non-dominant solution set and the dominant individuals in the population have the same fitness value, fitness assignment is respectively carried out on the individuals in the non-dominant solution set and the individuals in the population. The method comprises the following specific steps:

for each solution NP (i) e NP, an intensity value s (i) is assigned, representing the number of solution sets dominated:

in the formula: |. | is a cardinality representing a set; + represents a union set; is corresponding to Pareto dominance relation;

on the basis of the S value, calculating an original fitness r (i) of the individual i:

equation (109) indicates that the original fitness of the high-dimensional multi-objective parallel evolution optimization algorithm is determined by the individuals in the non-inferior solution set and the individuals in the population, and the genetic algorithm only considers the individuals in the population. In the process of evolutionary optimization, fitness needs to be minimized, i.e. r (i) ═ 0 corresponds to one non-dominant individual, while r (i) high value means that i is dominated by many species of population (and in turn dominates many species of population).

Original fitness assignment reflects the dominant and dominated information of the individuals, and k-order neighbor method density parameters are introduced to evaluate the individuals with the same fitness value.

In the formula:

the Euclidean distance between the individual i and the k-th individual after being sorted in ascending order.

Equation (110) adds 2 to the denominator to ensure that its value is greater than 0 and d (i) < 1. Then, adding d (i) to the original fitness value r (i) of the individual i to obtain its fitness f (i):

F(i)＝R(i)+D(i) (111)

environment variable selection

The selection of the environment variables is also a means for realizing the diversity of the group, and is an important guarantee for jumping out of the local optimum and searching the global optimum. Two differences exist between the environment variable selection process of the high-dimensional multi-target parallel evolution optimization algorithm and the genetic algorithm: the environment variable selection size is always a constant; the environment variable selection avoids boundary solutions from being removed from the population. The environment variable selection steps are as follows:

the method comprises the following steps: the population P _t And an external environment A _t Set replication all non-dominated solution set individuals evolve to

I.e., those with fitness below 1, are replicated to the next generation of environmental variables if

And if the size is equal to the size of the population N, the circulation is skipped.

Step two: if it is not

Then P will be _t And A _t The best is

An one is supportedAdding to the solution

In (1).

Step three: if it is not

Then it is removed from the environment according to the following principle, if the individual i satisfies the condition as in (112)

Removing until

Equation (112) indicates that i satisfies | P less than 0 for an individual _t+1 The Euclidean distance between the individual k and the individual i is equal to the Euclidean distance between the individual k and the individual j; or there is more than 0 and less than | P _t+1 If the Euclidean distance between the individual k and the individual i is less than the Euclidean distance between the individual k and the individual j, and if the Euclidean distance between the individual l and the individual i is less than the Euclidean distance between the individual l and the individual j, the individual i is selected from the group consisting of the individuals k, j, and j, f

And (5) removing.

7. Calculation examples and analysis

Taking a practical example as an example, based on the method disclosed in the application, intelligent regulation and control of unconventional reservoir fracturing multi-fracture non-equilibrium expansion are performed.

(1) Basic parameters

TABLE 1 basic parameters

According to 1462m fracture horizontal section vertical ground stress range (57-72MPa) of an X well, the minimum horizontal main stress range (50-70MPa), the maximum horizontal main stress (55-75MPa), the Young modulus (30000- ³ The viscosity of fracturing fluid (40 mPa & s at most), the diameter of an eyelet (8-16mm), the density of the eyelet (12-24 holes/m), the length of a perforation cluster (0.2-0.65m), and the design of an experimental scheme is developed by taking a 10 th fracturing section (3616 & 3546m) as an example; increasing the range of perforation positions (35-35 m corresponding to 3616-3546m and 30-30 m corresponding to 3362-3422 m). For vertical stress x ₁ Minimum horizontal principal stress x ₂ Maximum horizontal principal stress x ₃ Young's modulus x ₄ Poisson's ratio x ₅ Permeability x ₆ Porosity x ₇ Discharge x ₈ Viscosity x of fracturing fluid ₉ Diameter x of the hole ₁₀ Cell density x ₁₁ Length x of perforation cluster ₁₂ Perforation cluster position x ₁₃ And (3) adopting 13 independent variables as input parameters, adopting DOE (Design of Experiment), and calculating the crack length y of each Design point under the DOE test Design parameters according to the output parameters by using the output parameters and the non-planar crack propagation mechanism model established in the prior art ₁ Improved volume y ₂ Pressure y ₃ . Totally, 10812 groups of design points S (S) are designed ₁ ,s ₂ ,...,s ₁₀₈₁₂ ) Can reflect 17X 4 ¹² And (3) group testing, namely constructing a training sample big database X, and fully ensuring the precision of the training model by adopting a sample set. As shown in table 3, some DOE experimental design point input parameters and corresponding output parameters (where m is 10812 and q is 3) are shown, so as to create a mimo sample database covering the whole well section. And then based on the established unconventional reservoir multi-cluster perforation crack parameter multi-objective optimization algorithm, selecting a Sichuan basin unconventional reservoir gas well X well as a simulated parameter range, and carrying out unconventional reservoir multi-cluster perforation crack parameter multi-objective optimization by taking the 10 th fracturing segment (3616-3546m, figure 1) basic parameters as an example.

TABLE 3 training sample set

(2) Single target optimization analysis

Based on a sample big database, the established Gauss-Kriging agent model is used for training and predicting the seam length, the transformation volume and the construction pressure. The simulation basic parameters comprise that the vertical crustal stress is 67MPa, the minimum horizontal principal stress is 50MPa, the maximum horizontal principal stress is 60MPa, the Young modulus is 35000MPa, the Poisson ratio is 0.25, the permeability is 0.1mD, the porosity is 3 percent, and the construction displacement is 12.0m ³ The viscosity of the fracturing fluid is 5 mPas, the perforation aperture is 12.0mm, the perforation aperture is 20.0/m, the perforation position is 0m, and the number of the perforations is 0.

As can be seen in fig. 2, as the permeability increases from 0.01mD to 0.5mD and the porosity increases from 3% to 7%, the formation pore connectivity and permeability gradually increases to gradually dominate the increase in fracturing fluid loss, resulting in a decrease in propped hydraulic fracture fluid volume and a decrease in fracture length, the smaller the remodeling volume; the permeability of the reservoir is increased from 0.01 to 0.3mD, the injection speed of the fracturing fluid is less than the matching degree with the conductivity of the stratum, and the construction pressure is gradually reduced; when the permeability of the reservoir is increased from 0.3 to 0.5mD, the injection speed of the fracturing fluid is greater than the degree of matching with the conductivity of the stratum, and the construction pressure is gradually increased.

(3) Multi-objective evolutionary model verification and application

Based on stratum basic parameters of a 10 th fracturing section (3546-; the number of genetic iterations of the model is 15, each genetic iteration is 500, the genetic cross probability is 1, the mutation probability is the reciprocal of the lower bound, and the population size is 30.

The optimization process and the result of the invention are explained by the 10 th fracturing segment high-dimensional multi-objective parallel evolutionary optimization analysis

(ii) 2 Cluster optimization results

As can be seen from FIG. 3, as the number of different genetic iterations increases, the construction parameters of the 2-cluster optimization and the modification volume, construction pressure and the length variance of each cluster of fractures under perforation gradually reach the multi-objective optimization. The final optimization result is: the viscosity of the fracturing fluid is 19.7mPa & s; the 1 st to 2 nd cluster has the discharge capacity of 8m ³ /min、8m ³ Min; cluster 1-2-17.3 m, -27.8 m; the aperture of the 1 st-2 nd shower nozzle is 10.9mm and 12.4 mm; the density of the 1 st cluster of holes is 18 holes/m, and the density of the 2 nd cluster of holes is 17.0 holes/m; the length of the 1 st perforation cluster is 0.6m, and the length of the 1 st perforation cluster is 0.45 m; the length of the 1 st cluster of seams can reach 227.6m and the length of the 2 nd cluster of seams can reach 233.4m under the combination parameters; the reconstructed volume can reach 65.8 multiplied by 10 ⁴ m ³ (ii) a The construction pressure is 58.9 MPa.

6 clusters of optimized results

As can be seen from FIG. 4, as the number of different genetic iterations increases, the construction parameters of the 6-cluster optimization and the modification volume, construction pressure and the length variance of each cluster of fractures under perforation gradually reach the multi-objective optimization. The final optimization result is: the viscosity of the fracturing fluid is 5mPa & s; the 1 st to 6 th clusters have the discharge capacity of 3m in sequence ³ /min、3m ³ /min、3m ³ /min、2m ³ /min、2m ³ /min、3m ³ Min; the 1 st to 6 th shower holes are sequentially-22.75 m, -12.25m, -1.75m, 8.25m, 18.75m and 29.25 m; the aperture of the 1 st to 6 th shower holes is 10mm, 8mm, 10mm, 8mm and 10mm in sequence; the density of the 1 st to 6 th clusters of perforation holes is 16 holes/m, 20 holes/m, 16 holes/m, 20 holes/m and 16 holes/m in sequence; the lengths of the 1 st-6 th perforation clusters are 0.5m, 0.5m and 0.5m in sequence; the lengths of 1-6 clusters of seams are 92.8m, 83.8m, 60.4m, 70m, 85m and 76.5m in sequence under the combination parameters; the reconstructed volume can reach 128.6 multiplied by 10 ⁴ m ³ (ii) a The construction pressure is 64.99 MPa.

Thirdly, performing high-dimensional multi-target parallel evolution optimization results on different cluster numbers of the 10 th fracturing section

TABLE 3 high-dimensional multi-objective parallel evolutionary optimization results of different cluster numbers of the 10 th fracturing segment

According to the optimization result (table 3) of high-dimensional multi-objective parallel evolution, 6 clusters are comprehensively optimized by the large total modification volume and the small construction pressure, and the construction parameter and the perforation parameter are fracturing fluid viscosity of 5mPa & s; the 1 st to 6 th clusters have the discharge capacity of 3m in sequence ³ /min、3m ³ /min、3m ³ /min、2m ³ /min、2m ³ /min、3m ³ Min; the positions of the 1 st to 6 th shower holes are-22.75 m, -12.25m, -1.75m, 8.25m, 18.75m and 29.25m (3558.25m, 3568.75m, 3579.25m, 3589.25m, 3599.75m and 3610.25m) in sequence; the aperture of the 1 st to 6 th shower holes is 10mm, 8mm, 10mm, 8mm and 10mm in sequence; the density of the 1 st to 6 th clusters of perforation holes is 16 holes/m, 20 holes/m, 16 holes/m, 20 holes/m and 16 holes/m in sequence; the lengths of the 1 st-6 th perforation clusters are 0.5m, 0.5m and 0.5m in sequence; further analyzing the penetration rate of the perforation clusters corresponding to the perforation positions selected by the multi-objective optimization to be 0.018mD, 0.01mD, 0.015mD, 0.013mD, 0.008mD and 0.017mD, thereby showing that the high-permeability reservoir adopts small-aperture low-pore density (16 pores/m and 8mm), the low-permeability reservoir adopts large-aperture high-pore density (20 pores/m and 10mm), low viscosity (5-10mPa · s) and large displacement (16 m/m and 10mm), and the high-permeability reservoir adopts large-aperture high-pore density (20 pores/m and 10mm) ³ The matching measures of/min), cluster number (6-8 clusters) and cluster spacing (9-12m) can realize the balanced expansion of the multi-cluster perforation non-planar cracks of the target block, obtain larger reconstruction volume and have low construction pressure.

Claims

1. The method for regulating and controlling the highly-accurate intelligent fracturing of the unconventional reservoir is characterized by comprising the following steps of:

collecting basic parameters required by calculation;

establishing a planar multi-cluster crack initiation-expansion coupling model;

considering the flow of a non-planar fracture network fluid, coupling multi-fracture and bifurcation fracture flow dynamic distribution, and establishing a non-planar fracture initiation-expansion coupling model based on the interactive judgment criterion of hydraulic fractures and natural fractures considering anisotropic induced stress;

according to the change of a reservoir pore pressure field caused by the induction of the non-planar multi-cluster fractures, establishing a non-planar multi-cluster fracture reconstruction volume model based on reservoir rock shear slip and tensile failure mechanical conditions and a fracture network permeability model;

based on a non-planar crack initiation and propagation mechanism model, adopting DOE design to establish a sample database, training by using a Gauss-Kriging machine learning model, and establishing a multi-input multi-output intelligent agent model for simulating the mechanism model;

the method is characterized in that a multi-objective function is formed by taking the large modification volume and low construction pressure of each cluster of cracks matched with the geological characteristics of a reservoir stratum, combining constraint conditions and applying a genetic algorithm to comprehensively solve the optimal solution of the whole objective function, and performing high-precision intelligent fracturing regulation and control of unconventional reservoir stratum fracturing multi-crack unbalanced expansion.

2. The unconventional reservoir high-precision intelligent fracturing regulating and controlling method according to claim 1, wherein basic parameters required for calculation comprise reservoir parameters, mechanical parameters, fracturing parameters and wellbore parameters.

3. The unconventional reservoir high-precision intelligent fracturing regulation and control method of claim 1, wherein the establishing of the non-planar multi-cluster fracture reformation volume model further comprises:

acquiring the pore pressure after the pressure of the non-planar fracture is affected,

in the formula: p (x, y) is the pore pressure at plane (x, y) after the formation of the non-planar fracture pressure; p is a radical of _e Is the original formation pressure; p is a radical of _fi,j (x, y) is the pressure in the jth section of the ith perforation cluster fracture; p is a radical of _fi,j (x, y) is the fracture internal pressure at the plane (x, y) of the jth section of the ith natural fracture; m, s, u (i), v (k) represent numbers;

and (3) activating the modified volume of the natural fracture by utilizing the pore pressure after the non-planar fracture pressure is responded and combining with the opening shear failure criterion of the planar weak point and utilizing the multi-cluster fracture modified volume calculation process to stack and activate the modified volume of the natural fracture, and establishing a non-planar fracture modified volume simulation model.

4. The unconventional reservoir high-precision intelligent fracturing regulation and control method of claim 1, wherein the constraints comprise:

(1) the sum of the flow of each cluster of cracks is equal to the total construction displacement;

(2) the bottom hole pressure of each cluster is equal;

(3) the cluster spacing is greater than the set minimum cluster spacing;

(4) and the formation parameters, construction parameters and perforation parameters corresponding to each cluster are positioned between the minimum value and the maximum value set by the fracturing well.

5. The unconventional reservoir high-precision intelligent fracturing regulation and control method of claim 1, wherein the comprehensive solution of the overall objective function by applying the genetic algorithm further comprises: and (3) comprehensively optimizing one or more parameters of the number, the position, the cluster length, the pore density, the pore diameter, the construction discharge capacity and the viscosity of the fracturing fluid corresponding to the optimal solution of the whole objective function by applying an improved genetic algorithm.

6. The utility model provides a high accurate intelligent fracturing of unconventional reservoir regulates and control device which characterized in that includes following step:

the planar multi-cluster crack initiation-expansion coupling model module is used for establishing a planar multi-cluster crack initiation-expansion coupling model;

the non-planar fracture initiation-expansion coupling model is used for considering the flow of a non-planar fracture network fluid, coupling multi-fracture and bifurcation fracture flow dynamic distribution, and establishing the non-planar fracture initiation-expansion coupling model based on the interactive judgment criterion of hydraulic fractures and natural fractures considering anisotropic induced stress;

the non-planar multi-cluster fracture reconstruction volume prediction module is used for establishing a non-planar multi-cluster fracture reconstruction volume model based on the mechanical conditions of shear slip and tensile destruction of reservoir rock and a fracture network permeability model according to the change of a reservoir pore pressure field caused by the induction of the non-planar multi-cluster fracture;

the intelligent regulation and control optimized perforation and construction parameter acquisition module is used for establishing a sample database by adopting DOE (design of engineering) design based on a non-planar fracture initiation and propagation mechanism model, training by utilizing a Gauss-Kriging machine learning model and establishing a multi-input multi-output intelligent agent model for initiating the simulation of a mechanism model;

the high-precision intelligent fracturing regulation and control module comprehensively solves the optimal solution of the whole objective function by combining constraint conditions and applying a genetic algorithm by taking the large modification volume and the low construction pressure of each cluster of cracks matched with the geological characteristics of the reservoir as a multi-objective function, and carries out high-precision intelligent fracturing regulation and control of unconventional reservoir fracturing multi-crack unbalanced expansion.

7. The unconventional high-precision intelligent fracturing regulating and controlling device for reservoirs according to claim 6, wherein basic parameters required for calculation comprise reservoir parameters, mechanical parameters, fracturing parameters and wellbore parameters.

8. The unconventional reservoir high-precision intelligent fracturing regulation and control device of claim 6, wherein the establishing a non-planar multi-cluster fracture reformation volume model further comprises:

9. The unconventional reservoir high-precision intelligent fracturing regulating and controlling device of claim 6, wherein the constraints comprise:

(2) the bottom hole pressure of each cluster is equal;

(3) the cluster spacing is greater than the set minimum cluster spacing;

10. The unconventional reservoir high-precision intelligent fracturing regulation and control device of claim 6, wherein the comprehensive solution of the overall objective function by applying the genetic algorithm further comprises: and (3) comprehensively optimizing one or more parameters of the number, the position, the cluster length, the pore density, the pore diameter, the construction discharge capacity and the viscosity of the fracturing fluid corresponding to the optimal solution of the whole objective function by applying an improved genetic algorithm.