CN111259595B - Coal-sand interbedded through-layer fracturing perforation position optimization method - Google Patents

Coal-sand interbedded through-layer fracturing perforation position optimization method Download PDF

Info

Publication number
CN111259595B
CN111259595B CN202010100167.0A CN202010100167A CN111259595B CN 111259595 B CN111259595 B CN 111259595B CN 202010100167 A CN202010100167 A CN 202010100167A CN 111259595 B CN111259595 B CN 111259595B
Authority
CN
China
Prior art keywords
equation
fracture
formula
phase field
mpa
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
Application number
CN202010100167.0A
Other languages
Chinese (zh)
Other versions
CN111259595A (en
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 CN202010100167.0A priority Critical patent/CN111259595B/en
Publication of CN111259595A publication Critical patent/CN111259595A/en
Application granted granted Critical
Publication of CN111259595B publication Critical patent/CN111259595B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • E—FIXED CONSTRUCTIONS
    • E21—EARTH DRILLING; MINING
    • E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • E21B43/006—Production of coal-bed methane
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH DRILLING; MINING
    • E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • E21B43/25—Methods for stimulating production
    • E21B43/26—Methods for stimulating production by forming crevices or fractures
    • 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Abstract

The invention provides a coal-sand interbedded through-layer fracturing perforation position optimization method, which comprises the following steps of: (1) collecting input parameters; (2) establishing a stress calculation equation set; (3) establishing a fluid flow control equation set; (4) establishing a fracture phase field evolution equation set; (5) establishing a longitudinal extension numerical calculation model of the hydraulic fracture in the porous elastic medium by combining the steps (2) to (4); (6) and (3) inputting the parameters obtained in the step (1) into the model established in the step (5), and comparing the different formation parameters with the fracture trajectories under the conditions of the perforation positions, thereby optimizing the perforation positions. The fracture extension track and conditions are self-adaptive, the defect that the track prediction criterion needs to be additionally established to judge the fracture extension direction in the prior art is overcome, and the fluid loss coefficient does not need to be introduced to describe the fluid loss phenomenon of the fracturing fluid.

Description

Coal-sand interbedded through-layer fracturing perforation position optimization method
Technical Field
The invention relates to the field of yield increase transformation of oil and gas fields, in particular to a method for optimizing a coal-sand interbedded fracturing perforation position.
Background
The Chinese coal bed gas field contains rich natural gas resources, the gas field contains rich compact sandstone gas in the longitudinal direction besides coal bed gas, a coal-sand-rock interbed is formed, and the coal rock and sandstone reservoirs are communicated in the longitudinal direction through a hydraulic fracturing technology to form a mode for efficiently developing the gas field. However, due to the influence of the interlayer interface and the lithology and ground stress difference of different layers, if the perforation position is not proper, the hydraulic fracture can only extend in one layer, and therefore, the hydraulic fracture has important significance in communicating a plurality of reservoirs through the manual optimization of the perforation position.
The core of the technology is to establish a longitudinal extension model of the hydraulic fracture in multiple layers. A commonly used numerical simulation method for simulating the longitudinal extension of a crack in multiple layers comprises (1) a Palmer model and an improved model thereof, (2) a Displacement Discontinuity Method (DDM); (3) a finite element method based on the ABAQUS platform; (4) the extended finite element method. However, when the numerical simulation method is used for researching the extension track of the hydraulic fracture after encountering the interlayer interface, an intersection criterion needs to be established to judge whether the hydraulic fracture passes through the interface or opens the interface. In addition, most of the existing fracture extension models need to introduce a Carter fluid loss coefficient to describe the fluid loss phenomenon of the fracturing fluid.
Disclosure of Invention
Aiming at the technical problems, the invention establishes a numerical calculation model of the longitudinal extension of the hydraulic fracture in the porous elastic medium by comprehensively applying multidisciplinary knowledge such as a Biot porous elastic theory, a finite element theory, a phase-field method theory, a nonlinear equation numerical solving method and the like, and forms a coal-sand interbedded fracturing perforation position optimization method based on the model.
A coal-sand interbedded fracturing perforation position optimization method comprises the following steps:
(1) collecting input parameters;
(2) establishing a stress calculation equation set;
(3) establishing a fluid flow control equation set;
(4) establishing a fracture phase field evolution equation set;
(5) establishing a longitudinal extension numerical calculation model of the hydraulic fracture in the porous elastic medium by combining the steps (2) to (4);
(6) and (3) inputting the parameters obtained in the step (1) into the model established in the step (5), and comparing the different formation parameters with the fracture trajectories under the conditions of the perforation positions, thereby optimizing the perforation positions.
Further, the input parameters collected in step (1) include: the method comprises the following steps of geostress parameters, elastic modulus and Poisson ratio of different layers of rock, critical tensile stress of different layers of rock, permeability of different layers of rock, critical tensile stress of interlayer interface, permeability of interlayer interface, fracturing discharge capacity, fracturing fluid injection time and fracturing fluid viscosity.
The step (2) establishes a stress calculation equation set, which comprises the following contents:
(2.1) stress balance equation establishment
The stress balance equation of the porous elastic rock is
Figure BDA0002386666510000021
In the formula: sigmaeffThe effective stress, MPa, can be calculated by formula (2); α is the Biot coefficient; i is the unit tensor, in the two-dimensional case [ 110]T(ii) a p is fluid pressure, MPa;
Figure BDA0002386666510000022
substituting equation (2) into equation (1), the stress balance equation can be rewritten as:
Figure BDA0002386666510000023
in the formula: ΨeffIs the elastic strain energy density, MPa, stored in the rock skeleton; epsilon is the strain tensor; g (c) is an attenuation function, and the attenuation function is defined as shown in formula (4); c is a fracture phase field, c is 1 to represent that the rock is completely fractured, and c is 0 to represent that the rock is intact; psi+ effThe tensile elastic strain energy density can be calculated by formula (5), and is MPa; psi- effThe compression elastic strain energy density can be calculated by formula (5), MPa; sigma+ effTensile stress, MPa; sigma- effCompressive stress, MPa.
g(c)＝(1-c)2 (4)
Figure BDA0002386666510000024
In the formula: λ and G are Lame constants, MPa; epsiloni(i ═ 1, 2, 3) as the main strain; function(s)<x>+＝(|x|+x)/2，<x>-＝(|x|-x)/2。
(2.2) stress balance equation corresponding to boundary conditions
The above-mentioned stress balance equation (3) can be solved by combining the boundary condition (6):
Figure BDA0002386666510000025
in the formula (I), the compound is shown in the specification,
Figure BDA0002386666510000026
as Dirichlet boundaries
Figure BDA0002386666510000027
Upper fixed displacement, MPa; t is the function on Neumann boundary
Figure BDA0002386666510000028
Stress above, MPa; n is Neumann boundary
Figure BDA0002386666510000029
The direction vector of (2).
The step (3) establishes a fluid flow control equation set, which comprises the following steps:
(3.1) establishment of equation of continuity of fluid flow in porous Medium
The equation for continuity of fluid flow in porous media is:
Figure BDA00023866665100000210
in the formula: t is time, s; ζ is the increment of the fluid volume content, which can be calculated by equation (8); v is the fluid flow velocity, m/s, which can be calculated by equation (9):
Figure BDA0002386666510000031
Figure BDA0002386666510000032
substituting equations (8) and (9) into equation (7), the fluid flow continuity equation is rewritten as:
Figure BDA0002386666510000033
in the formula: m is Biot modulus, MPa; epsiloniiVolume strain of the rock skeleton; k is the anisotropic permeability tensor; μ is the fluid viscosity, MPa.s.
(3.2) equation for permeability calculation
For the two-dimensional case, the permeability tensor can be expressed as:
Figure BDA0002386666510000034
wherein
Figure BDA0002386666510000035
In the formula: k is a radical ofxAnd kyPermeability in x and y directions, respectively, m2；k0Is the initial permeability of the rock matrix, m2；WcThe permeability weighting coefficient represents the contribution of the hydraulic fracture or the interlayer interface to the permeability of the calculation unit, and the invention adopts a simple permeability weighting coefficient calculation formula, namely Wc＝w/heW is the crack width, and since the cracks are transformed and distributed in the whole calculation region in the phase field method, the crack width of all the cells needs to be calculated, as shown in formula (13), heThe grid size of the finite element unit; k is a radical offIs hydraulic fracture or interlaminar interface permeability, m2Can be calculated by formula (14); θ is the normal direction angle or the maximum principal strain direction angle of the crack surface, and can be calculated by the formula (15).
w＝<ε1-εc>+he (13)
Figure BDA0002386666510000036
Figure BDA0002386666510000037
In the formula: epsilon1Is the maximum principal strain; eta is the shape parameter of the crack surface; gamma rayxyIs shear strain; epsilonyIs the y-direction strain; epsiloncThe critical tensile strain at the onset of rock fracture is represented by the critical tensile strain ε at the onset of rock fracture in the phase field methodcCritical tensile stress sigmacAnd fracture critical energy release rate GcThe three can be related by the formula (16):
Figure BDA0002386666510000038
in the formula: sigmacCritical tensile stress, MPa; gcThe crack critical energy release rate is MPa.m; l0Is a length scale parameter for controlling the diffusion crack region width, m, as shown in fig. 1; e is the rock elastic modulus, MPa.
(3.3) boundary conditions corresponding to the equation of continuity for fluid flow
The fluid flow continuity equation (10) can be solved in conjunction with the boundary conditions (17):
Figure BDA0002386666510000041
in the formula (I), the compound is shown in the specification,
Figure BDA0002386666510000042
to act on Dirichlet boundaries
Figure BDA0002386666510000043
Upper pressure, MPa; q is from Neumann boundary
Figure BDA0002386666510000044
Displacement of fracturing fluid injected upwards, m2/s。
The step (4) of establishing a fracture phase field evolution equation set comprises the following contents:
(4.1) phase field method approximation of sharp cracks
In the phase field method, the sharp crack Γ is converted into a diffusive crack Γ by an auxiliary phase field c (see fig. 1)c(c) And a diffusion crack gammac(c) Can be described by equation (18):
Figure BDA00023866665100000411
wherein the content of the first and second substances,
Figure BDA0002386666510000045
as a function of fracture density, as shown in equation (19):
Figure BDA0002386666510000046
(4.2) free energy Density build-up as Hydraulic fractures extend in porous media
When a hydraulic fracture extends in the porous medium, the total free energy density psi of the porous medium is determined by the elastic strain energy density psi stored in the rock skeletoneffEnergy density of reservoir in fluid ΨfluidEnergy density at break psifracIs composed of three parts, i.e.
Figure BDA0002386666510000047
Wherein:
Figure BDA0002386666510000048
Figure BDA0002386666510000049
Figure BDA00023866665100000410
(4.3) establishing a fracture phase field evolution equation based on variational principle
Substituting equations (21) - (23) into equation (20) yields an expression for the total free energy density Ψ. Furthermore, the evolution equation of the fracture phase field can be determined by the variation derivative of psi, and the evolution criterion can be obtained by the Kuhn-Tucker equation under the condition of irrelevant to the velocity, namely
Figure BDA0002386666510000051
Then the evolution equation of the fracture phase field is obtained as follows:
Figure BDA0002386666510000052
to satisfy the property of irreversible rock damage, the phase field evolution equation (25) can be rewritten as follows:
Figure BDA0002386666510000053
wherein, H (epsilon, t) is the maximum value of the tensile elastic strain energy density in the whole process and can be calculated by a formula (27):
Figure BDA0002386666510000054
(4.4) boundary conditions corresponding to the evolution equation of the phase field
The boundary conditions for the phase field evolution equation (26) are shown in equation (28):
Figure BDA0002386666510000055
in the formula:
Figure BDA0002386666510000056
to calculate the regional outer boundary.
Further, establishing a longitudinal extension numerical calculation model of the hydraulic fracture in the porous elastic medium by combining the steps (2) to (4); the method comprises the following steps:
equations (3), (10), (26) form a nonlinear system of equations, which is discretized by the finite element method, and the backward Euler method is used to discretize the time-dependent terms in equation (10). The crack phase field control equation (26) is solved by a Picard iteration method, and seepage-stress coupling equation sets (3) and (10) are solved by Newton-Raphson (NR) iteration. In each NR iteration step, the fracture phase field is a fixed value, and the NR iteration format of the percolation-stress coupling equation set can be written as:
Figure BDA0002386666510000057
in the formula: ruAnd RpThe margins for the stress balance equation and the fluid flow continuity equation, respectively, as shown in equations (30) and (31); j. the design is a squareuu、Jup、JpuAnd JppThe four components of the Jacobian matrix can be obtained by calculation according to a formula (32); delta uhIs the displacement increment of the ith iteration step, m; delta PhIs the pressure increment of the ith iteration step, MPa.
Figure BDA0002386666510000058
Figure BDA0002386666510000059
Figure BDA0002386666510000061
Figure BDA0002386666510000062
Figure BDA0002386666510000063
Figure BDA0002386666510000064
In equations (29) to (32), the superscript h represents the value at the computational grid node, and the subscript n represents the value at the nth time step; Δ t is the time step, s; n is a radical ofc、NuAnd NpRespectively carrying out finite element interpolation shape functions of a fracture phase field, displacement and pressure; b isuAnd Bu volRespectively strain matrix and volume strain matrix, BpThe matrix of the derivative of the shape function is interpolated for pressure.
The displacement increment delta u of the ith iteration step is obtained by equation (29)hAnd pressure increase deltaPhLater, the displacement and pressure for the (i + 1) th iteration can be expressed as:
Figure BDA0002386666510000065
furthermore, the strain energy history function H (epsilon, t) of the (i + 1) th iteration step can be obtained, and the heuristic solution of the phase field of the (i + 1) th iteration step can be obtained through an equation (34)
Figure BDA0002386666510000066
Wherein
Figure BDA0002386666510000067
Figure BDA0002386666510000068
And when the displacement, the pressure and the fracture phase field all meet the convergence condition shown in the formula (37), ending the iteration, and entering the calculation of the next time step, otherwise, continuing the iteration.
||Ru||≤tol||Ru0||，||Rp||≤tol||Rp0||，||ci+1-ci||≤tol||Rc0|| (37)。
The method adopts a phase field method and is established based on the minimization of system energy of a variation principle, so that the fracture extension track is automatically determined.
The fracture extension track and conditions are self-adaptive, the defect that the track prediction criterion needs to be additionally established to judge the fracture extension direction in the prior art is overcome, and the fluid loss coefficient does not need to be introduced to describe the fluid loss phenomenon of the fracturing fluid.
Drawings
FIG. 1 is a schematic representation of the conversion of sharp cracks into propagating cracks and boundary conditions for the examples;
FIG. 2 is a schematic diagram of a calculated zone and boundary conditions for coal seam perforation fracturing according to an embodiment;
FIG. 3 is a fracture phase field profile at the end of an example coal seam perforation fracture calculation;
FIG. 4 is a schematic diagram of a calculated zone and boundary conditions for top plate sandstone perforation fracturing in accordance with an embodiment;
figure 5 is a fracture phase field profile at the end of the example roof sandstone perforation fracture calculation.
Detailed Description
The invention will be described in further detail with reference to a certain well in Shanxi of China, but the invention is not limited to any one, wherein the formation parameters are shown in Table 1.
Table 1 table of formation base parameters used in the calculation of example 1
Figure BDA0002386666510000071
The first step is as follows: scheme 1: assuming that perforation fracturing is carried out in the center of a coal seam, a calculation area is shown in figure 2, the calculation area is uniformly divided into 80 x 80 square units, and a length scale parameter l00.5m, and a fracturing flow q of 2.2X 10-3m2And the viscosity of the fracturing fluid is 1mPa.s, the injection time is 36s, the time step length adopted by simulation is 3s, and the parameters in the table 1 are substituted into the equation set established by the invention for simulation. The simulation results are shown in fig. 3, and it was found that under this formation condition, the perforation fracture hydraulic fractures in the coal seam were unable to communicate between the upper and lower sandstone reservoirs.
The second step is that: scheme 2: because the perforation fracturing can not penetrate upper and lower sandstone layers in the coal bed, the perforation position is changed, the perforation fracturing is carried out in the top plate sandstone, the calculation area is shown in figure 4, the calculation area is uniformly divided into 80 multiplied by 80 square units, and the length scale parameter l00.5m, and a fracturing flow q of 2.2X 10-3m2And the viscosity of the fracturing fluid is 1mPa.s, the injection time is 36s, the time step length adopted by simulation is 3s, and the parameters in the table 1 are substituted into the equation set established by the invention for simulation. The simulation results are shown in fig. 5, and it can be seen from fig. 5 that when the roof is fractured by perforation, the hydraulic fracture can penetrate through the upper interface and continue to extend into the coal seam, but when the hydraulic fracture reaches the lower interface, the hydraulic fracture will extend along the interface and cannot break through the lower interface and extend to the bottom plate.
The third step: comparing scheme 1 and scheme 2, it can be seen that the reservoir thickness for scheme 2 was greater than for scheme 1, and therefore scheme 2 was selected for perforating.

Claims (1)

1. A coal-sand interbedded fracturing perforation position optimization method is characterized by comprising the following steps:
(1) collecting input parameters;
(2) establishing a stress calculation equation set;
(3) establishing a fluid flow control equation set;
(4) establishing a fracture phase field evolution equation set;
(5) establishing a longitudinal extension numerical calculation model of the hydraulic fracture in the porous elastic medium by combining the steps (2) to (4);
(6) inputting the parameters obtained in the step (1) into the model established in the step (5), and comparing the different formation parameters with the fracture trajectories under the conditions of the perforation positions, so as to optimize the perforation positions;
the input parameters collected in the step (1) comprise: the method comprises the following steps of (1) carrying out geostress parameters, elastic modulus and Poisson ratio of different layers of rock, critical tensile stress of different layers of rock, permeability of different layers of rock, critical tensile stress of interlayer interface, permeability of interlayer interface, fracturing discharge capacity, fracturing fluid injection time and fracturing fluid viscosity;
the step (2) establishes a stress calculation equation set, which comprises the following contents:
(2.1) stress balance equation establishment
The stress balance equation of a multi-elastic rock is as follows:
Figure FDA0002946451740000011
in the formula: sigmaeffThe effective stress is MPa and is obtained by calculation according to a formula (2); α is the Biot coefficient; i is the unit tensor, in the two-dimensional case [ 110]T(ii) a p is fluid pressure, MPa;
Figure FDA0002946451740000012
substituting equation (2) into equation (1), the stress balance equation can be rewritten as:
Figure FDA0002946451740000013
in the formula: ΨeffIs the elastic strain energy density, MPa, stored in the rock skeleton; epsilon is the strain tensor; g (c) is an attenuation function, and the attenuation function is defined as shown in formula (4); c is a fracture phase field, and c is 1 for rockComplete fracture, c ═ 0 means rock is intact; psi+ effThe tensile elastic strain energy density can be calculated by formula (5), and is MPa; psi- effIs the compressive elastic strain energy density, obtained by calculation according to formula (5), MPa; sigma+ effTensile stress, MPa; sigma- effCompressive stress, MPa;
g(c)＝(1-c)2 (4)
Figure FDA0002946451740000014
in the formula: λ and G are Lame constants, MPa; epsiloni(i ═ 1, 2, 3) as the main strain; function(s)<x>+＝(|x|+x)/2，<x>-＝(|x|-x)/2；
(2.2) stress balance equation corresponding to boundary conditions
The above-mentioned stress balance equation (3) can be solved by combining the boundary condition (6):
Figure FDA0002946451740000021
in the formula (I), the compound is shown in the specification,
Figure FDA0002946451740000022
as Dirichlet boundaries
Figure FDA0002946451740000023
Upper fixed displacement, MPa; t is the function on Neumann boundary
Figure FDA0002946451740000024
Stress above, MPa; n is Neumann boundary
Figure FDA0002946451740000025
The direction vector of (a);
the step (3) establishes a fluid flow control equation set, which comprises the following steps:
(3.1) establishment of equation of continuity of fluid flow in porous Medium
The equation for continuity of fluid flow in porous media is:
Figure FDA0002946451740000026
in the formula: t is time, s; ζ is the increment of the fluid volume content, which can be calculated by equation (8); v is the fluid flow velocity, m/s, which can be calculated by equation (9):
Figure FDA0002946451740000027
Figure FDA0002946451740000028
substituting equations (8) and (9) into equation (7), the fluid flow continuity equation is rewritten as:
Figure FDA0002946451740000029
in the formula: m is Biot modulus, MPa; epsiloniiVolume strain of the rock skeleton; k is the anisotropic permeability tensor; μ is the fluid viscosity, mpa.s;
(3.2) equation for permeability calculation
For the two-dimensional case, the permeability tensor can be expressed as:
Figure FDA00029464517400000210
wherein
Figure FDA00029464517400000211
In the formula: k is a radical ofxAnd kyPermeability in x and y directions, respectively, m2；k0Is the initial permeability of the rock matrix, m2；WcIs a permeability weighting coefficient representing the contribution of hydraulic fracture or interlayer interface to the permeability of the calculation unit, and adopts a permeability weighting coefficient calculation formula, namely Wc＝w/heW is the crack width, and since the cracks are transformed and distributed in the whole calculation region in the phase field method, the crack width of all the cells needs to be calculated, as shown in formula (13), heThe grid size of the finite element unit; k is a radical offIs hydraulic fracture or interlaminar interface permeability, m2Can be calculated by formula (14); θ is the normal direction angle or the maximum principal strain direction angle of the crack surface, and can be calculated by the formula (15):
w＝<ε1-εc>+he (13)
Figure FDA0002946451740000031
Figure FDA0002946451740000032
in the formula: epsilon1Is the maximum principal strain; eta is the shape parameter of the crack surface; gamma rayxyIs shear strain; epsilonyIs the y-direction strain; epsiloncThe critical tensile strain at the onset of rock fracture is represented by the critical tensile strain ε at the onset of rock fracture in the phase field methodcCritical tensile stress sigmacAnd fracture critical energy release rate GcThe three can be related by the formula (16):
Figure FDA0002946451740000033
in the formula: sigmacCritical tensile stress, MPa; gc is fracture critical energy release rateMPa.m; l0 is a length scale parameter for controlling the diffusion crack region width, m; e is the rock elastic modulus, MPa.
(3.3) boundary conditions corresponding to the equation of continuity for fluid flow
The fluid flow continuity equation (10) is solved in combination with the boundary conditions (17):
Figure FDA0002946451740000034
in the formula (I), the compound is shown in the specification,
Figure FDA0002946451740000035
to act on Dirichlet boundaries
Figure FDA0002946451740000036
Upper pressure, MPa; q is from Neumann boundary
Figure FDA0002946451740000037
Displacement of fracturing fluid injected upwards, m2/s；
The step (4) of establishing a fracture phase field evolution equation set comprises the following contents:
(4.1) phase field method approximation of sharp cracks
In the phase field method, the sharp crack Γ is passed through an auxiliary phase field c (converted into a diffusive crack Γ)c(c) And a diffusion crack gammac(c) Can be described by equation (18):
Figure FDA0002946451740000038
wherein the content of the first and second substances,
Figure FDA00029464517400000310
as a function of fracture density, as shown in equation (19):
Figure FDA0002946451740000039
(4.2) free energy Density build-up as Hydraulic fractures extend in porous media
When a hydraulic fracture extends in the porous medium, the total free energy density psi of the porous medium is determined by the elastic strain energy density psi stored in the rock skeletoneffEnergy density of reservoir in fluid ΨfluidEnergy density at break psifracIs composed of three parts, i.e.
Figure FDA0002946451740000041
Wherein:
Figure FDA0002946451740000042
Figure FDA0002946451740000043
Figure FDA0002946451740000044
(4.3) establishing a fracture phase field evolution equation based on variational principle
Substituting equations (21) - (23) into equation (20) can obtain an expression of the total free energy density Ψ; furthermore, the evolution equation of the fracture phase field can be determined by the variation derivative of psi, and under the condition of being independent of the velocity, the evolution criterion can be obtained by a Kuhn-Tucker equation, namely:
Figure FDA0002946451740000045
then the evolution equation of the fracture phase field is obtained as follows:
Figure FDA0002946451740000046
to satisfy the property of irreversible rock damage, the phase field evolution equation (25) can be rewritten as follows:
Figure FDA0002946451740000047
wherein, H (epsilon, t) is the maximum value of the tensile elastic strain energy density in the whole process and can be calculated by a formula (27):
Figure FDA0002946451740000048
(4.4) boundary conditions corresponding to the evolution equation of the phase field
The boundary conditions for the phase field evolution equation (26) are shown in equation (28):
Figure FDA0002946451740000049
in the formula:
Figure FDA00029464517400000410
to calculate the regional outer boundary;
establishing a longitudinal extension numerical calculation model of the hydraulic fracture in the porous elastic medium by combining the steps (2) to (4); the method comprises the following steps:
the equations (3), (10) and (26) form a nonlinear equation set, the nonlinear equation set is dispersed by adopting a finite element method, and a backward Euler method is adopted to disperse the terms related to time in the equation (10); solving a crack phase field control equation (26) by adopting a Picard iteration method, and iteratively solving seepage-stress coupling equation sets (3) and (10) by adopting Newton-Raphson (NR); in each NR iteration step, the fracture phase field is a fixed value, and the NR iteration format of the percolation-stress coupling equation set can be written as:
Figure FDA0002946451740000051
in the formula: ruAnd RpThe margins for the stress balance equation and the fluid flow continuity equation, respectively, as shown in equations (30) and (31); j. the design is a squareuu、Jup、JpuAnd JppThe four components of the Jacobian matrix can be obtained by calculation according to a formula (32); delta uhIs the displacement increment of the ith iteration step, m; delta PhIs the pressure increment of the ith iteration step, MPa;
Figure FDA0002946451740000052
Figure FDA0002946451740000053
Figure FDA0002946451740000054
Figure FDA0002946451740000055
Figure FDA0002946451740000056
Figure FDA0002946451740000057
in equations (29) to (32), the superscript h represents the value at the computational grid node, and the subscript n represents the value at the nth time step; Δ t is the time step, s; n is a radical ofc、NuAnd NpRespectively carrying out finite element interpolation shape functions of a fracture phase field, displacement and pressure; b isuAnd Bu volRespectively strain matrix and volume strain matrix, BpA derivative matrix that is a pressure interpolation shape function;
the displacement increment delta u of the ith iteration step is obtained by equation (29)hAnd pressure increase deltaPhLater, the displacement and pressure for the (i + 1) th iteration can be expressed as:
Figure FDA0002946451740000058
furthermore, the strain energy history function H (epsilon, t) of the (i + 1) th iteration step can be obtained, and the heuristic solution of the phase field of the (i + 1) th iteration step can be obtained through an equation (34)
Figure FDA0002946451740000059
Wherein
Figure FDA0002946451740000061
Figure FDA0002946451740000062
When the displacement, the pressure and the fracture phase field all meet the convergence condition shown in the formula (37), ending the iteration, and entering the calculation of the next time step, otherwise, continuing the iteration;
||Ru||≤tol||Ru0||，||Rp||≤tol||Rp0||，||ci+1-ci||≤tol||Rc0|| (37)。
CN202010100167.0A 2020-02-18 2020-02-18 Coal-sand interbedded through-layer fracturing perforation position optimization method Active CN111259595B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010100167.0A CN111259595B (en) 2020-02-18 2020-02-18 Coal-sand interbedded through-layer fracturing perforation position optimization method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010100167.0A CN111259595B (en) 2020-02-18 2020-02-18 Coal-sand interbedded through-layer fracturing perforation position optimization method

Publications (2)

Publication Number Publication Date
CN111259595A CN111259595A (en) 2020-06-09
CN111259595B true CN111259595B (en) 2021-04-27

Family

ID=70949326

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010100167.0A Active CN111259595B (en) 2020-02-18 2020-02-18 Coal-sand interbedded through-layer fracturing perforation position optimization method

Country Status (1)

Country Link
CN (1) CN111259595B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112036060B (en) * 2020-08-03 2022-04-01 武汉大学 Bilinear adaptive phase field method for simulating damage of brittle material
CN112487557B (en) * 2020-12-03 2022-10-14 上海交通大学 Method for predicting interface failure and microscopic crack propagation of composite material under hydraulic permeation load
CN116877067B (en) * 2023-07-18 2024-03-12 重庆地质矿产研究院 Method for predicting hydraulic fracturing generated cracks and swept area fluid pressure

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260543A (en) * 2015-10-19 2016-01-20 中国石油天然气股份有限公司 Double-pore model-based multi-medium oil gas flow simulation method and device
CN105807316A (en) * 2016-04-25 2016-07-27 吉林大学 Surface observation microseism speed model correcting method based on amplitude stack
CN108280275A (en) * 2018-01-09 2018-07-13 中国石油大学(华东) A kind of high prediction technique of tight sand hydraulic fracturing seam
CN108319756A (en) * 2017-12-29 2018-07-24 西安石油大学 A kind of compact reservoir volume fracturing seam net extended simulation and characterizing method
WO2019097326A1 (en) * 2017-11-17 2019-05-23 Universidad Pedagogica Y Tecnologica De Colombia Uptc Gasification of carbonaceous matter biomass mixture and coal using a cyclone-type forced flow furnace
CN110424939A (en) * 2019-08-12 2019-11-08 西南石油大学 A method of increasing gneiss oil-gas reservoir and stitches net volume fracturing effect

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105422071B (en) * 2015-12-07 2018-12-11 西南石油大学 Evaluate the rational method of low-permeable heterogeneous gas reservoir fracture parameters of fractured horizontal wells
CN108197358B (en) * 2017-12-20 2021-07-16 北京石油化工学院 Method for efficiently and quickly simulating hydraulic fracturing
JP7085261B2 (en) * 2018-01-25 2022-06-16 株式会社Ｉｈｉエアロスペース Manufacturing method of Hall thruster heater and Hall thruster heater
CN110298057B (en) * 2019-04-04 2022-04-05 西南石油大学 Supercritical carbon dioxide fracturing fracture extension calculation method

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260543A (en) * 2015-10-19 2016-01-20 中国石油天然气股份有限公司 Double-pore model-based multi-medium oil gas flow simulation method and device
CN105807316A (en) * 2016-04-25 2016-07-27 吉林大学 Surface observation microseism speed model correcting method based on amplitude stack
WO2019097326A1 (en) * 2017-11-17 2019-05-23 Universidad Pedagogica Y Tecnologica De Colombia Uptc Gasification of carbonaceous matter biomass mixture and coal using a cyclone-type forced flow furnace
CN108319756A (en) * 2017-12-29 2018-07-24 西安石油大学 A kind of compact reservoir volume fracturing seam net extended simulation and characterizing method
CN108280275A (en) * 2018-01-09 2018-07-13 中国石油大学(华东) A kind of high prediction technique of tight sand hydraulic fracturing seam
CN110424939A (en) * 2019-08-12 2019-11-08 西南石油大学 A method of increasing gneiss oil-gas reservoir and stitches net volume fracturing effect

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
临兴区块石盒子组致密砂岩气储层压裂缝高控制;吴锐 等;《煤炭学报》;20170915;第42卷(第9期);2393-2401 *

Also Published As

Publication number Publication date
CN111259595A (en) 2020-06-09

Similar Documents

Publication Publication Date Title
Wang Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using XFEM with cohesive zone method
Guo et al. Numerical investigation of hydraulic fracture propagation in a layered reservoir using the cohesive zone method
CN111259595B (en) Coal-sand interbedded through-layer fracturing perforation position optimization method
CA2739590C (en) Sand and fluid production and injection modeling methods
Wang et al. Comparison of consecutive and alternate hydraulic fracturing in horizontal wells using XFEM-based cohesive zone method
Sobhaniaragh et al. Three-dimensional investigation of multiple stage hydraulic fracturing in unconventional reservoirs
Ji et al. A novel hydraulic fracturing model fully coupled with geomechanics and reservoir simulation
Wang et al. Poroelastic and poroplastic modeling of hydraulic fracturing in brittle and ductile formations
CN111274731B (en) Fractured stratum fracturing fracture extension track prediction method
Zhu et al. Coupled flow-stress-damage simulation of deviated-wellbore fracturing in hard-rock
Feng et al. Parameters controlling pressure and fracture behaviors in field injectivity tests: a numerical investigation using coupled flow and geomechanics model
Feng et al. XFEM-based cohesive zone approach for modeling near-wellbore hydraulic fracture complexity
Cong et al. Numerical simulation of hydraulic fracture height layer-through propagation based on three-dimensional lattice method
Sobhaniaragh et al. The role of stress interference in hydraulic fracturing of horizontal wells
CN106285598A (en) A kind of shale seam net pressure break perforation cluster separation optimization method and system
Li et al. Modeling progressive breakouts in deviated wellbores
Liu et al. Well type and pattern optimization method based on fine numerical simulation in coal-bed methane reservoir
Li et al. Development of hydraulic fracture zone in heterogeneous material based on smeared crack method
Park et al. Rapid modeling of injection and production phases of hydraulically fractured shale wells using the fast marching method
Sobhaniaragh et al. Computational modelling of multi-stage hydraulic fractures under stress shadowing and intersecting with pre-existing natural fractures
Mukherjee et al. Successful control of fracture height growth by placement of artificial barrier
Liu et al. Numerical simulation of simultaneous multiple fractures initiation in unconventional reservoirs through injection control of horizontal well
Ma et al. Numerical simulation of progressive sand production of open-hole completion borehole in heterogeneous igneous formation
Rahmati et al. Validation of predicted cumulative sand and sand rate against physical-model test
CA2866156A1 (en) Screening potential geomechanical risks during waterflooding

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