CN111259595B - Optimization method of perforation position for coal-sand interlayer penetrating fracturing - Google Patents

Abstract

The invention provides a method for optimizing the perforation position of coal-sand interlayer penetration fracturing, comprising the following steps: (1) collecting input parameters; (2) establishing a stress calculation equation group; (3) establishing a fluid flow control equation group; ( 4) Establish a fracture phase field evolution equation system; (5) Combine steps (2) to (4) to establish a numerical calculation model for the longitudinal extension of hydraulic fractures in poroelastic media; (6) Input the parameters obtained in step (1) into step (5) ) to compare the fracture trajectories under different formation parameters and perforation location conditions, so as to optimize the perforation location. The fracture propagation trajectory and conditions of the present invention are adaptively solved, which overcomes the defect of the prior art that an additional trajectory prediction criterion needs to be established to determine the fracture propagation direction, and the present invention does not need to introduce a fluid loss coefficient to describe the fluid loss of fracturing fluid Phenomenon.

Description

Optimization method of perforation position for coal-sand interlayer penetrating fracturing

technical field

The invention relates to the field of oil and gas field stimulation and reconstruction, in particular to a method for optimizing the perforation position of coal-sand interlayer penetration fracturing.

Background technique

China's coalbed methane fields are rich in natural gas resources. In addition to coalbed methane, this type of gas field also contains rich tight sandstone gas vertically, forming coal-sandstone interbeds, and hydraulic fracturing technology is used to connect coalstone and sandstone reservoirs vertically. It has become a way of efficient development of such gas fields. However, due to the influence of the interlayer interface and the difference in lithology and in-situ stress of different layers, if the perforation position is not appropriate, the hydraulic fracture can only extend in one layer. Layers are important.

The core of this technology is to establish a model of the longitudinal extension of hydraulic fractures in multiple layers. The commonly used numerical simulation methods for simulating the longitudinal extension of cracks in multiple layers include: (1) Palmer model and its improved model (2) Displacement discontinuity method (DDM); (3) Finite element method based on ABAQUS platform; (4) Extended finite element method. However, when the above numerical simulation method studies the extension trajectory of hydraulic fractures after encountering the interlayer interface, it is necessary to establish an intersection criterion to judge whether the hydraulic fractures pass through the interface or open the interface. In addition, most of the existing fracture propagation models need to introduce the Carter filtration coefficient to describe the fracturing fluid filtration phenomenon.

SUMMARY OF THE INVENTION

In view of the above technical problems, the present invention establishes a numerical calculation model for the longitudinal extension of hydraulic fractures in a poroelastic medium by comprehensively applying Biot's poroelasticity theory, finite element theory, phase field method theory, numerical solution method of nonlinear equations and other multidisciplinary knowledge. This model forms an optimization method for the perforation position of coal-sand interbed penetration fracturing.

A method for optimizing the perforation position of coal-sand interlayer penetration fracturing, comprising the following steps:

(1) Collect input parameters;

(2) Establish a set of stress calculation equations;

(3) Establish fluid flow control equations;

(4) Establish the evolution equations of the fracture phase field;

(5) Combining steps (2) to (4) to establish a numerical calculation model for the longitudinal extension of hydraulic fractures in a poroelastic medium;

(6) The parameters obtained in step (1) are input into the model established in step (5), and the fracture trajectories under different formation parameters and perforation position conditions are compared, so as to optimize the perforation position.

Further, the input parameters collected in the step (1) include: in-situ stress parameters, elastic modulus and Poisson's ratio of rocks in different layers, critical tensile stress of rocks in different layers, rock permeability of different layers, critical tensile stress of interlayer interface, Interlayer interface permeability, fracturing displacement, fracturing fluid injection time, fracturing fluid viscosity.

The step (2) establishes a set of stress calculation equations, including the following contents:

(2.1) Establishment of stress balance equation

The stress balance equation of poroelastic rock is

In the formula: σ _eff is the effective stress, MPa, which can be calculated by formula (2); α is the Biot coefficient; I is the unit tensor, which is [1 1 0] ^T in the two-dimensional case; p is the fluid pressure, MPa;

Bringing equation (2) into equation (1), the stress balance equation can be rewritten as:

In the formula: Ψ _eff is the elastic strain energy density stored in the rock skeleton, MPa; ε is the strain tensor; g(c) is the attenuation function, which is defined in the present invention as shown in formula (4); c is the fracture phase field, c=1 means the rock is completely broken, c=0 means the rock is intact; ψ ⁺ _eff is the tensile elastic strain energy density, which can be calculated by formula (5), MPa; ψ ^- _eff is the compressive elastic strain energy density, It can be calculated by formula (5), MPa; σ ⁺ _eff is tensile stress, MPa; σ ^- _eff is compressive stress, MPa.

g(c)=(1-c) ² (4)

In the formula: λ and G are the Lame constants, MPa; ε _i (i=1, 2, 3) is the principal strain; the function <x> ₊ = (|x|+x)/2, <x>-=( |x|-x)/2.

(2.2) The corresponding boundary conditions of the stress balance equation

The above stress balance equation (3) can be solved only by combining the boundary conditions (6):

为Dirichlet边界

上的固定位移，MPa；t为作用于Neumann边界

上的应力，MPa；n为Neumann边界

的方向向量。In the formula,

for the Dirichlet boundary

Fixed displacement on , MPa; t is acting on the Neumann boundary

stress on , MPa; n is the Neumann boundary

direction vector.

The step (3) establishes a fluid flow control equation system, including the following contents:

(3.1) Establishment of fluid flow continuity equation in porous media

The fluid flow continuity equation in porous media is:

Where: t is time, s; ζ is the increment of fluid volume content, which can be calculated by equation (8); v is fluid velocity, m/s, which can be calculated by equation (9):

Substituting equations (8) and (9) into equation (7), the fluid flow continuity equation is rewritten as:

where M is the Biot modulus, MPa; ε _ii is the volumetric strain of the rock skeleton; k is the anisotropic permeability tensor; μ is the fluid viscosity, MPa.s.

(3.2) Equation for calculating permeability

For the two-dimensional case, the penetration tensor can be expressed as:

in

where k _x and _ky are the permeability in the x and y directions, respectively, m ² ; k ₀ is the initial permeability of the rock matrix, m ² ; W _c is the permeability weighting coefficient, which represents the hydraulic fracture or interlayer interface pair. To calculate the contribution of the unit permeability, the present invention adopts a simple permeability weighting coefficient calculation formula, namely W _c = ^w /he , and w is the fracture width. Since the fractures are transformed and distributed in the entire calculation area in the phase field method, Therefore, it is necessary to calculate the fracture width of all elements, as shown in formula (13), h ^e is the mesh size of the finite element element; k _f is the hydraulic fracture or interlayer interface permeability, m ² , which can be calculated by formula (14) ; θ is the normal direction angle of the crack surface or the direction angle of the maximum principal strain, which can be calculated by formula (15).

w=<ε ₁ -ε _c > ₊ h ^e (13)

In the formula: ε ₁ is the maximum principal strain; η is the shape parameter of the fracture surface; γ _xy is the shear strain; ε _y is the y-direction strain; ε _c is the critical tensile strain when the rock begins to fracture. The critical tensile strain ε _c , the critical tensile stress σ _c and the critical energy release rate G _c of the crack can be related by formula (16):

Where: σ _c is the critical tensile stress, MPa; G _c is the critical energy release rate of the crack, MPa.m; l ₀ is the length scale parameter, which is used to control the width of the diffusion crack region, m, as shown in Fig. 1; E is Elastic modulus of rock, MPa.

(3.3) Boundary conditions corresponding to the fluid flow continuity equation

The above fluid flow continuity equation (10) can only be solved by combining the boundary conditions (17):

为作用于Dirichlet边界

上的压力，MPa；q为从Neumann边界

上注入的压裂液排量，m²/s。In the formula,

to act on the Dirichlet boundary

pressure on , MPa; q is from the Neumann boundary

Displacement of fracturing fluid injected above, m ² /s.

The step (4) establishes a fracture phase field evolution equation system, including the following contents:

(4.1) Approximation of sharp cracks by phase field method

In the phase field method, the sharp crack Γ is transformed into a diffusive crack Γ _c (c) by an auxiliary phase field c (see Fig. 1), and the diffusive crack Γ _c (c) can be described by Equation (18):

为裂缝密度函数，如公式(19)所示：in,

is the fracture density function, as shown in formula (19):

(4.2) Free energy density is established when hydraulic fractures extend in porous media

When a hydraulic fracture extends in a porous medium, the total free energy density Ψ of the porous medium consists of three parts: the elastic strain energy density Ψ _eff stored in the rock skeleton, the energy density of the reservoir in the fluid Ψ _fluid , and the fracture energy density Ψ _frac ,Right now

in:

(4.3) Establishment of 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 variational derivative of Ψ, and under the condition of being independent of the rate, the evolution criterion can be obtained by the Kuhn-Tucker equation, that is,

Then the evolution equation of the fracture phase field is:

In order to satisfy the irreversible property of rock damage, the phase field evolution equation (25) can be rewritten as follows:

Among them, H(ε, t) is the maximum value of tensile elastic strain energy density in the whole process, which can be calculated by formula (27):

(4.4) Boundary conditions corresponding to the phase field evolution equation

The boundary conditions corresponding to the phase field evolution equation (26) are shown in equation (28):

为计算区域外边界。where:

is the outer boundary of the calculation area.

Further, a numerical calculation model for the longitudinal extension of hydraulic fractures in a poroelastic medium is established in combination with steps (2) to (4), including the following contents:

Equations (3), (10) and (26) form a nonlinear equation system. The present invention uses the finite element method to discretize the nonlinear equation system, and uses the backward Euler method to discretize the time-related equations in equation (10). item. The fracture phase field governing equation (26) is solved by the Picard iteration method, and the seepage-stress coupling equations (3) and (10) are solved by the Newton–Raphson (NR) iteration. In each NR iteration step, the fracture phase field is a fixed value, then the NR iteration format of the seepage-stress coupling equation system can be written as:

where R _u and R _p are the remainders of the stress balance equation and the fluid flow continuity equation, respectively, as shown in equations (30) and (31); ^{Ju uu} , J ^up , J ^pu and J ^pp are the Jacobian matrices The four components of , can be calculated by formula (32); δu ^h is the displacement increment of the ith iteration step, m; δP ^h is the pressure increment of the ith iteration step, MPa.

In formulas (29) to (32), the superscript h represents the value at the node of the computing grid, and the subscript n represents the value at the nth time step; Δt is the time step, s; N _c , _Nu and N _p are the fracture phase field, displacement and pressure finite element interpolation shape functions, respectively; B _u and B _u ^vol are the strain matrix and volume strain matrix, respectively, and B _p is the derivative matrix of the pressure interpolation shape function.

After the displacement increment δu ^h and pressure increment δP ^h of the i-th iteration step are obtained by equation (29), the displacement and pressure of the i+1-th iteration step can be expressed as:

Furthermore, the strain energy history function H(ε, t) of the i+1-th iteration step can be obtained, and the tentative solution of the phase field of the i+1-th iteration step can be obtained by equation (34)

in

When the displacement, pressure and fracture phase field all satisfy the convergence conditions shown in formula (37), the iteration ends, and the calculation of the next time step is entered, otherwise the iteration continues.

||R _u ||≤tol||R _u0 ||, ||R _p ||≤tol||R _p0 ||, ||c _i+1 -c _i ||≤tol||R _c0 || ( 37).

The present invention adopts the phase field method to establish the system energy minimization based on the variational principle, so the crack extension trajectory is automatically determined.

The fracture propagation trajectory and conditions of the present invention are adaptively solved, which overcomes the defect of the prior art that an additional trajectory prediction criterion needs to be established to determine the fracture propagation direction, and the present invention does not need to introduce a fluid loss coefficient to describe the fluid loss of fracturing fluid Phenomenon.

Description of drawings

FIG. 1 is a schematic diagram of the conversion of sharp cracks into diffusive cracks and boundary conditions in an embodiment;

Fig. 2 is a schematic diagram of the calculation area and boundary conditions of coal seam perforation fracturing in an embodiment;

Fig. 3 is a fracture phase field distribution diagram at the end of the calculation of perforating fracturing in the coal seam of the embodiment;

4 is a schematic diagram of the calculation area and boundary conditions for perforating and fracturing of roof sandstone in an embodiment;

FIG. 5 is a distribution diagram of the fracture phase field at the end of the calculation of the perforation and fracturing of the roof sandstone in the embodiment.

Detailed ways

The present invention is further described in detail below in conjunction with a well in Shanxi, my country, but does not constitute any limitation to the invention, and the formation parameters are shown in Table 1.

Table 1 Basic parameters of the formation used in the calculation of Example 1

Step 1: Scheme 1: Suppose perforating fracturing is performed in the center of the coal seam. The calculation area is shown in Figure 2. The calculation area is evenly divided into 80 × 80 square units. The length scale parameter l ₀ is 0.5m. The fracturing displacement q is 2.2×10 ^-3 m ² /s, the fracturing fluid viscosity is 1 mPa.s, the injection time is 36s, the time step used in the simulation is 3s, and the parameters in Table 1 are brought into the present invention. The established equations are simulated. The simulation results are shown in Figure 3. It is found that under this stratum condition, the hydraulic fractures in the coal seam cannot communicate with the upper and lower sandstone reservoirs.

Step 2: Option 2: Since the perforation and fracturing in the coal seam cannot penetrate the upper and lower sandstone layers, the perforation position is changed, and the perforation and fracturing are performed in the roof sandstone. The calculation area is shown in Figure 4, and the calculation area is evenly distributed. is divided into 80×80 square units, the length scale parameter l ₀ is 0.5m, the fracturing displacement q is 2.2×10 ^-3 m ² /s, the fracturing fluid viscosity is 1mPa.s, the injection time is 36s, The time step used in the simulation is 3s, and the parameters in Table 1 are brought into the equation system established by the present invention for simulation. The simulation results are shown in Figure 5. It can be seen from Figure 5 that when the roof is perforated and fracturing, the hydraulic fracture can pass through the upper interface and enter the coal seam and continue to extend, but when the hydraulic fracture reaches the lower interface, the hydraulic fracture will extend along the interface. , it cannot break through the lower interface and extend to the bottom plate.

Step 3: Comparing scheme 1 and scheme 2, it can be seen that the thickness of the reservoir opened in scheme 2 is greater than that of scheme 1, so scheme 2 is selected for perforation.

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:

in the formula: sigma_effThe 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;

substituting equation (2) into equation (1), the stress balance equation can be rewritten as:

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)² (4)

in the formula: λ and G are Lame constants, MPa; epsilon_i(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):

in the formula,

as Dirichlet boundaries

Upper fixed displacement, MPa; t is the function on Neumann boundary

Stress above, MPa; n is Neumann boundary

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:

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):

substituting equations (8) and (9) into equation (7), the fluid flow continuity equation is rewritten as:

in the formula: m is Biot modulus, MPa; epsilon_iiVolume 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:

wherein

In the formula: k is a radical of_xAnd k_yPermeability in x and y directions, respectively, m²；k₀Is the initial permeability of the rock matrix, m²；W_cIs 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 W_c＝w/h^eW 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), h^eThe grid size of the finite element unit; k is a radical of_fIs hydraulic fracture or interlaminar interface permeability, m²Can 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＝<ε₁-ε_c>₊h^e (13)

in the formula: epsilon₁Is the maximum principal strain; eta is the shape parameter of the crack surface; gamma ray_xyIs shear strain; epsilon_yIs the y-direction strain; epsilon_cThe 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 method_cCritical tensile stress sigma_cAnd fracture critical energy release rate G_cThe three can be related by the formula (16):

in the formula: sigma_cCritical 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):

in the formula,

to act on Dirichlet boundaries

Upper pressure, MPa; q is from Neumann boundary

Displacement of fracturing fluid injected upwards, m²/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 gamma_c(c) Can be described by equation (18):

wherein,

as a function of fracture density, as shown in equation (19):

(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 skeleton_effEnergy density of reservoir in fluid Ψ_fluidEnergy density at break psi_fracIs composed of three parts, i.e.

Wherein:

(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:

then the evolution equation of the fracture phase field is obtained as follows:

to satisfy the property of irreversible rock damage, the phase field evolution equation (25) can be rewritten as follows:

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):

(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):

in the formula:

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:

in the formula: r_uAnd R_pThe 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 square^uu、J^up、J^puAnd J^ppThe four components of the Jacobian matrix can be obtained by calculation according to a formula (32); delta u^hIs the displacement increment of the ith iteration step, m; delta P^hIs the pressure increment of the ith iteration step, MPa;

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 of_c、N_uAnd N_pRespectively carrying out finite element interpolation shape functions of a fracture phase field, displacement and pressure; b is_uAnd B_u ^volRespectively strain matrix and volume strain matrix, B_pA 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 deltaP^hLater, the displacement and pressure for the (i + 1) th iteration can be expressed as:

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)

Wherein

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;

||R_u||≤tol||R_u0||，||R_p||≤tol||R_p0||，||c_i+1-c_i||≤tol||R_c0|| (37)。

Images