CN113821962A - Crack form prediction method considering interlayer weak surface - Google Patents

Crack form prediction method considering interlayer weak surface Download PDF

Info

Publication number
CN113821962A
CN113821962A CN202111185234.4A CN202111185234A CN113821962A CN 113821962 A CN113821962 A CN 113821962A CN 202111185234 A CN202111185234 A CN 202111185234A CN 113821962 A CN113821962 A CN 113821962A
Authority
CN
China
Prior art keywords
fracture
fluid
flow
stress
unit
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.)
Pending
Application number
CN202111185234.4A
Other languages
Chinese (zh)
Inventor
罗志锋
曾秀权
赵立强
沈杰
孙松
李建斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Petroleum University
Original Assignee
Southwest Petroleum University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southwest Petroleum University filed Critical Southwest Petroleum University
Priority to CN202111185234.4A priority Critical patent/CN113821962A/en
Publication of CN113821962A publication Critical patent/CN113821962A/en
Pending legal-status Critical Current

Links

Images

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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • 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/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Data Mining & Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

The invention discloses a crack form prediction method considering an interlayer weak surface, which comprises the following steps: establishing a three-dimensional geological model according to the acquired stratum parameters of the actual interactive stratum; establishing a rock seepage-stress coupling model based on a stress balance equation, a continuity equation and boundary conditions; determining a fracture extension propagation criterion; inputting parameters of a fracturing well of a target well region, and establishing a fracturing three-dimensional fully-coupled model; and simulating the influence of different geological factors and engineering factors on the fracture morphology. The method is based on a finite element method, a full three-dimensional seepage-stress coupling model is established, 2 fracture surfaces are preset by using a Cohesive pore pressure unit, interlayer weak surfaces or bedding weak surfaces are considered, the influence of different geological factors and engineering factors on the fracture morphology is researched, the fracture morphology can be rapidly predicted when a reservoir containing the interlayer weak surfaces is fractured, and the fracture construction parameter design is guided.

Description

Crack form prediction method considering interlayer weak surface
Technical Field
The invention relates to the technical field of oil-gas exploration and development, in particular to a crack form prediction method considering an interlayer weak surface.
Background
Along with the gradual shift of global oil and gas resource exploration and development to unconventional oil and gas development, the hydraulic fracturing technology is widely applied in China at present as one of production increasing measures of unconventional oil and gas wells. And the unconventional oil and gas reservoirs have strong heterogeneity in the longitudinal direction, for example, compact sandstone exists in a sand-shale interbedded layer and a mudstone interlayer, shale exists in a shale-sandstone interbedded layer and a shale-carbonate interbedded layer, and the interbedded causes complex fracture morphology and difficult prediction in hydraulic fracturing. How to predict the fracture state according to the actual stratum interaction condition to guide the optimization of construction parameters is one of the problems faced at present.
In order to research the extension form of hydraulic fracture in interactive stratum, at present, a plurality of scholars mainly build a model from two-dimensional and three-dimensional angles for research.
Two-dimensional model research aspect: patent document No. 202010545763.X discloses a multi-cluster fracture form prediction method based on discrete element fluid-solid coupling, which can effectively consider the influence of interlayer weak surfaces on fracture forms, but does not consider the flow of liquid in a third-dimensional direction, has a certain difference with an actual three-dimensional stratum, and has a certain limitation.
In the aspect of three-dimensional model research, patent document No. 201910470372.3 discloses a cross-layer fracturing chart construction method, which is used for establishing a three-dimensional model and researching the height extension rule of a sand-shale interaction reservoir fracture, but does not consider the influence of interlayer weak surfaces on fracture forms.
The existing method can consider that the research of the fracture form of the interlayer weak face is limited to a two-dimensional model, a certain gap exists between the two-dimensional model and an actual three-dimensional stratum, and the influence of the interlayer weak face on the fracture form is difficult to consider by a three-dimensional model.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides a fracture form prediction method considering an interlayer weak surface, which carries out fracture form prediction by establishing a full three-dimensional seepage-stress coupling model and considering the influence of the interlayer weak surface on hydraulic fracturing.
The technical scheme adopted by the invention for solving the technical problems comprises the following contents:
a crack form prediction method considering interlayer weak planes comprises the following steps:
s1, establishing a three-dimensional geological model according to the acquired stratum parameters of the actual interactive stratum;
s2, establishing a rock seepage-stress coupling model based on a stress balance equation, a continuity equation and boundary conditions;
s3, determining a crack extension and expansion criterion;
s4, inputting fracturing well parameters of a target well area, and establishing a fracturing three-dimensional fully-coupled model;
and S5, simulating the influence of different geological factors and engineering factors on the fracture form.
Preferably, the method for establishing the rock seepage-stress coupling model in step S2 includes:
(1) based on the porous medium theory, according to the virtual work principle, a rock stress balance equation is established:
Figure BDA0003298902900000011
in the formula (I), the compound is shown in the specification,
Figure BDA0003298902900000012
effective stress born by a solid framework, Pa; pnwNon-wetting phase pressure, Pa; i is [1,1,1,0,0,0 ]]TAn identity matrix; δ ε is the imaginary strain rate matrix, s-1(ii) a t is a surface force matrix, N/m3(ii) a δ v is a virtual velocity matrix, m/s; f is a physical strength matrix, N/m3(ii) a dV is a unit body, m3(ii) a S is the integration space surface, m2；
(2) Based on the principle of conservation of materials, the fluid flow in the rock is considered to meet Darcy's law, and a fluid continuity equation is established:
Figure BDA0003298902900000021
in the formula, J is the volume change rate of the porous medium and has no dimension; n isωIs the ratio of the fluid volume to the total volume, dimensionless; rhoωFluid density, kg/m3(ii) a x is a space vector, m; dt is the time increment step, s; v. ofωIs the seepage velocity, m/s;
(3) respectively setting a flow boundary condition, a pore pressure boundary condition and a position boundary condition:
wherein, the flow boundary conditions are as follows:
Figure BDA0003298902900000022
wherein n is the unit normal direction of the flow boundary; k is permeability coefficient tensor; q. q.s0Is the total amount of liquid flowing through the boundary per unit time;
wherein the pore pressure boundary condition can be P0＝P0bIndicating that the pore pressure at the boundary is considered to be a constant value P0b；
Wherein the position boundary condition is used for restricting X, Y, Z direction boundary node displacement to make Ux＝0，Uy＝0，Uz＝0；
(4) By defining a shape function and adopting a Galerkin method to carry out finite element dispersion, a seepage-stress coupling equation can be obtained:
Figure BDA0003298902900000023
in the formula: k ═ loop-V BTDepBdV；
Figure BDA0003298902900000024
Figure BDA0003298902900000025
Figure BDA0003298902900000026
Figure BDA0003298902900000027
Wherein N isuB is a defined shape function vector matrix;
Figure BDA0003298902900000028
is the node displacement;
Figure BDA0003298902900000029
the cell node pore pressure is obtained; n is a radical ofpIs a shape function; depIs an elastic matrix; m is a conversion matrix; q. q.sobIs the fluid flow on the boundary.
Preferably, the rock seepage-stress coupling model comprises a dynamic rock physical parameter model, and the dynamic rock physical parameter model is as follows:
(1) multi-field coupled dynamic porosity model:
Figure BDA00032989029000000210
in the formula, omega is the thermal expansion coefficient of 1/k; epsilonvIs volume strain, dimensionless;
(2) dynamic model of permeability based on Kozeyn-canann capillary bundle model derivation:
Figure BDA00032989029000000211
in the formula, epsilonvIs volume strain, dimensionless; omega is thermal expansion coefficient, 1/k; Δ T is the temperature change value, K;
(3) dynamic model of pore compressibility:
Figure BDA0003298902900000031
in the formula, epsilonvIs volume strain, dimensionless; Δ T is the temperature change value, K; Δ p is the pressure change value, MPa.
Preferably, in step S3, the fracture propagation criterion includes using a viscoelastic Cohesive pore pressure unit to preset 2 fracture surfaces to simulate fracture initiation and propagation, where the fracture surface I is perpendicular to the direction of minimum principal stress, the fracture surface II is present at a formation layer boundary, and the fracture surface I and the fracture II intersect.
Preferably, 4 pore pressure points exist at the intersection unit of the fracture surface I and the fracture surface II, and the continuity condition of the fluid pressure at the intersection point is met by replacing the 4 pore pressure nodes with a common pore pressure node.
Preferably, the fracture initiation conditions are as follows:
Figure BDA0003298902900000032
in the formula, σnThe applied stress in the normal direction of the Colesive unit is Mpa; tau iss、τtThe applied stress is two tangential directions of the unit, MPa;
Figure BDA0003298902900000033
the critical stress in the normal direction is MPa when the unit fails;
Figure BDA0003298902900000034
the critical stress in two directions, Mpa, is when the unit fails tangentially.
Preferably, the Cohesive unit enters a damage evolution stage after being cracked and damaged, namely the Cohesive unit is damaged continuously so that the unit is completely damaged when the fracture energy of the fracture is reached, the fracture is finally formed, and a damage factor calculation formula is as follows:
Figure BDA0003298902900000035
in the formula (I), the compound is shown in the specification,
Figure BDA0003298902900000036
m is the maximum displacement reached by the unit during loading;
Figure BDA0003298902900000037
displacement when the cell is completely destroyed, m;
Figure BDA0003298902900000038
is the displacement of the cell at initial damage, m;
wherein, the unit damage condition is judged by the damage factor D: when no damage has formed, the value of D is 0; when the cell was completely damaged, the D value was 1.
Preferably, the injected liquid is Newtonian fluid during fracturing, and the flow model of the simulated fracturing fluid in the fracture surface comprises tangential flow along the fracture surface and normal flow perpendicular to the fracture surface:
(1) the tangential flow along the fracture plane, defined by newtonian flow pressure conduction, is:
Q＝-ktΔp
wherein Q is the discharge capacity of the fracturing fluid, m3/s；ktIs the flow coefficient; p is the flow pressure, MPa;
wherein the flow coefficient ktCan be expressed as:
Figure BDA0003298902900000039
wherein d is the crack opening width; mu is fracturing fluid viscosity, Pa.s;
(2) the normal flow perpendicular to the fracture surface simulates seepage fluid loss in the normal direction of the fracture surface and is characterized by defining a fluid loss coefficient, and the normal flow is expressed as:
Figure BDA00032989029000000310
in the formula, qt、qbRespectively the seepage flow of the fluid on the upper and lower surfaces of the fracture, m3/s；ct、cbThe fluid loss coefficients of the fluid on the upper surface and the lower surface of the fracture are respectively; p is a radical oft、pbRespectively the pore pressure of the body on the upper and lower surfaces of the crack, Mpa; p is a radical ofiThe fluid pressure of the middle surface of the Colesive fracture unit is Mpa.
Preferably, in the step S5, the geological factor includes any one or more of reservoir thickness, interlayer thickness, number of layers, tensile strength of layer interface, petrophysical parameters, rock mechanical parameters and stress difference; the engineering factors include any one or more of construction displacement, perforation orientation, fracturing fluid viscosity, proppant concentration, and construction scale.
The invention has the beneficial effects that:
the method is based on a finite element method, a full three-dimensional seepage-stress coupling model is established, 2 crack surfaces are preset by using a Cohesive pore pressure unit, interlayer weak surfaces or bedding weak surfaces are considered, the influence of different geological factors and engineering factors on the crack form is researched, the design of fracturing construction parameters can be guided, and the problem that the influence of the interlayer weak surfaces or the bedding weak surfaces on crack expansion is difficult to research in a three-dimensional angle at present is solved preliminarily.
Drawings
FIG. 1 is a schematic flow chart of a crack shape prediction method considering an interlayer weak plane according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of a three-dimensional geological model of a shale-sandstone interaction formation, according to an embodiment of the invention;
FIG. 3 is a schematic illustration of a predetermined intersecting fracture plane I and fracture plane II according to one embodiment of the present invention;
FIG. 4 is a schematic view of the processing of merging the intersecting nodes of fracture surface I and fracture surface II;
FIG. 5 is a schematic diagram of a fracture three-dimensional fully-coupled model constructed according to an embodiment of the invention;
FIG. 6 is a schematic view of a calculated fracture morphology according to one embodiment of the present invention;
FIG. 7 is a graphical representation of calculated seam height versus layer interface tensile strength in accordance with one embodiment of the present invention;
FIG. 8 is a schematic illustration of a calculated relationship between seam height and construction displacement according to an embodiment of the present invention;
FIG. 9 is a graphical representation of calculated fracture height versus fracturing fluid viscosity according to one embodiment of the present invention.
Detailed Description
The following detailed description of the embodiments of the present invention will be provided with reference to the drawings and examples, so that how to apply the technical means to solve the technical problems and achieve the technical effects can be fully understood and implemented. It should be noted that, as long as there is no conflict, the embodiments and the features of the embodiments of the present invention may be combined with each other, and the technical solutions formed are within the scope of the present invention.
In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the embodiments of the invention. It will be apparent, however, to one skilled in the art that the present invention may be practiced without some of these specific details or with other methods described herein.
Additionally, the steps illustrated in the flow charts of the figures may be performed in a computer system such as a set of computer-executable instructions and, although a logical order is illustrated in the flow charts, in some cases, the steps illustrated or described may be performed in an order different than here.
The invention is further illustrated with reference to the following figures and examples.
A crack form prediction method considering interlayer weak planes comprises the following steps:
s1, establishing a three-dimensional geological model according to the acquired stratum parameters of the actual interactive stratum;
s2, establishing a rock seepage-stress coupling model based on a stress balance equation, a continuity equation and boundary conditions;
s3, determining a crack extension and expansion criterion;
s4, inputting fracturing well parameters of a target well area, and establishing a fracturing three-dimensional fully-coupled model;
and S5, simulating the influence of different geological factors and engineering factors on the fracture form.
Fig. 1 shows a schematic implementation flow diagram of the crack morphology prediction method considering an interlayer weak plane according to this embodiment.
In this embodiment, preferably, an Abaqus finite element simulation calculation tool is used as a platform to establish a three-dimensional geological model, establish a rock seepage-stress coupling model, preset 2 fracture surfaces by using a Cohesive pore pressure unit, and simulate fracture initiation and propagation according to a corresponding fracture extension rule. Of course, in other embodiments of the present invention, according to actual needs, the method may also use other reasonable tools or manners to perform fracture morphology prediction on the reservoir containing the interbed weak face.
In this embodiment, the method preferably first establishes a three-dimensional geological model of the shale-sandstone interaction formation in step S1. The shale-sandstone interaction stratum has the geological characteristics that: the three-dimensional geological model is built by arranging the perforation sections on the middle shale layer, the shale layers are shale (4m), sandstone (4m), shale (10m), sandstone (4m), shale (8m), sandstone (4m), shale (4m), sandstone (4m), shale (2m), sandstone (4m) and shale (2m) from top to bottom, and the three-dimensional geological model is shown in figure 2.
The method for establishing the rock seepage-stress coupling model in the step S2 comprises the following steps:
(1) based on the porous medium theory, according to the virtual work principle, a rock stress balance equation is established:
Figure BDA0003298902900000051
in the formula (I), the compound is shown in the specification,
Figure BDA0003298902900000052
effective stress born by a solid framework, Pa; pnwNon-wetting phase pressure, Pa; i is [1,1,1,0,0,0 ]]TAn identity matrix; δ ε is the imaginary strain rate matrix, s-1(ii) a t is a surface force matrix, N/m3(ii) a δ v is a virtual velocity matrix, m/s; f is a physical strength matrix, N/m3(ii) a dV is a unit body, m3(ii) a S is the integration space surface, m2；
(2) Based on the principle of conservation of materials, the fluid flow in the rock is considered to meet Darcy's law, and a fluid continuity equation is established:
Figure BDA0003298902900000053
in the formula, J is the volume change rate of the porous medium and has no dimension; n isωIs the ratio of the fluid volume to the total volume, dimensionless; rhoωFluid density, kg/m3(ii) a x is a space vector, m; dt is the time increment step, s; v. ofωIs the seepage velocity, m/s;
(3) respectively setting a flow boundary condition, a pore pressure boundary condition and a position boundary condition:
wherein, the flow boundary conditions are as follows:
Figure BDA0003298902900000054
wherein n is the unit normal direction of the flow boundary; k is permeability coefficient tensor; q. q.s0Is the total amount of liquid flowing through the boundary per unit time;
wherein the pore pressure boundary condition can be P0＝P0bIndicating that the pore pressure at the boundary is considered to be a constant value P0b；
Wherein the position boundary condition is used for restricting X, Y, Z direction boundary node displacement to make Ux＝0，Uy＝0，Uz＝0；
(4) By defining a shape function and adopting a Galerkin method to carry out finite element dispersion, a seepage-stress coupling equation can be obtained:
Figure BDA0003298902900000055
in the formula: k ═ loop-V BTDepBdV；
Figure BDA0003298902900000056
Figure BDA0003298902900000057
Figure BDA0003298902900000061
Figure BDA0003298902900000062
Wherein N isuB is a defined shape function vector matrix;
Figure BDA0003298902900000063
is the node displacement;
Figure BDA0003298902900000064
the cell node pore pressure is obtained; n is a radical ofpIs a shape function; depIs an elastic matrix; m is a conversion matrix; q. q.sobIs the fluid flow on the boundary.
Along with the injection of fracturing fluid, the change of fluid pressure can cause the change of petrophysical parameters, and a petrophysical parameter dynamic model needs to be established for updating the petrophysical parameter data in the rock seepage-stress coupling model, wherein the petrophysical parameter dynamic model is as follows:
(1) multi-field coupled dynamic porosity model:
Figure BDA0003298902900000065
in the formula, omega is the thermal expansion coefficient of 1/k; epsilonvIs volume strain, dimensionless;
(2) dynamic model of permeability based on Kozeyn-canann capillary bundle model derivation:
Figure BDA0003298902900000066
in the formula, epsilonvIs volume strain, dimensionless; omega is thermal expansion coefficient, 1/k; Δ T is the temperature change value, K;
(3) dynamic model of pore compressibility:
Figure BDA0003298902900000067
in the formula, epsilonvIs volume strain, dimensionless; Δ T is the temperature change value, K; Δ p is the pressure change value, MPa.
In step S3, the fracture propagation criterion preferably uses a viscoelastic Cohesive pore pressure unit to preset 2 fracture surfaces to simulate fracture initiation and propagation, the fracture surface I is perpendicular to the direction of the minimum principal stress, the fracture surface II exists at the formation layer interface, and the fracture surface I and the fracture II intersect, as shown in fig. 3.
As shown in FIG. 4, 4 pore pressure points exist at the unit where the fracture surface I and the fracture surface II intersect, and P in the intersected unit is obtained by replacing the 4 pore pressure nodes with a common pore pressure node1、P2、P3、P4Pore pressure node is replaced by common pore pressure node P0And the continuity condition of the fluid pressure at the intersection point is met.
A Cohesive pore pressure unit is prefabricated to simulate cracks, and the crack initiation conditions are as follows:
Figure BDA0003298902900000068
in the formula, σnThe applied stress in the normal direction of the Colesive unit is Mpa; tau iss、τtThe applied stress is two tangential directions of the unit, MPa;
Figure BDA0003298902900000069
the critical stress in the normal direction is MPa when the unit fails;
Figure BDA00032989029000000610
the critical stress in two directions, Mpa, is when the unit fails tangentially.
When the Cohesive unit is cracked and damaged, the Cohesive unit enters a damage evolution stage, namely the Cohesive unit is damaged continuously so that the Cohesive unit is completely damaged when the fracture energy of the fracture is reached, and finally the fracture is formed, wherein a damage factor calculation formula is as follows:
Figure BDA0003298902900000071
in the formula (I), the compound is shown in the specification,
Figure BDA0003298902900000072
m is the maximum displacement reached by the unit during loading;
Figure BDA0003298902900000073
displacement when the cell is completely destroyed, m;
Figure BDA0003298902900000074
is the displacement of the cell at initial damage, m;
wherein, the unit damage condition is judged by the damage factor D: when no damage has formed, the value of D is 0; when the cell was completely damaged, the D value was 1.
In this example, the fracturing fluid is regarded as a newtonian fluid, and the model for simulating the flow of the fracturing fluid in the fracture surface comprises a tangential flow along the fracture surface and a normal flow perpendicular to the fracture surface:
(1) the tangential flow along the fracture plane, defined by newtonian flow pressure conduction, is:
Q＝-ktΔp
wherein Q is the discharge capacity of the fracturing fluid, m3/s；ktIs the flow coefficient; p is the flow pressure, MPa;
wherein the flow coefficient ktCan be expressed as:
Figure BDA0003298902900000075
wherein d is the crack opening width; mu is fracturing fluid viscosity, Pa.s;
(2) the normal flow perpendicular to the fracture surface simulates seepage fluid loss in the normal direction of the fracture surface and is characterized by defining a fluid loss coefficient, and the normal flow is expressed as:
Figure BDA0003298902900000076
in the formula, qt、qbRespectively the seepage flow of the fluid on the upper and lower surfaces of the fracture, m3/s；ct、cbThe fluid loss coefficients of the fluid on the upper surface and the lower surface of the fracture are respectively; p is a radical oft、pbRespectively the pore pressure of the body on the upper and lower surfaces of the crack, Mpa; p is a radical ofiThe fluid pressure of the middle surface of the Colesive fracture unit is Mpa.
And step S4, inputting fracturing well parameters of the target well area, and establishing a fracturing three-dimensional fully-coupled model.
As shown in fig. 5, statistical parameters of the fractured well of the target well zone, including the elastic modulus of rock, poisson ratio, tensile strength of rock body, shear strength of rock body, tensile strength of prefabricated fracture surface, ground stress, injection displacement, viscosity of fracturing fluid, initial permeability and initial porosity, are input, and a fractured three-dimensional fully-coupled model is established.
TABLE 1 parameters table for fracturing well of target well
Figure BDA0003298902900000077
Figure BDA0003298902900000081
And (3) solving the fractured three-dimensional fully-coupled model by using an Abaqus implicit solver, wherein as shown in FIG. 6, the weak surface between layers cracks, so that an I-shaped seam is obviously formed, the length of the half seam of the fracture is predicted to be 44m, and the height of the seam is predicted to be 10 m.
And step S5, simulating the influence of different geological factors and engineering factors on the fracture morphology.
(1) Geological factors are changed, the tensile strength of the shale-sandstone layer interface is changed, the crack propagation conditions under the tensile strengths of 0.5, 1.0, 1.5, 2.0 and 3.0MPa are simulated respectively, and the results are shown in FIG. 7:
simulation results show that the lower the tensile strength of the layer interface is, the more easily the layer interface is broken, the more easily the crack extends along the direction of the layer interface, and when the tensile strength of the layer interface exceeds 2MPa, the layer interface is broken in a small range, and the crack is high and easily penetrates through the layer interface.
(2) The engineering factors are changed, the construction discharge capacity is changed, and the construction discharge capacity is simulated to be 8, 10, 12, 14 and 16m respectively3Fracture propagation at/min, results are shown in FIG. 8:
simulation results show that as the construction displacement is increased, cracks are continuously communicated with weak surfaces between layers, and the height of the cracks is increased; the discharge capacity is increased by 1 time, the seam height is increased by 25 percent, and the discharge capacity has certain influence on the seam height.
(3) The engineering factors are changed, the viscosity of the fracturing fluid is changed, and the fracture propagation conditions of the fracturing fluid with the viscosity of 1, 10, 20, 50, 100 and 150mPa.s are simulated respectively, and the results are shown in FIG. 9.
Simulation results show that the viscosity is increased, the crack height is increased, the viscosity is too high, the crack height breaks through the layer interface, and meanwhile, the layer interface is cracked; the fracture fluid viscosity has a large influence on the fracture morphology.
Reference in the specification to "one embodiment" or "an embodiment" means that a particular feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment of the invention. Thus, the appearances of the phrase "one embodiment" or "an embodiment" in various places throughout this specification are not necessarily all referring to the same embodiment.
Although the present invention has been described with reference to a preferred embodiment, it should be understood that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (9)

1. A crack form prediction method considering interlayer weak planes is characterized by comprising the following steps:
s1, establishing a three-dimensional geological model according to the acquired stratum parameters of the actual interactive stratum;
s2, establishing a rock seepage-stress coupling model based on a stress balance equation, a continuity equation and boundary conditions;
s3, determining a crack extension and expansion criterion;
s4, inputting fracturing well parameters of a target well area, and establishing a fracturing three-dimensional fully-coupled model;
and S5, simulating the influence of different geological factors and engineering factors on the fracture form.
2. The method for predicting the fracture morphology by considering the interlaminar weak plane as claimed in claim 1, wherein the method for establishing the rock seepage-stress coupling model in the step S2 is as follows:
(1) based on the porous medium theory, according to the virtual work principle, a rock stress balance equation is established:
Figure FDA0003298902890000011
in the formula (I), the compound is shown in the specification,
Figure FDA0003298902890000012
effective stress born by a solid framework, Pa; pnwNon-wetting phase pressure, Pa; i is [1,1,1,0,0,0 ]]TAn identity matrix; δ ε is the imaginary strain rate matrix, s-1(ii) a t is a surface force matrix, N/m3(ii) a δ v is a virtual velocity matrix, m/s; f is a physical strength matrix, N/m3(ii) a dV is a unit body, m3(ii) a S is the integration space surface, m2；
(2) Based on the principle of conservation of materials, the fluid flow in the rock is considered to meet Darcy's law, and a fluid continuity equation is established:
Figure FDA0003298902890000013
in the formula, J is the volume change rate of the porous medium and has no dimension; n isωIs the ratio of the fluid volume to the total volume, dimensionless; rhoωFluid density, kg/m3(ii) a x is a space vector, m; dt is the time increment step, s; v. ofωIs the seepage velocity, m/s;
(3) respectively setting a flow boundary condition, a pore pressure boundary condition and a position boundary condition:
wherein, the flow boundary conditions are as follows:
Figure FDA0003298902890000014
wherein n is the unit normal direction of the flow boundary; k is permeability coefficient tensor; q. q.s0Is the total amount of liquid flowing through the boundary per unit time;
wherein the pore pressure boundary condition can be P0＝P0bIndicating that the pore pressure at the boundary is considered to be a constant value P0b；
Wherein the position boundary condition is used for restricting X, Y, Z direction boundary node displacement to make Ux＝0，Uy＝0，Uz＝0；
(4) By defining a shape function and adopting a Galerkin method to carry out finite element dispersion, a seepage-stress coupling equation can be obtained:
Figure FDA0003298902890000015
in the formula: k ═ loop-VBTDepBdV；
Figure FDA0003298902890000016
Figure FDA0003298902890000017
Figure FDA0003298902890000018
Figure FDA0003298902890000021
Wherein N isuB is a defined shape function vector matrix;
Figure FDA0003298902890000022
is the node displacement;
Figure FDA0003298902890000023
the cell node pore pressure is obtained; n is a radical ofpIs a shape function; depIs an elastic matrix; m is a conversion matrix; q. q.sobIs the fluid flow on the boundary.
3. The method of claim 1 or 2, wherein the rock seepage-stress coupling model comprises a petrophysical parameter dynamic model of:
(1) multi-field coupled dynamic porosity model:
Figure FDA0003298902890000024
in the formula, omega is the thermal expansion coefficient of 1/k; epsilonvIs volume strain, dimensionless;
(2) dynamic model of permeability based on Kozeyn-canann capillary bundle model derivation:
Figure FDA0003298902890000025
in the formula, epsilonvIs volume strain, dimensionless; omega is thermal expansion coefficient, 1/k; Δ T is the temperature change value, K;
(3) dynamic model of pore compressibility:
Figure FDA0003298902890000026
in the formula, epsilonvIs volume strain, dimensionless; Δ T is the temperature change value, K; Δ p is the pressure change value, MPa.
4. The method of claim 1, wherein the fracture propagation criterion in step S3 includes using a viscoelastic Cohesive pore pressure cell to preset 2 fracture planes to simulate fracture initiation and propagation, the fracture plane I is perpendicular to the direction of least principal stress, the fracture plane II is present at a formation level boundary, and the fracture plane I and the fracture plane II intersect.
5. The method of claim 4, wherein 4 pore pressure points exist at the unit where the fracture surface I and the fracture surface II intersect, and the continuity condition of the fluid pressure at the intersection point is satisfied by replacing the 4 pore pressure nodes with one common pore pressure node.
6. The method according to any one of claims 1 to 5, wherein the fracture initiation conditions are:
Figure FDA0003298902890000027
in the formula, σnThe applied stress in the normal direction of the Colesive unit is Mpa; tau iss、τtThe applied stress is two tangential directions of the unit, MPa;
Figure FDA0003298902890000028
the critical stress in the normal direction is MPa when the unit fails;
Figure FDA0003298902890000029
the critical stress in two directions, Mpa, is when the unit fails tangentially.
7. The method according to any one of claims 1 to 6, wherein the Cohesive unit enters a damage evolution stage after initiating damage, namely the Cohesive unit is continuously damaged so that the unit is completely damaged when reaching fracture energy to finally form a fracture, and a damage factor calculation formula is as follows:
Figure FDA00032989028900000210
in the formula (I), the compound is shown in the specification,
Figure FDA00032989028900000211
m is the maximum displacement reached by the unit during loading;
Figure FDA00032989028900000212
displacement when the cell is completely destroyed, m;
Figure FDA00032989028900000213
is the displacement of the cell at initial damage, m;
wherein, the unit damage condition is judged by the damage factor D: when no damage has formed, the value of D is 0; when the cell was completely damaged, the D value was 1.
8. The method as claimed in any one of claims 1 to 7, wherein the injected liquid is Newtonian fluid during fracturing, and the flow model of the simulated fracturing fluid in the fracture surface comprises tangential flow along the fracture surface and normal flow perpendicular to the fracture surface:
(1) the tangential flow along the fracture plane, defined by newtonian flow pressure conduction, is:
Q＝-ktΔp
wherein Q is the discharge capacity of the fracturing fluid, m3/s；ktIs the flow coefficient; p is the flow pressure, MPa;
wherein the flow coefficient ktCan be expressed as:
Figure FDA0003298902890000031
wherein d is the crack opening width; mu is fracturing fluid viscosity, Pa.s;
(2) the normal flow perpendicular to the fracture surface simulates seepage fluid loss in the normal direction of the fracture surface and is characterized by defining a fluid loss coefficient, and the normal flow is expressed as:
Figure FDA0003298902890000032
in the formula, qt、qbRespectively the seepage flow of the fluid on the upper and lower surfaces of the fracture, m3/s；ct、cbThe fluid loss coefficients of the fluid on the upper surface and the lower surface of the fracture are respectively; p is a radical oft、pbRespectively the pore pressure of the body on the upper and lower surfaces of the crack, Mpa; p is a radical ofiThe fluid pressure of the middle surface of the Colesive fracture unit is Mpa.
9. The method for predicting the fracture morphology according to claim 1, wherein in the step S5, the geological factors comprise any one or more of reservoir thickness, interlayer thickness, number of layers, tensile strength of layer interface, petrophysical parameters, rock mechanical parameters and stress difference; the engineering factors include any one or more of construction displacement, perforation orientation, fracturing fluid viscosity, proppant concentration, and construction scale.
CN202111185234.4A 2021-10-12 2021-10-12 Crack form prediction method considering interlayer weak surface Pending CN113821962A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111185234.4A CN113821962A (en) 2021-10-12 2021-10-12 Crack form prediction method considering interlayer weak surface

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111185234.4A CN113821962A (en) 2021-10-12 2021-10-12 Crack form prediction method considering interlayer weak surface

Publications (1)

Publication Number Publication Date
CN113821962A true CN113821962A (en) 2021-12-21

Family

ID=78916449

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111185234.4A Pending CN113821962A (en) 2021-10-12 2021-10-12 Crack form prediction method considering interlayer weak surface

Country Status (1)

Country Link
CN (1) CN113821962A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116956769A (en) * 2023-06-30 2023-10-27 中国地质大学（北京） Numerical simulation method for hydraulic fracturing crack propagation in random fracture stratum

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110134984A (en) * 2019-03-06 2019-08-16 中国石油化工股份有限公司江汉油田分公司石油工程技术研究院 The analysis method of complex fracture extension influence factor in a kind of shale fracturing process
CN110348031A (en) * 2018-04-08 2019-10-18 中国石油化工股份有限公司 The nearly pit shaft crack distorted configurations method for numerical simulation of fractured horizontal well
US20200387650A1 (en) * 2019-06-10 2020-12-10 Southwest Petroleum University Fracturing fluid flow-back simulation method for fractured horizontal well in shale gas reservoir
US20210003731A1 (en) * 2019-07-04 2021-01-07 Chengdu University Of Technology Method for determining favorable time window of infill well in unconventional oil and gas reservoir

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110348031A (en) * 2018-04-08 2019-10-18 中国石油化工股份有限公司 The nearly pit shaft crack distorted configurations method for numerical simulation of fractured horizontal well
CN110134984A (en) * 2019-03-06 2019-08-16 中国石油化工股份有限公司江汉油田分公司石油工程技术研究院 The analysis method of complex fracture extension influence factor in a kind of shale fracturing process
US20200387650A1 (en) * 2019-06-10 2020-12-10 Southwest Petroleum University Fracturing fluid flow-back simulation method for fractured horizontal well in shale gas reservoir
US20210003731A1 (en) * 2019-07-04 2021-01-07 Chengdu University Of Technology Method for determining favorable time window of infill well in unconventional oil and gas reservoir

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
卢忠华;杨传礼;权秀鹃;: "节理发育地层水力裂缝扩展数值研究", 应用力学学报, no. 03, 15 June 2020 (2020-06-15) *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116956769A (en) * 2023-06-30 2023-10-27 中国地质大学（北京） Numerical simulation method for hydraulic fracturing crack propagation in random fracture stratum

Similar Documents

Publication Publication Date Title
Ma et al. A numerical gas fracturing model of coupled thermal, flowing and mechanical effects
Zou et al. 3-D numerical simulation of hydraulic fracturing in a CBM reservoir
CN108280275B (en) Compact sandstone hydraulic fracture height prediction method
CN108468538B (en) Shale hydraulic fracture propagation prediction method
CN108412477B (en) Method for making seam in intermittent partial-sealing and plugging seam in volume fracturing
US10161235B2 (en) Hydraulic fracturing in highly heterogeneous formations by resisting formation and/or sealing micro-fractures
Wang et al. Numerical simulation of fracture initiation, propagation and fracture complexity in the presence of multiple perforations
Cong et al. Numerical simulation of hydraulic fracture height layer-through propagation based on three-dimensional lattice method
CN110348031B (en) Numerical simulation method for horizontal well fracturing near-wellbore fracture distortion form
CN112541287A (en) Loose sandstone fracturing filling sand control production increase and profile control integrated design method
CN111259595B (en) Coal-sand interbedded through-layer fracturing perforation position optimization method
CN114722682A (en) Shale reservoir horizontal well temporary plugging fracturing multi-fracture competition fracture initiation prediction method
CA3062854A1 (en) Method for predicting of hydraulic fracturing and associated risks
CN114580100B (en) Method and device for calculating full wellbore pressure of fractured horizontal well and computer readable storage medium
Li et al. Effects of fluid and proppant properties on proppant transport and distribution in horizontal hydraulic fractures of coal under true-triaxial stresses
CN113821962A (en) Crack form prediction method considering interlayer weak surface
Hou et al. Vertical fracture propagation behavior upon supercritical carbon dioxide fracturing of multiple layers
CN112100707A (en) Construction method of through-layer fracturing plate
CN113987965B (en) Prediction method and device for temporary plugging steering crack
CN112012710A (en) Horizontal well fracturing three-dimensional fracture propagation simulation method for sand-shale interactive stratum
Haddad et al. Cohesive Fracture Analysis to Model Multiple-Stage Fracturing in Quasibrittle Shale Formations
CN114580315A (en) Hydraulic fracturing fracture extension and multiphase fluid flow simulation method
Zhu et al. Numerical simulation of complex fractures propagation in thin sand-shale interbeded formation: A case study of shale oil reservoir in Shengli Oilfield
CN115828440A (en) Construction method of hydraulic fracture extension model of shale oil horizontal well
Yuan et al. Numerical simulation of fracture morphology in Jurassic reservoir

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