Numerical simulation method for low-permeability fault block oil reservoir based on embedded discrete fracture

Abstract

The invention provides a numerical simulation method for a low-permeability fault block oil reservoir based on embedded discrete fractures, which comprises the following steps of: S1, constructing basic units of an angular point grid: establishing a three-dimensional rectangular coordinate system by taking a central point of a matrix grid as a circle center, and dividing each vertex of the matrix grid into eight different quadrants according to coordinate values of each vertex, determining the position of each grid block through the eight angular points; s2, crack definition and intersection point calculation: defining crack positions and shapes through coordinate points, and performing grid conductivity calculation; and S3, constructing a low-speed nonlinear seepage mathematical model of the low-permeability fault block oil reservoir, and solving the model to obtain simulation parameters. According to the method, the complex geological structure of the fault block oil reservoir can be accurately described, in addition, the nonlinear seepage characteristics of a low-permeability reservoir matrix and various complex mechanisms of oil reservoir dynamic development are considered, and dynamic development of the low-permeability fault block oil reservoir can be accurately and efficiently simulated.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil and gas field development engineering, and particularly relates to a numerical simulation method for low-permeability fault block reservoirs based on embedded discrete fractures. BACKGROUND

[0002] With the increasing proportion of offshore low-permeability reserves in China year by year, the fracturing technology is continuously developed, so that the low-permeability reservoirs become the key to realize the goal of increasing reserves and production in the future. At present, the fault blocks of offshore reservoirs are relatively developed, and the conventional block center grid is difficult to depict complex structures. Therefore, it is of great significance to develop a numerical simulation method that can accurately describe complex geological structures and consider various complex factors of low-permeability reservoir seepage for the development and production dynamic prediction of offshore low-permeability fault block reservoirs. SUMMARY

[0003] In order to overcome the problems in the prior art, the present application provides a numerical simulation method for low-permeability fault block reservoirs based on embedded discrete fractures, to solve the problems in the prior art that the block center grid subdivision method cannot accurately describe complex fault block structures and does not fully consider the special seepage mechanism of low-permeability reservoirs. In order to achieve the above purpose, the present application provides the following solutions:

[0004] In a first aspect, the embodiments of the present application provide a numerical simulation method for low-permeability fault block reservoirs based on embedded discrete fractures, comprising the following steps:

[0005] S1: constructing a basic unit of an angle point grid: establishing a three-dimensional rectangular coordinate system with the center point of the matrix grid as the center, dividing each vertex of the matrix grid into eight different quadrants according to its coordinate value, and determining the position of each grid block through eight angle points;

[0006] S2: defining fractures and calculating intersection points: defining the position and shape of the fractures through coordinate points, and calculating the grid conductivity;

[0007] S3: constructing a low-speed nonlinear seepage mathematical model of a low-permeability fault block reservoir, and solving the model to obtain simulation parameters.

[0008] Optionally, the step S2 further comprises: judging whether the fracture vertex is in a certain matrix grid according to the distance between the fracture vertex and the center point of the matrix grid.

[0009] Optionally, the step S2 further comprises: if there is a deformed grid, determining the deformed angle point coordinates based on the eight quadrants, and fusing the deformed angle point coordinates.

[0010] Optionally, the grid conductivity calculation in the step S2 further comprises: calculating the half-conductivity and full-conductivity of the grid based on the two-point flow approximation method of the control volume.

[0011] Optionally, the total conductivity is determined by the harmonic average of the semiconductivity.

[0012] Optionally, the oil-water two-phase flow equation in the low-speed nonlinear seepage mathematical model in step S3 is:

[0013]

[0014] In the formula, v is the flow velocity; μ is the fluid viscosity; k is the permeability; p is the fluid pressure; G is the actual starting pressure gradient; δ rD δ is the thickness of the dimensionless residual boundary layer. rD =h s / r0, where h s The boundary layer thickness is given by r0 under a large pressure gradient, where r0 is the radius of the microtube; δ iD Let δ be the dimensionless initial boundary layer thickness. iD =h0 / r0, where h0 is the initial boundary layer thickness; It represents the low-velocity nonlinear coefficient in porous media.

[0015] Optionally, the low-velocity nonlinear seepage mathematical model in step S3 may consider one or more mechanisms, such as matrix initiation pressure gradient, nonlinear seepage, matrix stress sensitivity, or crack dynamic closure.

[0016] Optionally, step S3 may further include: solving the model to obtain simulation parameters and verifying them with commercial software.

[0017] Optionally, step S3 further includes:

[0018] The partial differential equations are discretized using automatic differentiation techniques to obtain the Jacobian matrix;

[0019] For each time step, the nonlinear equations are solved using Newton-Raphson iteration.

[0020] Optionally, the steps also include: dynamically updating the conductivity, permeability, and flow capacity based on the pressure field calculated at each time step.

[0021] Secondly, embodiments of the present invention also provide an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements a numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures as described in any of the embodiments of the present invention.

[0022] Thirdly, embodiments of the present invention also provide a storage medium for storing computer-executable instructions, which, when executed by a computer processor, are used to perform a numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures as described in any of the embodiments of the present invention.

[0023] Compared with the prior art, the advantages of the present invention are: the reservoir numerical simulation method provided by the present invention can accurately describe the complex geological structure of fault-block reservoirs. In addition, considering the nonlinear seepage characteristics of low-permeability reservoir matrix and various complex mechanisms of reservoir dynamic development, it can accurately and efficiently simulate the dynamic development of low-permeability fault-block reservoirs. Attached Figure Description

[0024] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0025] Figure 1 This is a schematic diagram of the method flow provided by the present invention;

[0026] Figure 2 Compare the vertex labels of the center grid and the corner grid;

[0027] Figure 3 This is a schematic diagram of the center point of the corner grid;

[0028] Figure 4 The algorithm is for "star-to-triangle" conversion.

[0029] Figure 5 The image shows the grid division of a single injection-production model for a low-permeability oil reservoir.

[0030] Figure 6 A model for the distribution of porosity (a) and permeability (b) in a positive rhythmic reservoir;

[0031] Figure 7 A comparison of daily water injection volume of injection wells and daily oil and water production of production wells;

[0032] Figure 8 Comparison of oil saturation fields in each sublayer;

[0033] Figure 9 This is the actual model mesh for a certain well area;

[0034] Figure 10 The actual model properties and pressure field diagram of a certain well area;

[0035] Figure 11 This refers to the oil pressure and daily oil production in a certain well area. Detailed Implementation

[0036] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0037] The purpose of this invention is to provide a numerical simulation method for low-permeability fault-block reservoirs based on embedded discrete fractures. To make the above-mentioned objectives, features and advantages of this invention more apparent and understandable, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] Reservoir numerical simulation technology is an important technique for studying the dynamic development of oil and gas reservoirs. In existing technologies, when studying fractured wells using the embedded discrete fracture modeling (EDFM) method, regular block-centered grids (such as...) are often employed. Figure 2 As shown in the figure, however, the block center grid is difficult to depict complex geological structures such as fault blocks and intercalation, and therefore cannot accurately simulate the dynamic development of oil reservoirs under complex geological conditions.

[0039] Furthermore, micron-sized pore throats are developed in low-permeability reservoirs, and the nonlinear permeability characteristics of the reservoir matrix cannot be ignored. Moreover, during the dynamic development of the reservoir, changes in formation pressure and stress conditions will lead to changes in the seepage characteristic parameters of the matrix and fractures, which also brings difficulties to the simulation of low-permeability reservoir development.

[0040] In order to solve the technical problems existing in the prior art, such as Figure 1 As shown, this invention provides a numerical simulation method for low-permeability fault-block reservoirs based on embedded discrete fractures, comprising the following steps:

[0041] S1: Basic unit for constructing corner grids: Establish a three-dimensional rectangular coordinate system with the center point of the matrix grid as the center, divide each vertex of the matrix grid into eight different quadrants according to its coordinate values, and determine the position of each grid block through the eight corner points;

[0042] S2: Crack Definition and Intersection Calculation: Define the crack location and shape using coordinate points, and calculate the mesh conductivity;

[0043] S3: Construct a low-velocity nonlinear seepage mathematical model for low-permeability fault-block reservoirs and solve the model to obtain simulation parameters.

[0044] In step S1, in order to accurately depict the complex geological structure of the fault block, a matrix model is established using the corner grid method. Specifically, a three-dimensional rectangular coordinate system is established with the center point of the matrix grid as the center, and each vertex of the matrix grid is divided into eight different quadrants according to its coordinate values. The position of each grid block is determined by the eight corner points.

[0045] In step S2, as Figure 3 As shown, the location and shape of the crack are defined by coordinate points, and the distance between the crack vertex and the center point of the matrix mesh is used to determine whether the crack vertex is in a certain matrix mesh.

[0046] In addition, the relative positional relationship between the crack vertex and each face of the mesh can be considered, as well as the angle between the vector direction of the vertex and the mesh boundary, in order to more accurately determine whether the crack vertex is truly located within a specific matrix mesh.

[0047] In calculating the grid conductivity, the two-point flow approximation (TPFA) method based on control volume is used to calculate the grid semi-conductivity and total conductivity.

[0048] In a preferred embodiment, the total conductivity T ij It is calculated using the harmonic average of the semiconductivity:

[0049]

[0050] In the formula, T i T is the semiconductivity from grid i to the interface ij; j Let be the semiconductivity from grid j to the interface ij.

[0051]

[0052] In the formula, k i m is the absolute permeability of grid cell i. 2 A ij m is the interface area between grid cells i and j. 2 ;d i Let m be the distance from the center of grid cell i to the center of the interface; n be the distance from the center of grid cell i to the center of the interface. i c is the unit normal vector pointing into the mesh element i on the interface; i Let i be the unit direction vector pointing from the grid cell i to the center of the interface.

[0053] In step S3, considering the low-velocity nonlinear seepage characteristics of the low-permeability pore throat, a seepage model is established.

[0054] Considering the multiphase flow characteristics in low-permeability reservoirs, the fluid flow in the reservoir satisfies the law of mass conservation, where equations (3) and (4) represent the mass conservation equations in the matrix and fractures, respectively.

[0055]

[0056] In the formula, t represents time, in seconds; and Porosity of the matrix and fractures, respectively; subscript κ indicates fluid type, i.e., oil, gas, or water; S κm and S κf These represent the saturation levels of the κ fluid in the matrix and the fracture, respectively; v κ,m and v κ,f These represent the flow velocities of oil or water in the matrix and cracks, respectively, in m·s. -1 ;q w κ,m and q w κ,f These represent the source and sink terms of the well in the matrix or fracture, respectively, m 3 ·s -1 ;q κ,mf m represents the fluid exchange term between the matrix and the fracture. 3 ·s -1 ;ρ κ The density of κ fluid is expressed in kg·m³. -3 .

[0057] In a preferred embodiment, the dynamic change of fluid density with pressure is considered, i.e., the fluid density ρ κ It is the fluid pressure p κ The function is shown in equation (5):

[0058]

[0059] In the formula, ρ κ,re f This indicates that the κ fluid is at a reference pressure p κ,ref Reference density under the given conditions, kg·m -3 ;p κ For κ fluid, the pressure is Pa; C κ The compressibility coefficient of the κ fluid is represented by Pa. -1 .

[0060] In a preferred embodiment, the nonlinear seepage mathematical model for the oil-water two-phase flow in the matrix is ​​as follows:

[0061]

[0062] In the formula, v is the flow velocity; μ is the fluid viscosity; k is the permeability; p is the fluid pressure; G is the actual starting pressure gradient; δ rD δ is the thickness of the dimensionless residual boundary layer. rD =h s / r0, where h s The boundary layer thickness is given by r0 under a large pressure gradient, where r0 is the radius of the microtube; δ iD Let δ be the dimensionless initial boundary layer thickness. iD =h0 / r0, where h0 is the initial boundary layer thickness; It represents the low-velocity nonlinear coefficient in porous media.

[0063] In a preferred embodiment, considering the dynamic changes in fracture conductivity, a dynamic attenuation model for fracture conductivity is constructed as follows:

[0064]

[0065] In the formula, F c For dynamic fracture conductivity; F ci d represents the initial conductivity of the fracture; d represents the average particle size of the proppant; ρ represents the sand concentration; p i This is the initial pressure.

[0066] In a preferred embodiment, considering the matrix stress-sensitive effect, the dynamic change model of matrix permeability is constructed as follows:

[0067]

[0068] In the formula, k m Indicates the dynamic permeability of the matrix; k mi γ represents the initial permeability of the matrix; γ is the matrix permeability stress sensitivity coefficient.

[0069] In a preferred embodiment, considering the dynamic change of fluid density with fluid pressure, the oil-water density change model is constructed as follows:

[0070]

[0071] In the formula, C represents the compressibility coefficient of oil or water; ρ β,re f This indicates that the β fluid is at the reference pressure p β,ref Reference density under the given conditions, kg·m -3 The density of the gas was obtained using the Peng-Robinson equation of state.

[0072] For source and sink terms q in the grid w κ,m and q w κ,f The calculation formulas are shown in equations (10) and (11), respectively.

[0073]

[0074] In the formula, p w The bottom hole pressure is in Pa. and These represent the pressure in the fracture or matrix mesh where the well is located, in Pa; WI m and WI f These are the well indices for wells falling within matrix grids or fracture grids, respectively, and m 3 .

[0075] The calculation formulas for the two well indices are shown in equations (12) and (13), respectively.

[0076]

[0077] In the formula, h is the reservoir thickness in meters (m); r w r is the radius of the wellbore, in meters (m). em and r ef where represents the equivalent radial flow radius (m) when the wellbore falls on the matrix or fracture grid, respectively. The calculation formula is shown below.

[0078]

[0079] In the formula, Δx and Δy represent the length and width of the matrix mesh, respectively, in meters; L and W represent the length and height of the crack mesh, respectively, in meters.

[0080] In a preferred embodiment, the mesh conductivity is corrected to take into account non-Darcy flow and stress sensitivity effects.

[0081] For the matrix mesh, non-Darcy flow and stress sensitivity effects are considered. Non-Darcy flow is characterized by equation (6), and the matrix stress sensitivity effect is characterized by an exponential equation. After correcting the half-conductivity of the matrix mesh, we obtain:

[0082]

[0083] In the formula, T' mi Non-Darcy half-conductivity; T mi γ is Darcy semiconductivity; γ is matrix stress sensitivity coefficient.

[0084] In a preferred embodiment, a fully implicit method is used in the solution process, where p is defined as an automatic differentiation variable. Each iteration of automatic differentiation at each time step reassembles the Jacobian matrix, thereby achieving a fully implicit solution at each time step.

[0085] Among them, the equivalent transformation from node-connected to node-unconnected in a crack network, such as Figure 4 As shown, the “star-triangle” conversion algorithm from the field of electrical engineering is introduced, and combined with the two-point flow approximation (TPFA) method based on control volume, a crack semi-conductivity correction model is constructed to address the dynamic influence of crack stress sensitivity on conductivity.

[0086]

[0087] In the formula, ρ is the sand concentration; d is the average particle size of the proppant; σ h For stress.

[0088] In a preferred embodiment, the partial differential equation is discretized using automatic differentiation technique to obtain the Jacobian matrix during the solution process; the operation rules of the automatic differentiation technique are shown in equation (18).

[0089]

[0090] In the formula, a and b are vectors.

[0091] After obtaining the Jacobian matrix, for each time step, the Newton iteration method is used to solve the nonlinear equation system. The Newton iteration format is shown in equation (19).

[0092]

[0093] In the formula, J is the Jacobian matrix, and V n V is an unknown variable at time n. n+1 The unknown variables at time n+1 include the pressure p and saturation S of each grid, as shown in Equation (20).

[0094]

[0095] In a preferred embodiment, the method further includes comparing the simulation parameters obtained by the method provided by the present invention with the simulation parameters of commercial software to verify the method provided by the present invention.

[0096] In an application instance, such as Figure 5 As shown, a numerical simulation model was established based on the actual conditions of the study block. The reservoir size is 500×500×50m, and the number of grids in each direction is 53×53×5. Two fractured vertical wells are symmetrically located at the left and right ends of the reservoir, with the production well on the left and the water injection well on the right, with a well spacing of 300m. Compared with onshore fracturing, offshore fracturing is smaller in scale, with displacement, fluid volume, and proppant addition amounting to approximately 20-30% of that used in onshore fracturing. Therefore, the half-length of the hydraulic fracturing fractures in both wells was set to 70m, the fracture aperture to 3mm, the fracture porosity to 0.5, and the fracture permeability to 200μm. 2 The fracture extends across the entire reservoir thickness, reaching a height of 50m.

[0097] To address the thin-layered characteristics of low-permeability reservoirs, a positive rhythmic reservoir model is constructed, such as... Figure 6 As shown, five sublayers, each 10m thick, were set up, with permeability ranging from 30×10⁻⁶ at the top to the bottom. -3 μm 2 40×10 -3 μm 2 50×10 -3 μm 2 60×10 -3 μm 2 70×10 -3μm 2 The porosity of the sublayers from top to bottom is 0.09, 0.095, 0.1, 0.105, and 0.11, respectively.

[0098] The original formation pressure of the reservoir was 40 MPa, and the densities of water, oil, and gas were 1000 kg / m³. 3 920kg / m 3 160kg / m 3 The viscosities of water, oil, and gas are 1 mPa·s, 5 mPa·s, and 0.01 mPa·s, respectively. The compressibility coefficient of the aqueous phase is 1 × 10⁻⁶. -9 Pa -1 The oil phase compressibility coefficient is 2×10. -9 Pa -1 The rock compressibility coefficient is 1×10 -10 Pa -1 The bound water saturation is 0.2, the residual oil saturation is 0.1, and the initial water cut saturation is 0.4. Water injection wells will be injected at a constant pressure of 42.5 MPa, and production wells will be produced at a constant bottomhole flowing pressure of 37.5 MPa. The predicted simulation period is 300 days.

[0099] like Figure 7 As shown, the simulation models described above were constructed using the EDFM method of this invention and the commercial software CMG-IMEX to predict production dynamics. It can be seen that both methods show good fit for the daily water injection volume of injection wells and the daily oil and water production of production wells. Furthermore, Figure 8 Comparing the saturation field diagrams of each sublayer after 300 days of development, EDFM also shows good consistency with the commercial software CMG-IMEX. In summary, this verifies the reliability of the method provided by this invention.

[0100] To address the complex reservoir characteristics of this block, the actual geological model was imported into a three-dimensional embedded discrete fracture model based on corner grids. The actual model mesh for the well area is shown below. Figure 9 As shown, the actual model properties and pressure field diagram of the well area are as follows: Figure 10 As shown, historical data fitting and production dynamics prediction were performed. For the actual well group, the first half was used for historical data fitting, and the second half was used for predicting simulation parameters. The results are as follows. Figure 11 As shown.

[0101] Therefore, the method provided by this invention uses a corner grid method to establish a matrix model, accurately describing the complex geological structure of fault-block reservoirs, and can accurately and efficiently simulate the dynamic development of low-permeability fault-block reservoirs. Furthermore, the method provided by this invention can also consider the nonlinear seepage characteristics of the low-permeability reservoir matrix, realizing the simulation of various complex mechanisms in the seepage simulation of low-permeability reservoirs, which is of great significance for the dynamic prediction of development and production of offshore low-permeability reservoirs.

[0102] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A numerical simulation method for low-permeability fault-block reservoirs based on embedded discrete fractures, characterized in that, Includes the following steps: S1: Basic unit for constructing corner grids: Establish a three-dimensional rectangular coordinate system with the center point of the matrix grid as the center, divide each vertex of the matrix grid into eight different quadrants according to its coordinate values, and determine the position of each grid block through the eight corner points; S2: Crack Definition and Intersection Calculation: Define the crack location and shape using coordinate points, and calculate the mesh conductivity; S3: Construct a low-velocity nonlinear seepage mathematical model for low-permeability fault-block reservoirs and solve the model to obtain simulation parameters.

2. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 1, wherein step S2 further includes: Determine whether a crack vertex is within a matrix mesh based on the distance between the crack vertex and the center point of the matrix mesh.

3. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 1, wherein step S2 further includes: If a distorted grid exists, the coordinates of the distorted corner points are determined based on the eight quadrants, and the coordinates of the distorted corner points are fused.

4. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 1, wherein the grid conductivity calculation in step S2 further includes: The two-point flow approximation method based on control volume is used to calculate the grid semi-transmissivity and total transmissivity.

5. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 4, wherein the total conductivity is determined by the harmonic average of the semi-conductivity.

6. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 1, wherein the oil-water two-phase flow equation in the low-velocity nonlinear seepage mathematical model in step S3 is: In the formula, v is the flow velocity; μ is the fluid viscosity; k is the permeability; p is the fluid pressure; G is the actual starting pressure gradient; δ rD δ is the thickness of the dimensionless residual boundary layer. rD =h s / r0, where h s The boundary layer thickness is given by r0 under a large pressure gradient, where r0 is the radius of the microtube; δ iD Let δ be the dimensionless initial boundary layer thickness. iD =h0 / r0, where h0 is the initial boundary layer thickness; It represents the low-velocity nonlinear coefficient in porous media.

7. The numerical simulation method for low-permeability faulted block reservoirs with embedded discrete fractures according to claim 1, wherein the low-velocity nonlinear seepage mathematical model in step S3 considers one or more mechanisms among matrix initiation pressure gradient, nonlinear seepage, matrix stress sensitivity, or fracture dynamic closure.

8. The numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures according to claim 1 further includes the step of: dynamically updating the conductivity, permeability, and flow capacity based on the pressure field calculated at each time step.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements a numerical simulation method for low-permeability fault-block reservoirs with embedded discrete fractures as described in any one of claims 1-8.

10. A storage medium storing computer-executable instructions, which, when executed by a computer processor, are used to perform a numerical simulation method for a low-permeability fault-block reservoir with embedded discrete fractures as described in any one of claims 1-8.

Images