CN115034161B - Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration - Google Patents
Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration Download PDFInfo
- Publication number
- CN115034161B CN115034161B CN202210771868.6A CN202210771868A CN115034161B CN 115034161 B CN115034161 B CN 115034161B CN 202210771868 A CN202210771868 A CN 202210771868A CN 115034161 B CN115034161 B CN 115034161B
- Authority
- CN
- China
- Prior art keywords
- hydraulic fracture
- calculation
- time step
- grid
- tip
- 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
- 238000004364 calculation method Methods 0.000 title claims abstract description 90
- 230000001133 acceleration Effects 0.000 title claims abstract description 11
- 239000011159 matrix material Substances 0.000 claims abstract description 60
- 239000011435 rock Substances 0.000 claims abstract description 60
- 239000012530 fluid Substances 0.000 claims abstract description 44
- 230000035699 permeability Effects 0.000 claims abstract description 40
- 238000002347 injection Methods 0.000 claims abstract description 25
- 239000007924 injection Substances 0.000 claims abstract description 25
- 238000000034 method Methods 0.000 claims abstract description 11
- 230000000087 stabilizing effect Effects 0.000 claims abstract description 9
- 230000003044 adaptive effect Effects 0.000 claims description 18
- 238000006073 displacement reaction Methods 0.000 claims description 5
- 238000004088 simulation Methods 0.000 description 4
- 230000004048 modification Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 239000002253 acid Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000000243 solution Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computing Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
The invention discloses a self-adaptive time step calculation method for stabilizing three-dimensional hydraulic fracture expansion calculation and acceleration, which comprises the following steps: calculating the size of a cutting area and the number of reserved grids in each direction of a matrix rock mass area, and determining a corresponding grid reserved interval; calculating the width of the expanded grid and the equivalent permeability in the hydraulic fracture area; calculating the fluid pressure distribution in four time steps of future injection; calculating hydraulic fracture expansion critical pressure; judging whether the pre-selected time step meets the requirement or not, if not, re-calculating the pre-selected time step until the pre-selected time step meets the requirement; determining whether the next time step requires recalculation of the time step; the self-adaptive time step calculation accelerated under the condition of guaranteeing the stability of the three-dimensional hydraulic fracture expansion calculation can be realized by repeating the steps. According to the hydraulic fracture expansion method, based on a hydraulic fracture expansion criterion, the hydraulic fracture tip pressure is used as a judgment basis, and the calculation time steps are corrected in real time in each calculation step, so that the time step value is the maximum time step length capable of accurately calculating the hydraulic fracture expansion.
Description
Technical Field
The invention relates to the technical field of three-dimensional hydraulic fracture expansion numerical simulation in the hydraulic fracturing and acidizing fracturing processes, in particular to a self-adaptive time step calculation method for stabilizing and accelerating the three-dimensional hydraulic fracture expansion calculation.
Background
With the progress of exploration and development, hydraulic fracturing and acidizing fracturing have become indispensable technical means in the field of oil and gas exploration. Because of the complexity of the oil and gas reservoir conditions, the hydraulic fracture expansion numerical simulation technology is an important means for hydraulic fracturing and acidizing fracturing optimization design. With the progress of technology, a three-dimensional hydraulic fracture propagation model has become a main development direction of hydraulic fracture propagation numerical simulation. Meanwhile, in the hydraulic fracture expansion process, injected fluid is subjected to massive fluid loss by multi-scale flowing media such as hydraulic fracture wall surfaces, natural fractures and the like, so that the fluid fracture efficiency is reduced, and accurate calculation of the fluid loss is important to hydraulic fracture expansion numerical simulation. Early studies showed that in order to ensure stability and accuracy of hydraulic fracture expansion calculation, a very small time step (typically 0.01-0.1 s) is required at the initial stage of expansion, and the hydraulic fracture expansion calculation can be properly relaxed at the later stage of expansion. The whole hydraulic fracturing and acid fracturing construction time is often counted in hours or even days, and under the condition that a model matrix domain is a three-dimensional model, the small time step can lead to very high calculation cost in the calculation process of the crack width and the fluid filtration loss. Meanwhile, in order to ensure the conservation of mass of injected fluid in the crack propagation process, iterative calculation needs to be further carried out on the flow field in the crack, the matrix seepage field, the fluid loss and the width of the crack (influenced by the calculation of the pressure of the flow field in the crack). The above process makes the calculation cost too high, and it is often difficult (or takes a very long time) for the ordinary calculation resources to implement the above process.
Disclosure of Invention
In order to overcome the problems in the prior art, the invention provides a self-adaptive time step calculation method for stabilizing and accelerating the calculation of the three-dimensional hydraulic fracture expansion.
The technical scheme provided by the invention for solving the technical problems is as follows: the self-adaptive time step calculation method for calculating and accelerating the expansion of the stable three-dimensional hydraulic fracture comprises the following steps:
Step S1, acquiring a hydraulic fracture area omega hf, a hydraulic fracture length L hf, a hydraulic fracture height H hf and a hydraulic fracture width distribution w (x, z) when the current time step number t i is acquired;
S2, calculating the size d cut of a cutting area and the reserved grid number of the matrix rock mass area in each direction, determining a corresponding grid reserved interval, and cutting the matrix rock mass area of the three-dimensional hydraulic fracture expansion model;
s3, after cutting of the matrix rock mass area is completed, expanding a plane where the hydraulic fracture is located into a grid, and calculating the width of the expanded grid and the equivalent permeability in the hydraulic fracture area;
s4, selecting a preselected time step t ini according to the average permeability of the reservoir matrix;
step S5, calculating fluid pressure distribution P hf (x, z) of the hydraulic fracture in four time steps of future injection according to a preselected time step t ini;
Step S6, finding a grid corresponding to the maximum fluid pressure in the grid at the tip of the hydraulic fracture, defining the grid pressure as P tip,max, and injecting four time steps into the grid in the future to obtain calculation results of the maximum fluid pressure at the tip of the fracture, wherein the calculation results are P tip,max(1)、Ptip,max(2)、Ptip,max(3)、Ptip,max (4);
Step S7, calculating a hydraulic fracture expansion critical width, and calculating a corresponding hydraulic fracture expansion critical pressure P ic according to the hydraulic fracture expansion critical width;
Step S8, judging whether the preselected time step meets the requirements according to the hydraulic fracture expansion critical pressures P ic and P tip,max(3)、Ptip,max (4), if not, calculating a new preselected time step, and repeating the steps S5-S8 until the preselected time step meets the requirements, and then carrying out the next step;
Step S9, determining whether the next time step needs to recalculate the self-adaptive time step according to the hydraulic fracture area increment caused by the hydraulic fracture expansion in the current time step; repeating the steps S1-S8 to obtain a new self-adaptive time step if the self-adaptive time step needs to be recalculated;
And step S10, repeating the steps S1-S9, so that the accelerated self-adaptive time step calculation under the condition of guaranteeing the stability of the three-dimensional hydraulic fracture expansion calculation can be realized.
The further technical scheme is that the calculation formula of the cutting area size d cut in the step S2 is as follows:
Wherein: d cut is the size of the cutting area, m; a hf is the hydraulic fracture area, m 2;km is the average permeability of the reservoir matrix, m 2,Phf (x, z) is the fluid pressure in the hydraulic fracture area Ω hf, pa; p e is reservoir fluid pressure, pa; mu is the viscosity of the injected fluid and Pa.s; q is the injection displacement, m 3/s;tmax is the upper limit of the adaptive time step range, s.
The further technical scheme is that in the S2, the calculation formula of reserving the grid number in each direction by taking the injection point as the origin in the matrix rock mass area is as follows:
Rx＝Π(dcut/Δx)
Rz＝∏(dcut/Δz)
Wherein: d cut is the size of the cutting area, m; r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δy is the grid length in the y direction, m, where in order to ensure the calculation accuracy, the y direction generally adopts a logarithmic grid; Δz is the grid length in the z direction, m;
If the injection point (hydraulic fracture center point) is taken as an origin, and the grid position number of the injection point is set to be (0, 0), the grid reservation intervals in different directions after cutting are as follows:
x direction:
[-∏(Lhf/2Δx)-Rx,∏(Lhf/2Δx)+Rx]
y direction:
[-Ry,Ry]
z direction:
[Π(Hhf/2Δz)-Rz,∏(Hhf/2Δz)+Rz]
Wherein: r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; l hf is the hydraulic fracture length, m; h hf is the hydraulic fracture height, m.
The further technical scheme is that the calculation formula of the expanded grid width in the step S3 is as follows:
wherein: y kc is the expanded grid width, m; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; a hf is the hydraulic fracture area, m 2; w (x, z) is the hydraulic fracture width distribution, m.
The further technical scheme is that in the step S3, a calculation formula of the equivalent permeability in the hydraulic fracture area is as follows:
Wherein: y kc is the expanded grid width, m; k x,hf (x, z) is the equivalent permeability in the x direction in the hydraulic fracture region Ω hf, and m 2;ky,hf (x, z) is the equivalent permeability in the y direction in the hydraulic fracture region Ω hf, m 2; w (x, z) is the hydraulic fracture width distribution, m.
The further technical scheme is that the selection of the pre-selection time step t ini in the step S4 is as follows:
Wherein: k m is the reservoir matrix average permeability; t ini is a preselected time step, s.
The further technical scheme is that the calculation formula of the fluid pressure distribution in the step S5 is as follows:
Wherein: mu is the viscosity of the injected fluid and Pa.s; k x is the model x-direction permeability, m 2;ky is the model y-direction permeability, m 2;kz is the model z-direction permeability, m 2; phi is porosity, dimensionless; p m is the fluid pressure, pa; t is the calculation time, s; q f is the source term due to injection, s -1.
The further technical scheme is that in the step S7, a calculation formula of the critical width of the hydraulic fracture propagation is as follows:
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; k IC is fracture toughness of reservoir rock, pa.m 1/2;wtip is critical width of hydraulic fracture propagation, and m; Δx is the grid length in the x direction, m;
The calculation formula of the hydraulic fracture expansion critical pressure P ic is as follows:
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; w tip is the critical width of hydraulic fracture propagation, m; p ic is hydraulic fracture propagation critical pressure, pa; Δz is the grid length in the z direction, m; σ oh is the minimum horizontal principal stress, pa.
According to a further technical scheme, if the relation between P ic and P tip,max(3)、Ptip,max (4) in the step S8 meets the following formula, the current pre-selected time step is judged to meet the requirement, and t ini is the self-adaptive time step obtained through calculation:
Ptip,max(3)≤Pic≤Ptip,max(4)
If not, a new preselected time step is calculated based on the relative relationship of P tip,max (3) or P tip,max (4) to P ic,
If P ic≤Ptip,max (3), then take:
if P tip,max(4)＜Pic, then take:
Wherein: p ic is hydraulic fracture propagation critical pressure, pa; t ini is a preselected time step, s; t' ini is a new preselected time step, s.
The further technical scheme is that the condition of recalculating the self-adaptive time step in the step S9 is that the area change of the hydraulic fracture caused by the expansion of the hydraulic fracture meets the following conditions:
Ahf(ti)/Ahf(tn)≥1.1
Wherein: a hf is the hydraulic fracture area, and m 2;ti is the current time step number; t n is the number of time steps to obtain the current adaptive time step.
The invention has the following beneficial effects: according to the hydraulic fracture expansion method, based on a hydraulic fracture expansion criterion, the hydraulic fracture tip pressure is used as a judgment basis, and the calculation time step is corrected in real time in a required calculation step, so that the time step value is the maximum time step length for accurately calculating the hydraulic fracture expansion, the stability of hydraulic fracture expansion calculation is ensured, and the calculation acceleration is realized.
Drawings
FIG. 1 is a schematic diagram of a three-dimensional hydraulic fracture propagation model;
FIG. 2is a schematic view of a matrix rock mass region cut;
FIG. 3 is a schematic view of a fracture plane expanded into a grid;
FIG. 4 is a hydraulic fracture width distribution plot;
FIG. 5 is a hydraulic fracture zone equivalent permeability profile;
FIG. 6 is a graph of fluid pressure distribution in a hydraulic fracture at a fourth time step;
fig. 7 is a graph of the result of the adaptive time-step calculation.
Detailed Description
The following description of the embodiments of the present invention will be made apparent and fully in view of the accompanying drawings, in which some, but not all embodiments of the invention are shown. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
It is assumed that hydraulic fracture propagation is calculated in a three-dimensional matrix rock mass region of the size l×h×w and couples fluid flow within the hydraulic fracture, fluid seepage within the matrix rock mass and fluid loss from the hydraulic fracture to the matrix rock mass. Its main control equations include:
① Hydraulic fracture internal flow control equation:
② Matrix rock internal flow control equation:
Wherein: w is the width of the hydraulic fracture, m; p hf is the fluid pressure in the hydraulic fracture, pa; v nf is fluid loss velocity in hydraulic fracture, m/s; mu is the flow velocity of the fluid in the y direction, m/s; t is acid injection time, s; c t is the comprehensive compression coefficient of the oil reservoir, pa -1;kx is the permeability in the x direction of the model, m 2;ky is the permeability in the y direction of the model, m 2;kz is the permeability in the z direction of the model, and m 2; phi is porosity, dimensionless; p m is the fluid pressure, pa; q f is the source term due to injection, s -1.
③ The control equation for fluid loss from hydraulic fracture to matrix rock mass:
Wherein: v hf is the fluid loss rate from the hydraulic fracture to the matrix rock mass, m/s; k y is the model y-direction permeability, m 2; mu is the viscosity of the injected fluid, pa.s.
④ Hydraulic fracture width calculation equation:
x direction:
y direction:
z direction:
Wherein: u x、uy、uz is the displacement of the matrix rock mass in the x, y and z directions, m/s respectively; c is the elastic constant of the matrix rock mass, and subscript represents the positive and tangential directions of stress without dimension; alpha is the elastic coefficient of the hole, and is generally 0.6-0.8; p s is pore pressure, MPa.
⑤ Hydraulic fracture propagation criteria:
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; k IC is the fracture toughness of the reservoir rock, pa.m 1/2;wtip is the hydraulic fracture tip width, m; Δx is the length of the hydraulic fracture length direction cell, m;
As shown in fig. 1, the self-adaptive time step calculation method for calculating and accelerating the stable three-dimensional hydraulic fracture expansion comprises the following steps:
Step 1, assume that the current calculated time step number t i is that the area occupied by the crack is Ω hf, the crack length is L hf, the crack height is H hf, and the crack width distribution is w=f (x, z);
As shown in fig. 1, the three-dimensional hydraulic fracture propagation model determines that the number of grids in the length L direction is numX, the number of grids in the width W direction is numY, and the number of grids in the height H direction is numZ;
Step 2, calculating the size d cut of the cutting area and the reserved grid number of the matrix rock mass area in each direction, determining a corresponding grid reserved interval, and cutting the matrix rock mass area of the three-dimensional hydraulic fracture expansion model;
The calculation formula of the cutting area size d cut is as follows:
Wherein: d cut is the size of the cutting area, m; a hf is the hydraulic fracture area, m 2;km is the average permeability of the reservoir matrix, m 2,Phf (x, z) is the fluid pressure in the hydraulic fracture area Ω hf, pa; p e is reservoir fluid pressure, pa; mu is the viscosity of the injected fluid and Pa.s; q is injection displacement, m 3/s;tmax is the upper limit of the adaptive time step range, s;
the calculation formula of reserving the grid number in each direction by taking an injection point as an origin in a matrix rock mass area is as follows:
Rx＝Π(dcut/Δx) (9)
Rz＝Π(dcut/Δz) (11)
Wherein: r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δy is the grid length in the y direction, m, where in order to ensure the calculation accuracy, the y direction generally adopts a logarithmic grid; Δz is the grid length in the z direction, m; l hf is the hydraulic fracture length, m; h hf is the hydraulic fracture height, m;
If the injection point (hydraulic fracture center point) is taken as an origin, and the grid position number of the injection point is set to be (0, 0), the grid reservation intervals in different directions after cutting are as follows:
x direction:
[-∏(Lhf/2Δx)-Rx,∏(Lhf/2Δx)+Rx]
y direction:
[-Ry,Ry]
z direction:
[Π(Hhf/2Δz)-Rz,∏(Hhf/2Δz)+Rz]
Wherein: r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; l hf is the hydraulic fracture length, m; h hf is the hydraulic fracture height, m;
Step 3, after cutting of the matrix rock mass area is completed, in order to avoid iteration of the flow field, the matrix seepage field and the fluid loss in the crack, expanding a plane where the hydraulic crack is positioned into a layer of grid, representing the flow capacity of the hydraulic crack by adopting equivalent permeability and porosity of the hydraulic crack, and uniformly calculating three processes of flow in the crack, matrix seepage and fluid loss into three-dimensional equivalent seepage so as to calculate the width of the expanded grid and the equivalent permeability in the hydraulic crack area;
The average width of the water force crack is generally obtained by expanding the grid width, namely, the calculation formula of the expanded grid width is as follows:
wherein: y kc is the expanded grid width, m; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; a hf is the hydraulic fracture area, m 2; w (x, z) is the width distribution of the hydraulic fracture, m;
At the moment, ignoring the width change of the crack flowing in a short time (less than or equal to t max), and integrally calculating the hydraulic crack, the matrix rock and the fluid flow between the two by adopting a three-dimensional seepage equation in the matrix rock. Namely:
Wherein: mu is the viscosity of the injected fluid and Pa.s; k x is the model x-direction permeability, m 2;ky is the model y-direction permeability, m 2;kz is the model z-direction permeability, m 2; phi is porosity, dimensionless; p m is the fluid pressure, pa; t is the calculation time, s; q f is the source term due to injection, s -1;
the equivalent permeability in the hydraulic fracture area is calculated according to cube law and hydraulic fracture width, namely:
Wherein: y kc is the expanded grid width, m; k x,hf (x, z) is the equivalent permeability in the x direction in the hydraulic fracture region Ω hf, and m 2;ky,hf (x, z) is the equivalent permeability in the y direction in the hydraulic fracture region Ω hf, m 2; w (x, z) is the width distribution of the hydraulic fracture, m;
Step 4, setting self-adaptive time step optimizing calculation time according to the average permeability of the reservoir matrix, wherein the lower limit t min of the step length range is generally 0.01s, and the upper limit of the step length range is generally t max; the method comprises the following steps:
Wherein: k m is the reservoir matrix average permeability; t ini is a preselected time step, s;
Step 5, selecting a proper preselected time step according to the established cut area model, and solving the formula (16) to obtain the fluid pressure distribution P hf (x, z) of the hydraulic fracture in four time steps of future injection;
Step 6, finding a grid corresponding to the maximum fluid pressure value in the grid at the tip of the hydraulic fracture, defining the grid pressure as P tip,max, and injecting four time steps into the grid in the future to obtain calculation results of the maximum fluid pressure value at the tip of the fracture, wherein the calculation results are P tip,max(1)、Ptip,max(2)、Ptip,max(3)、Ptip,max (4);
Step 7, in order to ensure the expansion precision of hydraulic fracture calculation and minimize the calculation cost, the upper limit of pressure change in the fracture tip grid in the time step must be smaller than and close to the critical pressure of hydraulic fracture expansion; calculating the hydraulic fracture expansion critical width based on the calculated hydraulic fracture expansion critical width in the formula (7), and calculating the corresponding hydraulic fracture expansion critical pressure P ic according to the hydraulic fracture expansion critical width;
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; w tip is the critical width of hydraulic fracture propagation, m; p ic is hydraulic fracture propagation critical pressure, pa; Δz is the grid length in the z direction, m; sigma oh is the minimum horizontal principal stress, pa;
step 8, judging whether the preselected time step meets the requirements according to the hydraulic fracture expansion critical pressures P ic and P tip,max(3)、Ptip,max (4), if not, calculating a new preselected time step, and repeating the steps 5-8 until the time step meets the requirements, and then carrying out the next step;
If the following formula is satisfied, the current pre-selected time step is judged to be in accordance with the requirement, and t ini is the calculated self-adaptive time step:
Ptip,max(3)≤Pic≤Ptip,max(4) (21)
If not, calculating a new preselected time step based on the relative relationship of P tip,max (3) or P tip,max (4) to P ic;
if P ic≤Ptip,max (3), then take:
If P tip,max(4)≤Pic, then take:
Wherein: p ic is critical pressure, pa; t ini is a preselected time step, s; t' ini is a new preselected time step, s;
Wherein the time step is adopted, and the main program iterative computation can be performed in combination with formulas (1) - (7);
Step 9, determining whether the next time step needs to recalculate the self-adaptive time step according to the hydraulic fracture area increment caused by the hydraulic fracture expansion in the current time step; if the self-adaptive time step needs to be recalculated, repeating the steps 1-8 to obtain a new self-adaptive time step;
The condition of the self-adaptive time step recalculation is that the area change of the hydraulic fracture caused by the expansion of the hydraulic fracture meets the following conditions:
Ahf(ti)/Ahf(tn)≥1.1 (24)
Wherein: a hf is the hydraulic fracture area, and m 2;ti is the current time step number; t n is the number of time steps for obtaining the current adaptive time step;
And step 10, repeating the steps S1-S9, so that the accelerated self-adaptive time step calculation under the condition of guaranteeing the stability of the three-dimensional hydraulic fracture expansion calculation can be realized.
Examples
1. Let the number of meshes in the length L direction be numX =200, the number of meshes in the width W direction be numY =21 (the W direction generally takes the form of logarithmic meshes), and the number of meshes in the height H direction be numZ =100 in a three-dimensional matrix rock mass region of 400×200×50m size. Hydraulic fracture length L hf = 72m, hydraulic fracture height H hf = 64m, hydraulic fracture width w (x, z) as shown in fig. 4. The time step number t i =24 at this time.
2. At this time, the hydraulic fracture area a hf＝4096m2, the average permeability k m＝10×10-15m2 of the reservoir matrix, the maximum fluid pressure P hf＝132×106 Pa in the hydraulic fracture, the reservoir fluid pressure P e＝60×106 Pa, the injection fluid viscosity μ=0.025 pa·s, and the injection displacement q=0.1 m 3/s.
The cutting area size d cut is calculated to be 2.36m according to equation (8).
Then the number of reserved grids in each direction of the matrix rock mass area is calculated according to the formulas (9) - (11): r x＝2,Ry＝7,Rz = 2.
In the original point grid number numX 0＝101,numY0＝11,numZ0 =51, according to equations (12) to (14), the interval of the reserved grids in different directions after cutting can be calculated as the x direction: [81, 121]; y direction: [4, 18]; z direction: [33, 69].
The number of grids after model cutting is 41×15×37=7770, which is reduced by 98% compared with the original model 200×100×21=42000.
3. After the cutting of the matrix rock mass region is completed, the plane in which the hydraulic fracture is located is expanded into a layer of grid, and the expanded grid width y kc is calculated to be 2.1 multiplied by 10 -3 m according to the formula (15).
Based on the hydraulic fracture width distribution w (x, z) and equations (17) to (18), the equivalent permeability in the hydraulic fracture region can be calculated, and the calculation result is shown in fig. 5.
4. The average permeability k m = 10mD of the reservoir matrix, a preselected time step t ini = 1s is selected,
5. The fluid pressure distribution P hf (x, z) of the hydraulic fracture over four time steps of future injection is solved according to equation (16) as shown in fig. 6.
6. Searching a hydraulic fracture tip grid, wherein the corresponding calculation result of the maximum fluid pressure of the fracture tip of four time steps injected in the future is that the grid length delta x=delta y=2m in the x/y direction in a Ptip,max(1)＝1.2221×108Pa,Ptip,max(2)＝1.2563×108Pa,Ptip,max(3)＝1.2796×108Pa,Ptip,max(4)＝1.2967×108Pa. model, the fracture toughness K IC＝7×105Pa·m-1 and the Young modulus E=50000× 6 Pa are respectively calculated, and w tip =0.0011 m is calculated according to a formula (7); further, the hydraulic fracture propagation critical pressure P ic＝1.2832×108 Pa is calculated according to the formula (20). The relationship between P tip,max(3)、Ptip,max (4) and P ic satisfies equation (21), and the present preselected time step t ini =1s.
7. The main routine iterative calculation can be performed by combining equations (1) to (7). Meanwhile, the adaptive time step should be calculated again according to steps 1 to 5 at an appropriate time, and n=1 should be calculated according to equation (24), and the adaptive time step should be calculated again for the next time step.
8. And (3) repeating the steps 1-7, so that the accelerated self-adaptive time step calculation under the condition of guaranteeing the stability of the three-dimensional hydraulic fracture expansion calculation can be realized.
The present invention is not limited to the above-mentioned embodiments, but is not limited to the above-mentioned embodiments, and any person skilled in the art can make some changes or modifications to the equivalent embodiments without departing from the scope of the technical solution of the present invention, but any simple modification, equivalent changes and modifications to the above-mentioned embodiments according to the technical substance of the present invention are still within the scope of the technical solution of the present invention.
Claims (10)
1. The self-adaptive time step calculation method for calculating and accelerating the expansion of the stable three-dimensional hydraulic fracture is characterized by comprising the following steps of:
Step S1, acquiring a hydraulic fracture area omega hf, a hydraulic fracture length L hf, a hydraulic fracture height H hf and a hydraulic fracture width distribution w (x, z) when the current time step number t i is acquired;
s2, calculating the size d cut of a cutting area and the reserved grid number of the matrix rock mass area in each direction, determining a corresponding grid reserved interval, and cutting the matrix rock mass area of the three-dimensional hydraulic fracture expansion model;
s3, after cutting of the matrix rock mass area is completed, expanding a plane where the hydraulic fracture is located into a grid, and calculating the width of the expanded grid and the equivalent permeability in the hydraulic fracture area;
s4, selecting a preselected time step t ini according to the average permeability of the reservoir matrix;
step S5, calculating fluid pressure distribution P hf (x, z) of the hydraulic fracture in four time steps of future injection according to a preselected time step t ini;
Step S6, finding a grid corresponding to the maximum fluid pressure in the grid at the tip of the hydraulic fracture, defining the grid pressure as P tip,max, and injecting four time steps into the grid in the future to obtain calculation results of the maximum fluid pressure at the tip of the fracture, wherein the calculation results are P tip,max(1)、Ptip,max(2)、Ptip,max(3)、Ptip,max (4);
Step S7, calculating a hydraulic fracture expansion critical width, and calculating a corresponding hydraulic fracture expansion critical pressure P ic according to the hydraulic fracture expansion critical width;
Step S8, judging whether the preselected time step meets the requirements according to the hydraulic fracture expansion critical pressures P ic and P tip,max(3)、Ptip,max (4), if not, calculating a new preselected time step, and repeating the steps S5-S8 until the preselected time step meets the requirements, and then carrying out the next step;
Step S9, determining whether the next time step needs to recalculate the self-adaptive time step according to the hydraulic fracture area increment caused by the hydraulic fracture expansion in the current time step; repeating the steps S1-S8 to obtain a new self-adaptive time step if the self-adaptive time step needs to be recalculated;
And step S10, repeating the steps S1-S9, so that the accelerated self-adaptive time step calculation under the condition of guaranteeing the stability of the three-dimensional hydraulic fracture expansion calculation can be realized.
2. The adaptive time-step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the calculation formula of the cutting area size d cut in step S2 is:
Wherein: d cut is the size of the cutting area, m; a hf is the hydraulic fracture area, m 2;km is the average permeability of the reservoir matrix, m 2,Phf (x, z) is the fluid pressure in the hydraulic fracture area Ω hf, pa; p e is reservoir fluid pressure, pa; mu is the viscosity of the injected fluid and Pa.s; q is the injection displacement, m 3/s;tmax is the upper limit of the adaptive time step range, s.
3. The adaptive time-step calculation method for stable three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the calculation formula for reserving the grid number in each direction in the matrix rock mass area in S2 with the injection point as the origin is as follows:
Rx＝Π(dcut/Δx)
Rz＝Π(dcut/Δz)
wherein: d cut is the size of the cutting area, m; r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δy is the grid length in the y direction, m, where the y direction adopts a logarithmic grid to ensure the calculation accuracy; Δz is the grid length in the z direction, m;
If the injection point is taken as an origin, the injection point is taken as a hydraulic fracture center point, and the grid position numbers of the injection point are (0, 0), the grid reserved intervals in different directions after cutting are as follows:
x direction:
[-Π(Lhf/2Δx)-Rx,Π(Lhf/2Δx)+Rx]
y direction:
[-Ry,Ry]
z direction:
[Π(Hhf/2Δz)-Rz,Π(Hhf/2Δz)+Rz]
Wherein: r x is the reserved grid number in the x direction in the matrix rock mass region; r y is the reserved grid number in the y direction in the matrix rock mass area; r z is the reserved grid number in the z direction in the matrix rock mass region; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; l hf is the hydraulic fracture length, m; h hf is the hydraulic fracture height, m.
4. The adaptive time-step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the calculation formula for expanding the mesh width in step S3 is:
wherein: y kc is the expanded grid width, m; Δx is the grid length in the x direction, m; Δz is the grid length in the z direction, m; a hf is the hydraulic fracture area, m 2; w (x, z) is the hydraulic fracture width distribution, m.
5. The adaptive time-step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the calculation formula of the equivalent permeability in the hydraulic fracture area in step S3 is as follows:
Wherein: y kc is the expanded grid width, m; k x,hf (x, z) is the equivalent permeability in the x direction in the hydraulic fracture region Ω hf, and m 2;ky,hf (x, z) is the equivalent permeability in the y direction in the hydraulic fracture region Ω hf, m 2; w (x, z) is the hydraulic fracture width distribution, m.
6. The adaptive time step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the selection of the preselected time step t ini in step S4 is as follows:
Wherein: k m is the reservoir matrix average permeability; t ini is a preselected time step, s.
7. The adaptive time-step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 1, wherein the fluid pressure distribution in step S5 is calculated by the following formula:
Wherein: mu is the viscosity of the injected fluid and Pa.s; k x is the model x-direction permeability, m 2;ky is the model y-direction permeability, m 2;kz is the model z-direction permeability, m 2; phi is porosity, dimensionless; p m is the fluid pressure, pa; t is the calculation time, s; q f is the source term due to injection, s -1.
8. The adaptive time-step calculation method for calculating and accelerating the propagation of the stable three-dimensional hydraulic fracture according to claim 1, wherein the calculation formula of the critical width of the propagation of the hydraulic fracture in step S7 is as follows:
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; k IC is fracture toughness of reservoir rock, pa.m 1/2;wtip is critical width of hydraulic fracture propagation, and m; Δx is the grid length in the x direction, m;
The calculation formula of the hydraulic fracture expansion critical pressure P ic is as follows:
Wherein: e is the elastic modulus of reservoir rock and Pa; v is the reservoir rock poisson ratio, dimensionless; w tip is the critical width of hydraulic fracture propagation, m; p ic is hydraulic fracture propagation critical pressure, pa; Δz is the grid length in the z direction, m; σ oh is the minimum horizontal principal stress, pa.
9. The method for calculating and accelerating the self-adaptive time step for calculating the propagation of the stable three-dimensional hydraulic fracture according to claim 1, wherein the relationship between P ic and P tip,max(3)、Ptip,max (4) in the step S8 is determined to be satisfactory if the relationship satisfies the following equation, and t ini is the calculated self-adaptive time step:
Ptip,max(3)≤Pic≤Ptip,max(4)
If not, a new preselected time step is calculated based on the relative relationship of P tip,max (3) or P tip,max (4) to P ic,
If P ic≤Ptip,max (3), then take:
if P tip,max(4)＜Pic, then take:
Wherein: p ic is hydraulic fracture propagation critical pressure, pa; t ini is a preselected time step, s; t i'ni is a new preselected time step, s.
10. The adaptive time-step calculation method for stabilizing three-dimensional hydraulic fracture propagation calculation and acceleration according to claim 9, wherein the condition for the adaptive time-step recalculation in step S9 is that the hydraulic fracture area change caused by hydraulic fracture propagation meets the following condition:
Ahf(ti)/Ahf(tn)≥1.1
Wherein: a hf is the hydraulic fracture area, and m 2;ti is the current time step number; t n is the number of time steps to obtain the current adaptive time step.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210771868.6A CN115034161B (en) | 2022-06-30 | 2022-06-30 | Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210771868.6A CN115034161B (en) | 2022-06-30 | 2022-06-30 | Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115034161A CN115034161A (en) | 2022-09-09 |
| CN115034161B true CN115034161B (en) | 2024-05-03 |
Family
ID=83129673
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210771868.6A Active CN115034161B (en) | 2022-06-30 | 2022-06-30 | Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115034161B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117494601B (en) * | 2023-11-06 | 2024-05-07 | 西南石油大学 | Fracture-cavity type reservoir acid fracturing effect evaluation method based on embedded discrete fracture |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2505318A1 (en) * | 2004-04-26 | 2005-10-26 | Schlumberger Canada Limited | Method and apparatus and program storage device for front tracking in hydraulic fracturing simulators |
| CA2914348A1 (en) * | 2015-12-10 | 2016-12-14 | Fanhua Zeng | Method of modelling hydrocarbon production from fractured unconventional formations |
| CN108280275A (en) * | 2018-01-09 | 2018-07-13 | 中国石油大学(华东) | A kind of high prediction technique of tight sand hydraulic fracturing seam |
| CN108680730A (en) * | 2018-06-19 | 2018-10-19 | 长安大学 | Simulator and analogy method are endangered in ground fissure place under a kind of seismic loading |
| CN112576240A (en) * | 2020-12-09 | 2021-03-30 | 中国石油大学(华东) | Method for monitoring hydraulic fracturing fracture based on closed wellbore pressure fluctuation |
| CN113836753A (en) * | 2021-11-26 | 2021-12-24 | 西南石油大学 | Temporary blocking steering ball throwing optimization method between cluster perforation gaps in horizontal well section |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8412500B2 (en) * | 2007-01-29 | 2013-04-02 | Schlumberger Technology Corporation | Simulations for hydraulic fracturing treatments and methods of fracturing naturally fractured formation |
| US20220127940A1 (en) * | 2018-03-21 | 2022-04-28 | ResFrac Corporation | Systems and methods for hydraulic fracture and reservoir simulation |
| CN110609974B (en) * | 2019-09-09 | 2022-09-16 | 西南石油大学 | Acid fracturing fracture dynamic fluid loss calculation method considering wormhole expansion |
-
2022
- 2022-06-30 CN CN202210771868.6A patent/CN115034161B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2505318A1 (en) * | 2004-04-26 | 2005-10-26 | Schlumberger Canada Limited | Method and apparatus and program storage device for front tracking in hydraulic fracturing simulators |
| CA2914348A1 (en) * | 2015-12-10 | 2016-12-14 | Fanhua Zeng | Method of modelling hydrocarbon production from fractured unconventional formations |
| CN108280275A (en) * | 2018-01-09 | 2018-07-13 | 中国石油大学(华东) | A kind of high prediction technique of tight sand hydraulic fracturing seam |
| CN108680730A (en) * | 2018-06-19 | 2018-10-19 | 长安大学 | Simulator and analogy method are endangered in ground fissure place under a kind of seismic loading |
| CN112576240A (en) * | 2020-12-09 | 2021-03-30 | 中国石油大学(华东) | Method for monitoring hydraulic fracturing fracture based on closed wellbore pressure fluctuation |
| CN113836753A (en) * | 2021-11-26 | 2021-12-24 | 西南石油大学 | Temporary blocking steering ball throwing optimization method between cluster perforation gaps in horizontal well section |
Non-Patent Citations (3)
| Title |
|---|
| An innovative concept on deep carbonate reservoir stimulation: Three-dimensional acid fracturing technology;Jianchun Guo;Natural Gas Industry;20201030;第7卷(第5期);全文 * |
| 基于离散裂缝网络的网络裂缝酸化井产能预测模型;任冀川;中国优秀硕士学位论文全文数据库工程科技Ⅰ辑;20150215(第2期);全文 * |
| 有限导流垂直裂缝气井产量计算与数值模拟;齐宗耀;中国优秀硕士学位论文全文数据库工程科技I辑;20160115(第1期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115034161A (en) | 2022-09-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11542801B2 (en) | Optimized design method for temporary blocking agent to promote uniform expansion of fractures produced by fracturing in horizontal wells | |
| CN115034161B (en) | Self-adaptive time step calculation method for stable three-dimensional hydraulic fracture expansion calculation and acceleration | |
| CN106372341B (en) | A kind of modification method of reservoir range method water level storage-capacity curve | |
| CN110765695B (en) | Simulation calculation method for obtaining crack propagation path of concrete gravity dam based on high-order finite element method | |
| CN112069654B (en) | Carbonate acidizing numerical simulation method | |
| CN111456709A (en) | Horizontal well multistage fracturing and staged clustering method based on logging curve | |
| CN107145671A (en) | A kind of numerical reservoir simulation method and system | |
| CN115114834B (en) | Fracturing well test simulation method under complex condition | |
| CN113486556B (en) | Improved efficient automatic history fitting method for oil and gas reservoir | |
| CN110516407B (en) | Method for calculating complexity of multiple clusters of fractured fractures in horizontal well section of fractured reservoir | |
| CN114372428B (en) | Multi-cluster fracturing crack extension trans-scale simulation method in horizontal well section of sandstone reservoir | |
| CN109558614B (en) | Simulation method and system for gas flow in shale gas reservoir multi-scale fracture | |
| CN110105085B (en) | Water cooling method for roller compacted concrete dam in construction period | |
| CN109783991B (en) | Landslide sliding process simulation method without known bottom sliding surface | |
| CN103336907A (en) | Method for quickly calculating static storage capacity of reservoir based on DSI technology | |
| CN109356577B (en) | Compact gas reservoir reserve determination method based on gas layer drilling rate | |
| CN104747154A (en) | Method for improving research precision of steam flooding residual oil by utilizing oil displacement efficiency ratio | |
| CN114547998B (en) | Method for determining fracturing modification volume of horizontal well through coupling reservoir fluid | |
| CN109389684B (en) | Numerical simulation method for equivalence of zonal weighting media of fracture-cave oil reservoir | |
| CN112241593B (en) | Fractured reservoir fluid loss calculation method based on multiple time steps | |
| CN110991084A (en) | Reservoir permeability calculation method based on streamline numerical well testing | |
| CN111188613A (en) | Method and system for determining well control radius of tight gas reservoir gas well | |
| CN117421939B (en) | Shale oil fracture system simulation agent method based on track piecewise linearization | |
| CN116861714B (en) | Method for determining water flooding sweep degree of fracture-cavity oil reservoir | |
| CN109145502A (en) | A kind of Weibull type cell life estimation of distribution parameters method |
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 |