CN110863810B - Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process - Google Patents
Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process Download PDFInfo
- Publication number
- CN110863810B CN110863810B CN201911151497.6A CN201911151497A CN110863810B CN 110863810 B CN110863810 B CN 110863810B CN 201911151497 A CN201911151497 A CN 201911151497A CN 110863810 B CN110863810 B CN 110863810B
- Authority
- CN
- China
- Prior art keywords
- fracture
- crack
- natural
- stress
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 75
- 238000004519 manufacturing process Methods 0.000 title claims abstract description 52
- 238000004088 simulation Methods 0.000 title claims abstract description 40
- 230000008878 coupling Effects 0.000 title claims abstract description 17
- 238000010168 coupling process Methods 0.000 title claims abstract description 17
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 17
- 239000011159 matrix material Substances 0.000 claims abstract description 51
- 230000008569 process Effects 0.000 claims abstract description 35
- 239000012530 fluid Substances 0.000 claims abstract description 27
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 27
- 238000009826 distribution Methods 0.000 claims abstract description 16
- 238000011065 in-situ storage Methods 0.000 claims abstract description 11
- 230000015572 biosynthetic process Effects 0.000 claims abstract description 8
- 230000005514 two-phase flow Effects 0.000 claims abstract description 6
- 238000004364 calculation method Methods 0.000 claims description 17
- 230000035699 permeability Effects 0.000 claims description 15
- 238000002347 injection Methods 0.000 claims description 13
- 239000007924 injection Substances 0.000 claims description 13
- 238000004422 calculation algorithm Methods 0.000 claims description 12
- 239000011435 rock Substances 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 8
- 239000011541 reaction mixture Substances 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 4
- 230000003213 activating effect Effects 0.000 claims description 2
- 150000001875 compounds Chemical class 0.000 claims description 2
- 230000003044 adaptive effect Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 230000035945 sensitivity Effects 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000003776 cleavage reaction Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000007017 scission Effects 0.000 description 1
- 239000003079 shale oil Substances 0.000 description 1
- 239000000243 solution Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
- E21B43/25—Methods for stimulating production
- E21B43/26—Methods for stimulating production by forming crevices or fractures
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B49/00—Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Geology (AREA)
- Mining & Mineral Resources (AREA)
- Physics & Mathematics (AREA)
- Environmental & Geological Engineering (AREA)
- Fluid Mechanics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Geophysics (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
The invention discloses an integrated simulation method for coupling a shale gas reservoir hydraulic fracturing flowback production process, which comprises the following steps: (S1) setting a potential propagation path of the fracture based on the natural fracture distribution, in-situ stress magnitude and direction of the formation; (S2) taking all potential cracks as grid boundaries, and adopting triangular units to divide the space into unstructured grids; (S3) calculating the pressure and the saturation of the fractures and the matrix grid in the fracturing, flowback and production processes by adopting a gas-water two-phase discrete fracture model based on the unstructured grid; (S4) fracture propagation process simulation; (S5) judging the fracture intersection behavior when the fracture propagates forwards and meets the natural fracture; (S6) after the simulation of the fracturing and expanding process is finished, simulating the flowback and production process. The method provided by the invention combines a discrete fracture model and dynamic fractures, considers the two-phase flow of the fracturing fluid and shale gas in the stratum, and realizes the integrated simulation of the hydrofracture flowback production process of the shale gas reservoir.
Description
Technical Field
The invention relates to a shale gas reservoir hydraulic fracturing simulation method, in particular to an integrated simulation method for coupling a shale gas reservoir hydraulic fracturing flowback production process.
Background
Hydraulic fracturing has become an effective means for increasing production and improving the quality of oil and gas reservoirs, particularly shale oil and gas reservoirs, but for a long time, hydraulic fracturing simulation software mostly only considers the characteristic changes of rock fracture and crack expansion caused by the injection of a fracturing fluid into a stratum in the hydraulic fracturing process, does not consider the two-phase seepage relation between a fluid and the fracturing fluid in the stratum, and only considers the injection fracturing simulation of a single fluid, so that most fracturing simulations cannot reflect the influence of the injection of the fracturing fluid on a stratum matrix and oil and gas. If the interaction and the flow change between injected fracturing fluid and reservoir fluid in the fracturing process can be accurately described, the real formation pressure change and the permeation characteristic in the fracturing process are analyzed, and the coupling simulation is realized with the subsequent well closing and production processes, so that the method is not only beneficial to the research of the crack expansion characteristic in the fracturing process, but also has important significance for realizing the transformation production of oil and gas reservoirs and realizing the integrated optimization design of fracturing and production.
Disclosure of Invention
The invention aims to provide an integrated simulation method for a coupling shale gas reservoir hydraulic fracturing flowback production process, which solves the problem that the existing simulation method does not consider the two-phase seepage relation, combines a discrete fracture model and a dynamic fracture design, considers the two-phase flow of fracturing fluid and shale gas in a stratum and realizes the integrated simulation of the shale gas reservoir hydraulic fracturing-flowback-production process.
In order to achieve the above object, the present invention provides an integrated simulation method for coupling a shale gas reservoir hydraulic fracturing flowback production process, comprising:
(S1) setting all fractures to be perpendicular to the direction of the minimum horizontal principal stress based on the natural fracture distribution, the in-situ stress magnitude and direction of the formation, and setting propagation paths of the fractures, which are divided into: fracturing fractures, natural fractures;
(S2) taking all potential cracks as grid boundaries, adopting an unstructured grid based on triangular units to mesh the space, carrying out grid encryption on the peripheries of the cracks, enabling natural cracks and fracturing crack units to coincide with the adjacent matrix grid boundaries, and closing all crack units at the initial moment;
(S3) at the beginning of fracturing, setting the fracturing fracture unit connected with the perforation to be opened to the perforationThe water injection process of points replaces the injection process of fracturing fluid, and the opening degree w (x) and the permeability K of a fracture unit are calculated based on a gas-water two-phase discrete fracture model of a triangular unstructured gridfCalculating the conductivity T among the units based on the seepage characteristics of the fluid, and calculating the pressure and the saturation in the cracks and the matrix grids in the fracturing, flowback and production processes by adopting a conventional gas-water two-phase flow numerical simulation method after obtaining the conductivity T among various units;
(S4) simulation of crack propagation process: judging whether the crack unit is opened or not based on the pressure change characteristics and the stress intensity factor criterion of all potential crack units calculated by the gas-water two-phase discrete crack model: if the I-type stress concentration factor at the tip of the crack reaches a critical value, the crack begins to expand; when the crack begins to expand, activating a crack unit adjacent to the current crack, and simulating the crack expansion process;
(S5) when the fracture propagates forward to encounter the natural fracture, making a judgment on the fracture intersection behavior: judging whether the fracture passes through the met natural fracture or not based on the natural fracture attribute, the natural fracture orientation, the matrix mechanical attribute, the ground stress orientation and the size, and if so, judging whether the fracture passes through the met natural fracture along the critical radius r of the natural fracturecIf the circumferential stress borne by the matrix is greater than the tensile strength of the rock, new fracturing cracks are generated on the wall surfaces of the natural cracks; if at the critical radius rcThe position natural crack has slippage, the natural crack is damaged before a new crack is generated, and the fracturing crack cannot penetrate through the natural crack; if at the critical radius rcIf the natural crack does not slip, the fracturing crack penetrates through the natural crack and continues to expand along the original path;
(S6) after the simulation of the fracture propagation process is finished, closing the fractured well, and simulating the well-closing process by taking the stratum saturation and the pressure after fracturing as initial conditions; and after the well is closed, the water injection well is changed into a production well, and production is carried out in a constant pressure production mode so as to simulate a gas-water two-phase flowback process and a shale gas production process.
Preferably, in the step (S3), for the opened crack, with the crack center point as 0 point, the crack cell opening w along the crack length direction x is:
in the formula (1), p represents the fluid pressure in the fracture cell, σnThe in-situ normal stress acting on the fracture unit is shown, E and v respectively represent the Young modulus and Poisson ratio of the rock, and a represents the half-length of the fracture.
Preferably, in the step (S3), the permeability K of the fracture unit is based on a cubic criterionfComprises the following steps:
preferably, in the step (S3), the inter-cell conductivity T includes: matrix-matrix conductivity, matrix-fracture conductivity, and fracture-fracture conductivity.
Preferably, the matrix-matrix conductivity, for conductivity T between matrix meshes i, j of arbitrary shapei,jComprises the following steps:
in the formula (4), the reaction mixture is,
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the two matrix units, D is the distance from the central point of the matrix unit to the intersection surface,
is a normal unit vector of the intersecting surface,
subscripts i and j correspond to parameters of matrix grids i and j, which are unit vectors connecting the central points of the cells.
Preferably, the matrix-fracture conductivity is:
in the formula (11), the reaction mixture is,
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the grid unit, D is the distance from the center point of the grid unit to the intersection surface,
is a normal unit vector of the intersecting surface,
the unit vector is the unit vector connecting the center points of the cells, the subscript m represents the matrix property, and f represents the fracture property.
Preferably, the fracture-fracture conductivity, the conductivity between any two fracture grids i and j is:
in the formula (10), the compound represented by the formula (10),
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the crack grid unit, D is the distance from the center point of the crack grid unit to the intersection surface, and n represents the total number of all the intersected cracks at the point.
Preferably, in step (S4), the crack tip type I stress concentration factor is:
in the formula (11), KIRepresenting the I-type stress concentration factor of the crack tip, p representing the pressure intensity in the crack, sigmanIndicates the in situ normal stress acting on the fracture cell and a indicates the fracture half-length.
Preferably, the step size of the fracture propagationdt adopts an adaptive time step correction algorithm: fracture tip I-type stress concentration factor K at current time stepICritical stress concentration factor K associated with crack propagationICThe difference dK ═ KI-KICFor distinguishing the parameters, the absolute value of the difference is less than 0.1KICAdopting a dichotomy method to perform approximation calculation on the time step as the standard of stopping the algorithm;
if dK>0 and dK>0.1KICIf yes, updating the current time step dt to be the upper limit dt of the time stepmax(ii) a If dK<0 and dK>0.1KICIf yes, then update the current time step dt to the lower limit dt of the time stepmin(ii) a The updated time step is the average value dt of the upper limit and the lower limit of the time step is 0.5 (dt)max+dtmin) And continuously updating the upper limit and the lower limit of the time step until the algorithm stops, and obtaining the current time step for realizing the critical crack extension.
Preferably, in the step (S5), the critical radius rcComprises the following steps:
σ1(rc(θ))＝To(14)；
in formula (14), σ1Is maximum principal stress, ToThe tensile strength of the rock;
judging whether the crack slips according to a molar coulomb criterion:
|τβ|＜So-λσβn(15)；
in the formula (15), τβRepresents the critical distance rcTangential stress at (theta) along the natural fracture plane, S0Denotes the natural fracture cohesion, λ denotes the natural fracture friction coefficient, σβnRepresenting the normal stress acting on the natural fracture surface;
the stress distribution near the tip of the fracture can be calculated according to stress concentration factors, and is as follows:
in formula (16), σxx(θ)，σyy(theta) and tauxy(theta) represents the x and y directions of the crack tip, respectivelyNormal and tangential stress components of direction, KIAnd KIIRespectively representing stress concentration factors of the type I and the type II cracks, r represents the distance from any point in space to the crack tip, SHmaxAnd ShminRespectively representing the local maximum and minimum horizontal principal stress;
substituting formula (16) into formula (14), calculating critical length rc(β), β shows the angle between the current growth direction of the fracture and the natural fracture;
tangential stress τ of the natural fractureβStress σ normal toβnThe calculation formula is as follows:
the critical radius rc(β) bringing into formula (17) to obtain the tangential stress τ of the natural fractureβStress σ normal toβn；
If the tangential stress tau of the natural fractureβStress σ normal toβnIf equation (15) is not satisfied, the natural fracture begins to slip before the matrix is fractured, and the fractured fracture cannot penetrate through the natural fracture;
if the tangential stress tau of the natural fractureβStress σ normal toβnSatisfying equation (15), the fracture passes through the natural fracture.
The integrated simulation method for the coupling shale gas reservoir hydraulic fracturing flowback production process solves the problem that the existing simulation method does not consider the two-phase seepage relation, and has the following advantages:
(1) the method of the invention combines the design of the discrete fracture model and the dynamic fracture, solves the problem that the single-phase fluid fracturing simulation does not consider the two-phase seepage relation, considers the complexity of the fracture expansion, characterizes the fracture expansion characteristic based on the discrete fracture model and the dynamic fracture method, realizes the gas-water two-phase fracturing-flowback-production numerical simulation, considers the fracture expansion characteristic and the fluid change under the two-phase characteristic under the condition of formation fluid, can better reflect the real situation in the formation, and has important significance for knowing the formation fluid distribution in the fracturing process and improving the fracturing construction and flowback schemes;
(2) the method respectively calculates the conductivity of the seepage characteristics of the matrix and the cracks under different conditions, and is more accurate;
(3) according to the method, the crack propagation process is better simulated through judgment of different conditions according to the judgment of the intersection behavior of the crack and different crack propagation paths under different conditions;
(4) according to the method, in the crack expansion process, a self-adaptive time step correction algorithm is adopted, the set step length is ensured to be suitable, the difference between the stress concentration factor of the tip and the critical value after the current time step is finished is ensured to be small, and the requirement of critical expansion can be met.
Drawings
Fig. 1 is a flow chart of the integrated simulation method for the coupling shale gas reservoir hydraulic fracturing flowback production process.
FIG. 2 is a characteristic diagram of the fracture path and discrete fracture grid set by the present invention.
FIG. 3 is a graph illustrating the sensitivity of the present invention to guessing the initial time step without using the pressure length without adjusting the time step using the adaptive algorithm.
FIG. 4 is a fracture propagation signature after completion of hydraulic fracturing in accordance with the present invention.
FIG. 5 is a graph of matrix pressure at the end of fracturing, completion of stuffy well, and end of production according to the present invention.
FIG. 6 is a graph showing the pressure of the fracture at the end of fracturing, completion of well plugging and production.
FIG. 7 is a graph of the matrix saturation levels at the end of fracturing, completion of stuffy wells, and end of production according to the present invention.
FIG. 8 is a graph of fracture saturation levels corresponding to the end of fracturing, completion of stuffy wells, and end of production according to the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Example 1
An integrated simulation method for coupling a shale gas reservoir hydraulic fracturing flowback production process is shown in fig. 1, and is a flow chart of the integrated simulation method for coupling the shale gas reservoir hydraulic fracturing flowback production process, and comprises the following steps:
(S1) setting all the fracturing fractures to be vertical to the direction of the minimum horizontal principal stress based on the distribution of the natural fractures of the stratum and the magnitude and direction of the in-situ stress, wherein the fracturing fractures can pass through or not pass through the natural fractures, the possible propagation paths of all the fractures are preset under the setting, and the possible propagation paths are divided into two types of fracturing fractures and natural fractures, and as shown in FIG. 2, a characteristic diagram is drawn for the fracture paths and the discrete fracture grids set by the invention;
(S2) taking all potential cracks as grid boundaries, adopting an unstructured grid based on triangular units to mesh the space, carrying out grid encryption on the peripheries of the cracks, enabling natural cracks and fracturing crack units to coincide with the adjacent matrix grid boundaries, and closing all crack units at the initial moment;
(S3) when fracturing starts, setting the fracturing fracture units connected with the perforation to be opened, replacing the fracturing fluid injection process with the water injection process to the perforation point, and calculating the pressure and the saturation in the fractures and the matrix grids in the fracturing, flowback and production processes by adopting a gas-water two-phase discrete fracture model based on a triangular unstructured grid;
specifically, before calculating gas-water two-phase flow by using a discrete fracture model, calculation of the opening and conductivity of a fracture unit needs to be completed first, specifically as follows:
for the activated crack, the distribution of the opening degree w along the length direction x of the crack with the central point of the crack as a point 0 satisfies the following formula:
in the formula (1), p represents the fluid pressure in the fracture cell, σnThe in-situ normal stress acting on the fracture unit is shown, E and v respectively represent the Young modulus and Poisson ratio of the rock, and a represents the half-length of the fracture.
Permeability K of the fracture element based on cubic criterionfThe calculation formula is as follows:
after the permeability and the opening degree of the fracture unit are obtained, the conductivity between the fracture and the matrix unit in the discrete fracture model needs to be calculated, and the specific calculation mode is as follows:
for any control body, the flow formula between the units is as follows:
Qi,j＝Ti,jλ(pi-pj) (3)
in the formula (3), QijRepresents a flow rate from the lattice block i to the lattice block j in unit time; p represents a lattice block pressure; t isijConductivity, which is related to the properties of the grid and the porous medium only; l represents the fluid fluidity.
Conductivity calculations for the discrete fracture model are divided into: matrix-matrix (M-M) conductivity, matrix-fracture (M-F) conductivity, and fracture-fracture (F-F) conductivity.
The conductivity of the matrix-matrix (M-M) is calculated by taking the conductivity between matrix meshes i, j of arbitrary shape as a parameter independent of variables and related only to the geometry of the meshes and the properties of the porous medium, and is:
in the formula (4), the reaction mixture is,
a represents a grid cellK represents the harmonic mean of the permeabilities of the two matrix grid cells, D is the distance from the central point of the matrix grid cell to the intersection plane,
is a normal unit vector of the intersecting surface,
the subscripts i and j represent two matrix grid cells, respectively, for the unit vectors of the central points of the matrix grid cells connecting the lines.
The above calculation of the crack-crack (F-F) conductivity is not directly calculated by using a formula because there is generally no common plane between the two crack grids for the two connected crack units. To solve the problem, a virtual grid with infinitesimal volume, namely a '0 grid' is introduced for calculation, and the method is helpful for simulating the turning direction and the thickness change of the crack. The '0 grid' is orthogonally connected with two connected cracks, namely a connecting line 0-1 of the central points of the grids and a connecting line 0-2 of the central points of the grids are respectively orthogonal to the common plane of the grids. Assuming that all flow occurring from fracture grid 1 to fracture grid 2 is conducted through the "0 grid", then:
therefore, it is
In the formula (6), AiRepresenting the area of the intersection of the grid i and the grid 0, namely the width of the crack grid i; k is a radical ofiRepresenting the permeability of the fracture grid; diRepresenting the distance from the center point of the i-grid to the center point of the intersection.
Since the "0 grid" is set to infinitesimal, there are:
the conductivity between the connected fracture grids is then:
when the fractures intersect (including branch, merge, etc.), conductivity calculations between three or more intersecting fracture grids occur. Similar to the two connected fracture cells described above, a "0 grid" is also added at the intersection. Here, by analogy the flow problem of the fluid with the current problem, the fracture conductivity can be compared with the resistivity and transformed. From the star-to-corner transformation formula of kirchhoff's law, the star connection can be similarly converted into a triangular connection, so that:
in formula (9), the meanings of the variables are analogous to those of the above-mentioned connected cleavage units.
Similar to the star-to-delta transform, each of the above can be extended to the case where any n fracture grids intersect at the same node, where the conductivity formula between any two fracture grids i and j is as follows:
n represents the total number of all intersecting fractures at that point, and the amount of flow between the matrix and the fractures (M-F) can be expressed in terms of conductivity, in a manner similar to that for conductivity between the same media, although in the discrete fracture approach, the fractures are treated as a rectangle of thickness such that the fractures are perfectly orthogonal to the adjoining grid faces. The matrix-fracture conductivity calculation is essentially identical to the matrix-matrix conductivity calculation, as:
wherein the content of the first and second substances,
the subscript m represents the matrix properties and f represents the fracture properties.
After the conductivity among various units is obtained, the pressure and the saturation can be calculated according to a conventional gas-water two-phase flow numerical simulation method.
And (S4) judging whether the crack unit is opened or not based on the stress intensity factor criterion according to the pressure change characteristics in all the potential crack units calculated based on the discrete crack model. For any existing fracture, the fracture tip type I stress concentration factor can be calculated by:
in the formula (12), KIRepresenting the I-type stress concentration factor of the crack tip, p representing the pressure intensity in the crack, sigmanIndicates the in situ normal stress acting on the fracture cell and a indicates the fracture half-length. It is generally accepted that a fracture propagates in a critical state, i.e. at any one time:
KI＝KIC(13)
in formula (13), KICIs the critical stress concentration factor for crack propagation.
If the stress concentration factor reaches a critical value, the crack begins to propagate. If the crack begins to expand, a crack unit is activated forwards, and the crack expansion process is simulated. However, during actual calculation, the preset time step may be too short or too long, which results in too large difference between the stress concentration factor of the tip and the critical value after the current time step is finished, and thus the requirement of critical expansion is difficult to achieve. Therefore, the invention adopts a self-adaptive time step correction algorithm, and the difference dK between the stress concentration factor and the critical stress concentration factor under the current time step is equal to KI-KICFor distinguishing the parameters, the absolute value of the difference is less than 0.1KICAnd (3) adopting a dichotomy method to perform approximation calculation on the time step as a standard for stopping the algorithm. If dK is greater than 0 and the absolute value is greater than the tolerance, updating the current time step dt to be the upper limit dt of the time stepmax. Inverse directionIf dK is less than 0 and its absolute value is greater than the tolerance, then updating the current time step dt to be the lower limit dt of the time stepmin. The updated time step is the average value dt of the upper limit and the lower limit of the time step is 0.5 (dt)max+dtmin) Continuously updating the upper and lower limits of the time step until the algorithm stops, so as to obtain the current time step for realizing the critical crack extension, as shown in FIG. 3, which is a schematic diagram of the sensitivity of the pressure length of the unadjusted time step by the adaptive algorithm to the guess of the initial time step;
(S5) when the fracture propagates forward and encounters a natural fracture, a routing problem is encountered to determine whether the fracture propagates along or through the natural fracture, depending on the orientation of the natural fracture, the mechanical properties of the fracture and the matrix, and the in-situ stress magnitude and orientation of the formation. The specific judgment method is as follows:
according to the theory of linear elastic mechanics, the stress at the tip of the crack tends to be infinite, and the real rock is already yielding when the stress reaches a certain degree. The fracture process zone can be drawn by using the boundary that the rock stress reaches the critical tensile strength. The radius r of the process area of the crack tip at the theta angle position is determined according to the local coordinate system of the crack tipc(θ) is calculated as:
σ1(rc(θ))＝To(14)
in formula (14), σ1Maximum principal stress (pull positive), ToThe tensile strength of the rock.
If at the critical radius r along the natural fracturecAnd if the circumferential stress borne by the matrix is greater than the tensile strength of the rock, a new fracturing fracture is considered to be generated on the wall surface of the natural fracture. If the natural fracture has slipped at this critical radius, it is assumed that the fracture failure occurred before the new fracture was created and the fractured fracture cannot pass through the natural fracture. Otherwise, the fracture is considered to pass through the natural fracture and continue to propagate along the original path. Judging whether the crack slips according to a molar coulomb criterion:
|τβ|＜So-λσβn(15)
in the formula (15), τβRepresents the critical distance rcTangential stress at (theta) along the natural fracture plane, S0Denotes the natural fracture cohesion, λ denotes the natural fracture friction coefficient, σβnIndicating the normal stress acting on the natural fracture surface.
The stress distribution near the crack tip can be calculated according to the stress concentration factor, and the expression is as follows:
in formula (16), σxx(θ)，σyy(theta) and tauxy(theta) represents normal and tangential stress components of the fracture tip in the x and y directions, respectively, KIAnd KIIRespectively representing stress concentration factors of the type I and the type II cracks, r represents the distance from any point in space to the crack tip, SHmaxAnd ShminRespectively, the maximum and minimum horizontal principal stresses in place.
The critical length r can be calculated by substituting equation (16) for equation (14)c(β), β indicates the angle of the current growth direction of the fracture to the natural fracture, which is a known quantity.
The calculation formula of the tangential stress and the normal stress acting on the natural fracture is as follows:
calculating the critical radius rc(β) if the calculated stress does not satisfy equation (15), then the natural fracture has begun to slip before the matrix fractures and the fractured fracture cannot pass through the natural fracture.
(S6) after the simulation of the fracturing and expanding process is finished, closing the fracturing well, and simulating the well-closing process by taking the stratum saturation and the pressure after fracturing as initial conditions. And after the well is sealed, the water injection well is changed into a production well, and production is carried out in a constant pressure production mode so as to simulate a gas-water two-phase flowback process and a subsequent shale gas production process.
After fracturing is finished, the invention can obtain the path change of the fracture and the pressure distribution in the fracture in the fracturing process (see figure 4). From the calculation results, it can be seen that as the fracturing process progresses, the activated fracture unit increases and the pressure in the fracture unit gradually rises (1-4 in fig. 4). The present invention can also calculate the pressure distribution and water saturation distribution within the matrix (see a in fig. 5 and a in fig. 7) and the fracture unit (see a in fig. 6 and a in fig. 8) at the end of the fracture. And simulating the closed well and the flowback process by taking the pressure and saturation distribution after fracturing as initial conditions and taking the shut-in well and the transfer injection well as a constant pressure production well, and obtaining the pressure distribution and the saturation distribution in matrixes (see b and c in figure 5 and b and c in figure 7) and fracture units (see b and c in figure 6 and b and c in figure 8) at different time nodes.
As can be seen from FIGS. 5 and 6, at the end of fracturing, the pressure inside the fracture is higher, the pressure inside the fracture and its periphery after well closing is reduced but still higher than the initial pressure of the matrix, after the production process begins, the pressure inside the fracture is rapidly reduced, and at the end of production, the pressure inside the fracture is already significantly lower than the pressure of the periphery matrix. As can be seen from the water saturation profiles of fig. 7 and 8, at the end of the fracture, the water saturation in the fracture is close to 1, while the matrix water saturation is less due to slow fluid migration due to the short fracturing time and low matrix permeability. As the well-closing and production process progresses, the water saturation in the fracture is rapidly reduced, and the water saturation of the matrix is not obviously changed.
The shale gas reservoir hydraulic fracturing flowback production process integrated simulation method realizes numerical simulation of a gas-water two-phase fracturing-flowback-production process on the basis of representing fracture propagation characteristics by a discrete fracture model and a dynamic fracture method on the basis of considering the complexity of fracture propagation, can reflect the real situation in a stratum better by considering the fracture propagation characteristics and fluid changes under the two-phase characteristics under the condition of stratum fluid, and has important significance for knowing the distribution of the stratum fluid in the fracturing process and improving the fracturing construction and flowback schemes.
While the present invention has been described in detail with reference to the preferred embodiments, it should be understood that the above description should not be taken as limiting the invention. Various modifications and alterations to this invention will become apparent to those skilled in the art upon reading the foregoing description. Accordingly, the scope of the invention should be determined from the following claims.
Claims (8)
1. A method for integrally simulating a coupling shale gas reservoir hydraulic fracturing flowback production process is characterized by comprising the following steps of:
(S1) setting all fractures to be perpendicular to the direction of the minimum horizontal principal stress based on the natural fracture distribution, the in-situ stress magnitude and direction of the formation, and setting propagation paths of the fractures, which are divided into: fracturing fractures, natural fractures;
(S2) taking all potential cracks as grid boundaries, adopting an unstructured grid based on triangular units to mesh the space, carrying out grid encryption on the peripheries of the cracks, enabling natural cracks and fracturing crack units to coincide with the adjacent matrix grid boundaries, and closing all crack units at the initial moment;
(S3) when fracturing starts, setting the fracturing fracture unit connected with the perforation to be opened, replacing the fracturing fluid injection process with the water injection process to the perforation point, and calculating the opening degree w (x) and the permeability K of the fracture unit based on the gas-water two-phase discrete fracture model of the triangular unstructured gridfCalculating the conductivity T among the units based on the seepage characteristics of the fluid, and calculating the pressure and the saturation in the cracks and the matrix grids in the fracturing, flowback and production processes by adopting a conventional gas-water two-phase flow numerical simulation method after obtaining the conductivity T among various units;
(S4) simulation of crack propagation process: judging whether the crack unit is opened or not based on the pressure change characteristics and the stress intensity factor criterion of all potential crack units calculated by the gas-water two-phase discrete crack model: if the I-type stress concentration factor at the tip of the crack reaches a critical value, the crack begins to expand; when the crack begins to expand, activating a crack unit adjacent to the current crack, and simulating the crack expansion process;
(S5) when the fracture propagates forward to encounter the natural fracture, making a judgment on the fracture intersection behavior: based on natural fracture properties, orientationDetermining whether the fracture passes through the encountered natural fracture or not according to the mechanical property of the matrix, the ground stress direction and the size of the ground stress, and if the fracture passes through the encountered natural fracture, determining that the fracture passes through the encountered natural fracture along the critical radius r of the natural fracturecIf the circumferential stress borne by the matrix is greater than the tensile strength of the rock, new fracturing cracks are generated on the wall surfaces of the natural cracks; if at the critical radius rcThe position natural crack has slippage, the natural crack is damaged before a new crack is generated, and the fracturing crack cannot penetrate through the natural crack; if at the critical radius rcIf the natural crack does not slip, the fracturing crack penetrates through the natural crack and continues to expand along the original path;
(S6) after the simulation of the fracture propagation process is finished, closing the fractured well, and simulating the well-closing process by taking the stratum saturation and the pressure after fracturing as initial conditions; after the well is sealed, the water injection well is changed into a production well, production is carried out in a constant pressure production mode, and the method is used for simulating a gas-water two-phase flowback process and a shale gas production process;
in step (S3), for an opened fracture, the fracture cell opening w along the fracture length direction x with the fracture center point as 0 point is:
in the formula (1), p represents the fluid pressure in the fracture cell, σnThe in-situ normal stress acting on the fracture unit is shown, E and v respectively represent the Young modulus and the Poisson ratio of the rock, and a represents the half length of the fracture;
in step (S3), based on the cubic criterion, the permeability K of the fracture unitfComprises the following steps:
2. the integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process of claim 1, wherein in step (S3), the intercell conductivity T comprises: matrix-matrix conductivity, matrix-fracture conductivity, and fracture-fracture conductivity.
3. The integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process as claimed in claim 2, wherein the matrix-matrix conductivity is the conductivity T between matrix grids i, j of any shapei,jComprises the following steps:
in the formula (4), the reaction mixture is,
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the two matrix units, D is the distance from the central point of the matrix unit to the intersection surface,
is a normal unit vector of the intersecting surface,
subscripts i and j correspond to parameters of matrix grids i and j, which are unit vectors connecting the central points of the cells.
4. The integrated simulation method for the coupling shale gas reservoir hydraulic fracturing flowback production process of claim 2, wherein the matrix-fracture conductivity is:
in the formula (11), the reaction mixture is,
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the grid unit, D is the distance from the center point of the grid unit to the intersection surface,
is a normal unit vector of the intersecting surface,
the unit vector is the unit vector connecting the center points of the cells, the subscript m represents the matrix property, and f represents the fracture property.
5. The integrated simulation method for the coupling shale gas reservoir hydraulic fracturing flowback production process according to claim 2, wherein the fracture-fracture conductivity, the conductivity between any two fracture grids i and j is:
in the formula (10), the compound represented by the formula (10),
a represents the intersection surface of the grid unit, K represents the harmonic average value of the permeability of the crack grid unit, D is the distance from the center point of the crack grid unit to the intersection surface, and n represents the total number of all the intersected cracks at the point.
6. The integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process of claim 1, wherein in step (S4), the fracture tip type I stress concentration factor is:
in the formula (12), KIRepresenting the I-type stress concentration factor of the crack tip, p representing the pressure intensity in the crack, sigmanIndicates the in situ normal stress acting on the fracture cell and a indicates the fracture half-length.
7. The integrated simulation method for coupled shale gas reservoir hydraulic fracturing flowback production process as claimed in any one of claims 1-6, wherein the integrated simulation method isCharacterized in that the step dt of the crack propagation adopts a self-adaptive time step correction algorithm: fracture tip I-type stress concentration factor K at current time stepICritical stress concentration factor K associated with crack propagationICThe difference dK ═ KI-KICFor distinguishing the parameters, the absolute value of the difference is less than 0.1KICAdopting a dichotomy method to perform approximation calculation on the time step as the standard of stopping the algorithm;
if dK>0 and | dK ∞>0.1KICIf yes, updating the current time step dt to be the upper limit dt of the time stepmax(ii) a If dK<0 and | dK ∞>0.1KICIf yes, then update the current time step dt to the lower limit dt of the time stepmin(ii) a The updated time step is the average value dt of the upper limit and the lower limit of the time step is 0.5 (dt)max+dtmin) And continuously updating the upper limit and the lower limit of the time step until the algorithm stops, and obtaining the current time step for realizing the critical crack extension.
8. The integrated simulation method for coupled shale gas reservoir hydraulic fracturing flowback production process as claimed in claim 1, wherein in step (S5), the critical radius rcComprises the following steps:
σ1(rc(θ))＝To(14)；
in formula (14), σ1Is maximum principal stress, ToThe tensile strength of the rock;
judging whether the crack slips according to a molar coulomb criterion:
|τβ|＜So-λσβn(15)；
in the formula (15), τβRepresents the critical distance rcTangential stress at (theta) along the natural fracture plane, S0Denotes the natural fracture cohesion, λ denotes the natural fracture friction coefficient, σβnRepresenting the normal stress acting on the natural fracture surface;
the stress distribution near the tip of the fracture can be calculated according to stress concentration factors, and is as follows:
in formula (16), σxx(θ)，σyy(theta) and tauxy(theta) represents normal and tangential stress components of the fracture tip in the x and y directions, respectively, KIAnd KIIRespectively representing stress concentration factors of the type I and the type II cracks, r represents the distance from any point in space to the crack tip, SHmaxAnd ShminRespectively representing the local maximum and minimum horizontal principal stress;
substituting formula (16) into formula (14), calculating critical length rc(β), β shows the angle between the current growth direction of the fracture and the natural fracture;
tangential stress τ of the natural fractureβStress σ normal toβnThe calculation formula is as follows:
the critical radius rc(β) bringing into formula (17) to obtain the tangential stress τ of the natural fractureβStress σ normal toβn；
If the tangential stress tau of the natural fractureβStress σ normal toβnIf equation (15) is not satisfied, the natural fracture begins to slip before the matrix is fractured, and the fractured fracture cannot penetrate through the natural fracture;
if the tangential stress tau of the natural fractureβStress σ normal toβnSatisfying equation (15), the fracture passes through the natural fracture.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911151497.6A CN110863810B (en) | 2019-11-21 | 2019-11-21 | Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911151497.6A CN110863810B (en) | 2019-11-21 | 2019-11-21 | Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110863810A CN110863810A (en) | 2020-03-06 |
| CN110863810B true CN110863810B (en) | 2020-08-18 |
Family
ID=69655142
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911151497.6A Active CN110863810B (en) | 2019-11-21 | 2019-11-21 | Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110863810B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111879674B (en) * | 2020-07-15 | 2022-02-01 | 西南石油大学 | Testing device and method for determining reasonable well closing time based on shale imbibition permeability |
| CN113111607B (en) * | 2021-04-15 | 2022-04-15 | 西南石油大学 | Oil reservoir flowing full-coupling pressure production integrated numerical simulation method |
| CN114547998B (en) * | 2022-02-28 | 2024-03-22 | 西南石油大学 | Method for determining fracturing modification volume of horizontal well through coupling reservoir fluid |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104863560A (en) * | 2015-03-09 | 2015-08-26 | 东方宝麟科技发展(北京)有限公司 | Wide-net fracturing method for shale gas exploitation |
| CN107622165A (en) * | 2017-09-25 | 2018-01-23 | 西南石油大学 | A kind of shale gas horizontal well refracturing Productivity |
| CN109408859A (en) * | 2018-09-05 | 2019-03-01 | 中国石油集团川庆钻探工程有限公司 | Shale gas reservoir pressure break horizontal well two dimension treble medium numerical model method for building up |
| CN109472037A (en) * | 2017-09-08 | 2019-03-15 | 中国石油化工股份有限公司 | Shale gas reservoir man-made fracture parameter preferred method and system |
| CN109505576A (en) * | 2017-09-13 | 2019-03-22 | 中国石油化工股份有限公司 | Shale hydraulic fracturing Three-dimensional full coupling discrete fracture network analogy method and system |
| CN110206522A (en) * | 2019-06-10 | 2019-09-06 | 西南石油大学 | A kind of shale gas reservoir pressure break horizontal well fracturing fluid recovery (backflow) analogy method |
| CN110222477A (en) * | 2019-07-08 | 2019-09-10 | 西南石油大学 | The perforating parameter optimization method of maintenance level well staged fracturing crack equilibrium extension |
| CN110219630A (en) * | 2019-06-04 | 2019-09-10 | 西南石油大学 | A kind of fracturing fluid recovery calculation method of fractured sandstone gas reservoir pressure break horizontal well |
-
2019
- 2019-11-21 CN CN201911151497.6A patent/CN110863810B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104863560A (en) * | 2015-03-09 | 2015-08-26 | 东方宝麟科技发展(北京)有限公司 | Wide-net fracturing method for shale gas exploitation |
| CN109472037A (en) * | 2017-09-08 | 2019-03-15 | 中国石油化工股份有限公司 | Shale gas reservoir man-made fracture parameter preferred method and system |
| CN109505576A (en) * | 2017-09-13 | 2019-03-22 | 中国石油化工股份有限公司 | Shale hydraulic fracturing Three-dimensional full coupling discrete fracture network analogy method and system |
| CN107622165A (en) * | 2017-09-25 | 2018-01-23 | 西南石油大学 | A kind of shale gas horizontal well refracturing Productivity |
| CN109408859A (en) * | 2018-09-05 | 2019-03-01 | 中国石油集团川庆钻探工程有限公司 | Shale gas reservoir pressure break horizontal well two dimension treble medium numerical model method for building up |
| CN110219630A (en) * | 2019-06-04 | 2019-09-10 | 西南石油大学 | A kind of fracturing fluid recovery calculation method of fractured sandstone gas reservoir pressure break horizontal well |
| CN110206522A (en) * | 2019-06-10 | 2019-09-06 | 西南石油大学 | A kind of shale gas reservoir pressure break horizontal well fracturing fluid recovery (backflow) analogy method |
| CN110222477A (en) * | 2019-07-08 | 2019-09-10 | 西南石油大学 | The perforating parameter optimization method of maintenance level well staged fracturing crack equilibrium extension |
Non-Patent Citations (3)
| Title |
|---|
| Modelling fluid flow through a single fracture with different contacts using cellular automata;Peng-ZhiPand等;《Computers and Geotechnics》;20111231;第959-969页 * |
| 多尺度页岩气藏水平井压裂产能模拟研究;舒亮;《中国优秀硕士学位论文全文数据库工程科技I辑》;20150930;第22-39页 * |
| 页岩水力压裂的关键力学问题;庄茁等;《科学通报》;20160131;第61卷(第1期);第72-81页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110863810A (en) | 2020-03-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110863810B (en) | Integrated simulation method for coupling shale gas reservoir hydraulic fracturing flowback production process | |
| CN107545113B (en) | Method for simulating formation process of complex fracture network of hydraulic fracturing of unconventional oil and gas reservoir | |
| CN109505576B (en) | Shale hydraulic fracturing three-dimensional full-coupling discrete fracture network simulation method and system | |
| CN112036098A (en) | Method for simulating hydraulic fracture propagation numerical value of deep oil and gas reservoir | |
| CN108681635B (en) | Compact reservoir volume fracturing compressibility evaluation method | |
| CN110017135B (en) | Method for predicting crack propagation pressure of well wall of fractured stratum | |
| CN112012712B (en) | Numerical simulation method and device for water injection growth seam of embedded discrete seam | |
| Segura et al. | Coupled HM analysis using zero‐thickness interface elements with double nodes—Part II: Verification and application | |
| CN104695928A (en) | Method for evaluating volume transformation capacity of vertical well of fractured tight oil reservoir | |
| CN113076676B (en) | Unconventional oil and gas reservoir horizontal well fracture network expansion and production dynamic coupling method | |
| Hustedt et al. | Induced fracturing in reservoir simulations: Application of a new coupled simulator to a waterflooding field example | |
| CN110348031B (en) | Numerical simulation method for horizontal well fracturing near-wellbore fracture distortion form | |
| CN111577269B (en) | Multi-cluster fracturing fracture morphology prediction method based on discrete element fluid-solid coupling | |
| Liu et al. | Coupled flow/geomechanics modeling of interfracture water injection to enhance oil recovery in tight reservoirs | |
| CN116306385B (en) | Oil reservoir fracturing imbibition energy increasing numerical simulation method, system, equipment and medium | |
| CN112541287A (en) | Loose sandstone fracturing filling sand control production increase and profile control integrated design method | |
| CN115510778A (en) | Continental facies shale reservoir infinite stage fracturing process optimization method and system | |
| CN111125905B (en) | Two-dimensional fracture network expansion model for coupling oil reservoir fluid flow and simulation method thereof | |
| Zhang et al. | Phase field model for simulating hydraulic fracture propagation and oil–water two-phase flow in reservoir | |
| Miao et al. | An easy and fast EDFM method for production simulation in shale reservoirs with complex fracture geometry | |
| Wang et al. | Failure patterns and mechanisms of hydraulic fracture propagation behavior in the presence of naturally cemented fractures | |
| CN115288650A (en) | Method for parallel computing and simulating hydraulic fracturing in pore elastic medium | |
| CN113919201A (en) | Multi-scale expansion grid self-adaption method for hydraulic fracturing fracture | |
| Han et al. | An unequal fracturing stage spacing optimization model for hydraulic fracturing that considers cementing interface integrity | |
| CN112360448B (en) | Method for determining post-pressure soaking time by utilizing hydraulic fracture creep expansion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |