A rapid simulation method for unconventional oil and gas reservoir productivity based on EDFM

Abstract

The present invention discloses a method for rapidly simulating the productivity of unconventional oil and gas reservoirs based on EDFM. Based on the embedded discrete fracture model and combined with the differential method, after dividing each hydraulic fracture into several micro-elements, the number of hydraulic fracture micro-elements in each matrix grid is divided according to the actual matrix grid size, and parameters such as the length of the hydraulic fracture and the average orthogonal distance in the grid are further calculated, so as to calculate the volume coefficient when crossflow occurs between the hydraulic fracture grid and the matrix grid. By coupling the continuous medium model, the deficiencies of the traditional embedded discrete fracture model are overcome. The embedded discrete fracture model is used to characterize hydraulic fractures, and the continuous medium model is used to characterize natural micro-fractures. The obtained model can accurately characterize various types of fractures in the reservoir, and the calculation results of the numerical model obtained in this way are more reasonable and reliable. The method of the present invention helps to study the seepage law in the formation of complex hydraulic fracture networks.

Description

Technical Field

[0001] The present invention relates to the technical field of unconventional oil and gas exploration and development, and in particular to a rapid productivity simulation method for unconventional oil and gas reservoirs based on EDFM. Background Art

[0002] With the increasing annual dependence on foreign oil and gas in China and the gradual depletion of conventional oil and gas resources, increasing the production of unconventional oil and gas reservoirs is crucial for China's energy strategic security. Unconventional oil and gas reservoirs are characterized by ultra-low porosity and ultra-low permeability, and generally use multi-stage fractured horizontal well technology to obtain industrial productivity. Given the complex seepage mechanism of unconventional oil and gas reservoirs, it is essential to use numerical simulation methods to conduct production simulation research on them.

[0003] Compared with the discrete fracture model, the embedded discrete fracture model, i.e., EDFM, does not require grid refinement when simulating hydraulic fractures, and only sacrifices a small computational accuracy in exchange for a significant improvement in computational speed, and is widely used in the simulation of unconventional oil and gas reservoirs. However, when using EDFM to characterize hydraulic fractures, there is a problem that it is difficult to quickly calculate the geometric parameters between the hydraulic fracture grid and the matrix grid. Currently, scholars still manually calculate the geometric parameters according to the position and length of each hydraulic fracture in the matrix grid. This method has a huge workload, brings inconvenience to related research, and makes the productivity simulation of oil and gas reservoirs time-consuming and laborious. Therefore, it is necessary to carry out relevant research work to solve the problem of rapid calculation of geometric parameters in EDFM and realize the rapid simulation of the productivity of unconventional oil and gas reservoirs. Summary of the Invention

[0004] In order to solve the problem of rapid calculation of geometric parameters in EDFM, the present invention provides a rapid productivity simulation method for unconventional oil and gas reservoirs based on EDFM.

[0005] The rapid productivity simulation method for unconventional oil and gas reservoirs based on EDFM provided by the present invention comprises the following steps:

[0006] S1: Obtain the porosity, permeability, oil saturation, gas saturation, and water saturation of the matrix and fractures, and obtain the relative permeability curve, capillary pressure curve, fluid PVT data, and rock compressibility.

[0007] S2: Derive the corresponding seepage equation according to the fluid components of the target reservoir, and use the fully implicit difference method for solution. The specific method is as follows:

[0008] The corresponding seepage equation derived according to the fluid components of the target reservoir is as follows:

[0009]

[0010] In the formula, M is the medium type, the matrix is m, the natural micro-fracture is f, and the hydraulic fracture is F; l is the fluid type, oil, gas, water; ql The fluid exchange volume between other medium type grids connected to this grid and this grid, m 3 / d; p l is the pressure of phase l, MPa; k M is the absolute permeability of medium of type M, D; K rl is the relative permeability of phase l, dimensionless; u l is the fluid viscosity of phase l, mPa·s; B l is the fluid volume coefficient of phase l, dimensionless; S l is the saturation of phase l, dimensionless; φ M is the porosity of medium of type M, dimensionless.

[0011] The derived seepage equation is solved using the fully implicit difference method. The process of the fully implicit difference method is as follows:

[0012] In the fully implicit method, for any variable X, the difference between its n-th time step and (n + 1)-th time step is defined as:

[0013] δX = X n+1 - X n

[0014] The difference between the v-th iteration step and the (v + 1)-th iteration step is:

[0015]

[0016] When v = 0, That is:

[0017] When the iteration meets the accuracy requirement, the result is obtained. Since and are too cumbersome to express, in this invention, hereinafter they are respectively replaced by X (v+1) and X (v) That is:

[0018]

[0019] Therefore, in the fully implicit method, when the iteration accuracy requirement is met, the value at the (n + 1)-th time step is approximately equal to the value at the (v + 1)-th iteration step. For example, during the simulation, when the pressure p at the (n + 1)-th time step satisfies the relation p n+1 - p (v+1) < ε, then there is:

[0020]

[0021] where: p (0) = p n

[0022] For the derived seepage equation, the results after processing the fluid exchange terms are as follows:

[0023] Oil phase:

[0024]

[0025] In the formula: q op is the volume of oil-phase fluid, m 3 ; p o1 , p o2 are the oil-phase pressures of the main grid and the grid adjacent to the main grid, MPa; S w1 , S w2 are the water saturations of the main grid and the grid adjacent to the main grid, dimensionless; S g1 , S g2 are the gas saturations of the main grid and the grid adjacent to the main grid, dimensionless; T op = G·f p (p o )·f s (s w , s g ) = G·f p ·f s , where G is a geometric parameter, f p is a function related to pressure, and f s is a function related to saturation.

[0026] Water phase:

[0027]

[0028] In the formula: q wp is the volume of water-phase fluid, m 3 ; p w1 , p w2 are the water-phase pressures of the main grid and the grid adjacent to the main grid, MPa; S w1 , S w2 are the water saturations of the main grid and the grid adjacent to the main grid, dimensionless; T wp = G·f p (p w )·f s (s w ) = G·f p ·f s .

[0029] Gas phase:

[0030]

[0031] In the formula: q gp is the volume of gas-phase fluid, m3 ; p g1 , p g2 are the gas-phase pressures of the main grid and the grid adjacent to the main grid, in MPa; S g1 , S g2 are the gas saturations of the main grid and the grid adjacent to the main grid, dimensionless;

[0032] T gp = G·f p (p g )·f s (s g ) = G·f p ·f s .

[0033] The result after processing the mass accumulation term of the derived seepage equation is as follows:

[0034] Oil phase:

[0035]

[0036] In the formula: B o is the formation volume factor of crude oil, dimensionless; S w is the water saturation, dimensionless; S g is the gas saturation, dimensionless; φ is the porosity, dimensionless.

[0037] Gas phase:

[0038]

[0039] In the above formula, 1 - s g - s w is used to replace the oil saturation.

[0040] In the formula: B g is the formation volume factor of formation gas, dimensionless; R s is the dissolved gas-oil ratio, dimensionless.

[0041] Water phase:

[0042]

[0043] In the formula: B w is the formation volume factor of formation water, dimensionless.

[0044] S3: Obtain the fracture height, fracture length, and azimuth angle of each hydraulic fracture, and determine the size of the target reservoir and the matrix grid size.

[0045] S4: Determine the differential unit length according to the matrix grid size, and calculate the calibration point coordinates based on the differential unit length, the end point position of the hydraulic fracture, and the azimuth angle of the hydraulic fracture.

[0046] Differential unit length d l The calculation formula is as follows:

[0047] d l =(d x +d y ) / 200

[0048] In the formula, d l is the differential unit length, m; d x is the length of the matrix grid in the x direction, m; d y is the length of the matrix grid in the y direction, m.

[0049] When determining the calibration points, starting from any endpoint of the hydraulic fracture as the starting calibration point, each point at intervals of one differential unit length along the hydraulic fracture is a calibration point, and the other endpoint of the hydraulic fracture is the last calibration point.

[0050] S5: Combine the matrix grid size divided in step S3 and the set of calibration point coordinates obtained in step S4 to determine the starting and ending coordinates and the fracture length of the hydraulic fracture in each matrix grid. Specifically, when operating, take the boundaries of each matrix grid as the boundaries for dividing the hydraulic fracture grid, and the starting calibration point and the ending calibration point in the matrix grid are the two endpoints of the hydraulic fracture in this matrix grid.

[0051] S6: Determine the intersection conditions of different hydraulic fractures according to the fracture length, azimuth angle and position in the model of each hydraulic fracture. When judging the intersection conditions of different fractures, determine whether the hydraulic fractures are parallel according to the azimuth angle; if the two fractures are not parallel, calculate the intersection coordinates of the two straight lines where the hydraulic fractures are located and determine whether this point is on the hydraulic fracture.

[0052] S7: Use the embedded discrete fracture model to characterize the hydraulic fracture, use the continuous medium model to characterize the natural fracture, and establish a numerical model.

[0053] The connection between the matrix and the hydraulic fracture grid for the geometric parameter G mF is as follows:

[0054]

[0055] In the formula, β c is the conductivity conversion factor, β c =86.4×10 -6 ; A mF is the area of the intersection surface between the hydraulic fracture grid and the matrix grid, m 2 ; k mF is the permeability of the connection pair between the fracture grid and the matrix grid, k mF ≈k m , D; d mF is the average orthogonal distance, the equivalent distance between the matrix grid and the fracture grid, m;

[0056] d mF The calculation formula of

[0057]

[0058] is as follows: In the formula, dv is the volume microelement within the grid block, m 3 ; x n is the normal distance from the microelement body to the fracture surface, m; V is the volume of the grid block, m 3 .

[0059] For the geometric parameter G of the adjacent grid connection pairs of the same hydraulic fracture F it is as follows:

[0060]

[0061] In the formula, are the permeabilities of the adjacent fracture grids 1 and fracture grid 2, D; A F is the area of the common surface of the two fracture grids, m 2 ; is the average distance from fracture grid 1 and fracture grid 2 to the common surface, m.

[0062] The calculation formula of

[0063]

[0064] is as follows: In the formula, ds i is the area microelement in fracture grid i (i = 1, 2), m 2 ; S i is the area of fracture grid i (i = 1, 2), m 2 ; x n is the distance from the microelement surface to the common surface, m.

[0065] For the geometric parameter G of the connection pairs when different hydraulic fractures intersect FF it is as follows

[0066]

[0067] In the formula, are the permeabilities of fracture grid 1 and fracture grid 2, D; ω F1 , ω F2 are the apertures of fracture grid 1 and fracture grid 2, m; L int is the length of the intersection line of the two fracture grids, m;

[0068] When the well grid is only connected to the hydraulic fracture grid, the well index calculation formula between the hydraulic fracture grid and the well grid is as follows:

[0069]

[0070]

[0071] where: w is the aperture of the artificial fracture, m; L F , h F are the length and height of the artificial fracture section, m; Δθ is the central angle of the radial well included in the fracture, rad; r w is the well radius, m

[0072] The calculation formula for the crossflow coefficient between the matrix grid and the natural fracture grid is as follows:

[0073]

[0074] where: L x , L y , L z are the dimensions of the matrix grid where the natural fracture is located in the x, y, and z directions, m.

[0075] S8: Use the numerical model established in step S7 to conduct production simulation on the unconventional oil and gas reservoir.

[0076] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0077] The rapid productivity simulation method for unconventional oil and gas reservoirs based on EDFM provided by the present invention is based on the embedded discrete fracture model and combines the differential method. After dividing each hydraulic fracture into several micro-elements, the number of hydraulic fracture micro-elements in each matrix grid is divided according to the actual matrix grid size, and parameters such as the length and average orthogonal distance of the hydraulic fracture in the grid are further calculated, so as to calculate the volume coefficient when crossflow occurs between the hydraulic fracture grid and the matrix grid. Compared with manually calculating the geometric parameters when crossflow occurs between the hydraulic fracture and the matrix grid, combined with the differential idea, a large number of calibration points are selected on each hydraulic fracture and the coordinates of the two end points of the hydraulic fracture in each matrix grid are determined according to the grid size, and then the parameters required for calculating the geometric parameters when crossflow occurs between the hydraulic fracture and the matrix are determined. Coupling the continuous medium model, it overcomes the deficiencies of the traditional embedded discrete fracture model. The embedded discrete fracture model is used to characterize the hydraulic fracture, and the continuous medium model is used to characterize the natural micro-fracture. The obtained model can accurately characterize various fractures in the reservoir, and the calculation results of the numerical model obtained in this way are more reasonable and reliable. The present invention helps to study the seepage law in the formation of complex hydraulic fracture networks.

[0078] Other advantages, objectives, and features of the present invention will be partially reflected by the following description and partially understood by those skilled in the art through the research and practice of the present invention. Description of the Drawings

[0079] Figure 1It is the capillary pressure curve and relative permeability curve.

[0080] Figure 2 It is the verification diagram of the correctness of the proposed method using commercial software. Specific implementation manners

[0081] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0082] Embodiment 1

[0083] The data used in this embodiment are from the public literature (Liu Lingfu, 2019). The relative permeability curve and capillary pressure curve of the target area are as Figure 1 shown, the formation crude oil properties are shown in Table 1, and the model parameter settings are shown in Table 2.

[0084] Table 1 Formation crude oil properties under different pressures

[0085]

[0086] Table 2 Model parameter settings

[0087]

[0088]

[0089] The two-phase oil-water seepage differential equation based on the coupled embedded discrete fracture model and continuous medium model is derived. The volume coefficient of each hydraulic fracture grid is calculated by the method described in the present invention, and the relevant numerical simulator is written using the Matlab computer language. The reservoir size is set to 1200m×1200m×10m (length×width×height), the grid size is 40m×40m×10m, and there is one hydraulic fracture located at the center of the model. The commercial software Saphir is used to verify the written numerical simulator, and the verification results are as Figure 2 shown. It can be seen that the model established by using the method of the present invention has a good fitting effect with the commercial software, indicating that the model obtained by the present invention can accurately characterize various fractures in the reservoir, and the calculation results of the obtained numerical model are reasonable and reliable.

[0090] In summary, the rapid productivity simulation method for unconventional oil and gas reservoirs based on EDFM is based on the embedded discrete fracture model and combines the differential method. After dividing each hydraulic fracture into several micro-elements, the number of hydraulic fracture micro-elements in each matrix grid is divided according to the actual matrix grid size, and parameters such as the length of the hydraulic fracture and the average orthogonal distance in the grid are further calculated, so as to calculate the volume coefficient when cross-flow occurs between the hydraulic fracture grid and the matrix grid. Coupling the continuous medium model overcomes the deficiencies of the traditional embedded discrete fracture model. The embedded discrete fracture model is used to characterize hydraulic fractures, and the continuous medium model is used to characterize natural micro-fractures. The obtained model can accurately characterize various fractures in the reservoir, and the calculation results of the numerical model obtained in this way are more reasonable and reliable. It helps to carry out the research on the seepage mechanism when developing unconventional oil and gas reservoirs using the multi-stage fractured horizontal well technology.

[0091] The above are only the preferred embodiments of the present invention and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to the equivalent embodiments by using the disclosed technical content within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM, characterized in that, the steps are as follows: S1: Obtain the porosity, permeability, oil saturation, gas saturation, and water saturation of the matrix and fractures, and obtain the relative permeability curve, capillary pressure curve, fluid PVT data, and rock compressibility; S2: Derive the corresponding seepage equation according to the fluid components of the target reservoir and solve it using the fully implicit difference method; S3: Obtain the height, length, and azimuth angle of each hydraulic fracture, and determine the size of the target reservoir and the matrix grid size; S4: Determine the differential unit length based on the matrix grid size, and the differential unit length is d l The calculation formula is as follows: d l = (d x + d y ) / 200 where d l is the differential unit length, m; d x is the length of the matrix grid in the x direction, m; d y is the length of the substrate grid in the y direction, m; Then, calculate the coordinates of the calibration points according to the differential unit length, the end point position of the hydraulic fracture, and the azimuth angle of the hydraulic fracture; when determining the calibration points, start from any end point of the hydraulic fracture as the starting calibration point, and each point at an interval of one differential unit length along the hydraulic fracture is a calibration point, and the other end point of the hydraulic fracture is the last calibration point; S5: Combine the matrix grid size divided in step S3 and the set of calibration point coordinates obtained in step S4 to determine the starting and ending coordinates and the length of the hydraulic fracture in each matrix grid; among them, take the boundary of each matrix grid as the boundary for dividing the hydraulic fracture grid, and the starting calibration point and the ending calibration point in the matrix grid are the two end points of the hydraulic fracture in this matrix grid; S6: Determine the intersection situation of different hydraulic fractures according to the length, azimuth angle, and position of each hydraulic fracture in the model; S7: Use the embedded discrete fracture model to characterize the hydraulic fracture, use the continuous medium model to characterize the natural fracture, and establish a numerical model; S8: Use the numerical model established in step S7 to conduct production simulation on the unconventional oil and gas reservoir.

2. The rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM according to claim 1, characterized in that, in step S6, when judging the intersection situation of different fractures, determine whether the hydraulic fractures are parallel according to the azimuth angle; if the two fractures are not parallel, calculate the intersection coordinates of the straight lines where the two hydraulic fractures are located and determine whether this point is on the hydraulic fracture.

3. The rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM according to claim 1, characterized in that, the corresponding seepage equation derived according to the fluid components of the target reservoir in step S2 is as follows: In the formula, M is the medium type, where the matrix is m, the natural microfracture is f, and the hydraulic fracture is F; l is the fluid type, including oil, gas, and water; q l is the fluid exchange volume between the grid of other medium types connected to this grid and this grid, m 3 / d; p l is the pressure of phase l, MPa; k M is the absolute permeability of medium M, D; K rl is the relative permeability of phase l, dimensionless; u l is the fluid viscosity of phase l, mPa·s; B l is the fluid volume coefficient of phase l, dimensionless; S l is the saturation of phase l, dimensionless; φ M is the porosity of medium M, dimensionless.

4. The rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM according to claim 3, characterized in that, for the seepage equation derived in step S2, the result after processing the fluid exchange term is as follows: Oil phase: where: q op is the volume of the oil-phase fluid, m 3 ; p o1 and p o2 are the oil phase pressures of the main grid and the grid adjacent to the main grid, in MPa; S w1 and S w2 are the water saturation of the main grid and the grid adjacent to the main grid, dimensionless; S g1 and S g2 are the gas saturation of the main grid and the grid adjacent to the main grid, dimensionless; T op = G·f p (p o )·f s (s w , s g ) = G·f p ·f s , where G is a geometric parameter, f p is a function related to pressure, and f s is a function related to saturation; Water phase: Where: q wp is the volume of the aqueous phase fluid, m 3 ; p w1 and p w2 are the aqueous phase pressures of the main grid and the grid adjacent to the main grid, in MPa; S w1 and S w2 are the water saturation of the main grid and the grid adjacent to the main grid, dimensionless; T wp = G·f p (p w )·f s (s w ) = G·f p ·f s ; Gas phase: Where: q gp is the volume of the gas-phase fluid, m 3 ; p g1 、p g2 are the gas-phase pressures of the main grid and the grid adjacent to the main grid, respectively, in MPa; S g1 、S g2 are the gas saturation of the main grid and the grid adjacent to the main grid, respectively, dimensionless; T gp = G·f p (p g )·f s (s g ) = G·f p ·f s .

5. The rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM according to claim 4, characterized in that, for the seepage equation derived in step S2, the result after processing the mass accumulation term is as follows: Oil phase: where: B o is the crude oil volume factor, dimensionless; S w is the water saturation, dimensionless; S g is the gas saturation, dimensionless; φ is the porosity, dimensionless; Gas phase: 1-s is used in the above formula g -s w to replace the oil saturation; Where: B g is the formation volume factor of the gas, dimensionless; R s is the solution gas-oil ratio, dimensionless; Water phase: Where: B w is the formation water volume factor, dimensionless.

6. The rapid simulation method for the productivity of unconventional oil and gas reservoirs based on EDFM according to claim 1, characterized in that, In the step S7, the connection between the matrix and the hydraulic fracture network to the geometric parameter G mF is as follows: where β c is the conductivity conversion factor, and β c = 86.4×10 -6 ; A mF is the area of the intersection surface between the hydraulic fracture network and the matrix network, m 2 ; k mF is the permeability of the fracture network and matrix network connection pair, D; d mF is the average orthogonal distance, the equivalent distance between the matrix network and the fracture network, m; d mF The calculation formula is as follows: where $dv$ is the volume element within the grid block, in $m$ 3 ; $x$ n is the normal distance from the element body to the fracture surface, in $m$; $V$ is the grid block volume, in $m$ 3 ; Geometric parameter G of adjacent grid connections in the same hydraulic fracture F As shown below: In the formula, is the permeability of adjacent fracture grids 1 and 2, D; A F is the common surface area of the two fracture grids, m 2 ; is the average distance from fracture grids 1 and 2 to the common surface, m; The calculation formula is as follows: where, ds i is the area microelement in crack network i, i = 1, 2, ..., m 2 ; S i is the area of the fracture grid i, where i = 1, 2, m 2 ; x n is the distance from the infinitesimal surface to the common surface, m; Connection of different hydraulic fractures when intersecting to geometric parameter G FF As shown below: In the formula, are the permeabilities of fracture grid 1 and fracture grid 2, D; are the apertures of fracture grid 1 and fracture grid 2, m; L int is the length of the intersection line of the two fracture grids, m; when the well grid is only connected to the hydraulic fracture grid, the well index calculation formula between the hydraulic fracture grid and the well grid is as follows: where: w is the aperture of the artificial fracture, m; L F , h F are the length and height of the artificial fracture section, m; Δθ is the central angle of the radial well contained in the fracture, rad; r w is the well radius, m; The calculation formula for the cross-flow coefficient between the matrix grid and the natural fracture grid is as follows: Among them: L x , L y , L z are the dimensions of the matrix grid where the natural fractures are located in the x, y, and z directions, in m.

Images