Method and device for predicting gas-water two-phase productivity of fractured horizontal well

Abstract

The invention discloses a method and a device for predicting gas-water two-phase productivity of a fractured horizontal well. The method comprises the following steps: establishing a target block heterogeneous permeation Gaussian random field model and a target well reservoir geologic model; establishing a corresponding coarse scale grid, a fine scale grid and a dual control volume grid based on the target well reservoir geologic model; determining a dual basis function for mutual mapping of the coarse-scale grid and the fine-scale grid according to the coarse-scale grid, the fine-scale grid and the dual control volume grid; discretizing a control equation corresponding to the coarse-scale grid, and calculating pressure in cracks and a matrix on the coarse grid according to the heterogeneous permeation Gaussian random field model of the target block; according to the dual basis function, the pressure in the cracks and the matrix on the coarse grid is mapped to the fine-scale grid; and determining the gas-water two-phase productivity of the target well according to the cracks on the fine-scale grids and the pressure in the matrix.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical simulation of oil and gas field development, and in particular to a method and device for predicting gas-water two-phase productivity of fractured horizontal wells. BACKGROUND

[0002] Shale gas, as a kind of unconventional oil and gas reservoir, is abundant in reserves and widely distributed in China, and hydraulic fracturing has become the main technical means for shale gas development. Due to the very tight shale reservoirs, the permeability is low, so after fracturing, the well is usually closed for a period of time. Field practice has proved that the lower the fracturing fluid flowback rate after well closure, the higher the initial productivity of the gas well. Therefore, fracturing fluid flowback is a key link affecting the development effect of unconventional gas reservoirs, and therefore it is of great significance for the efficient development of shale gas in China to establish a reasonable numerical model for the flowback of fractured horizontal wells in shale gas reservoirs to predict the productivity of gas wells after flowback.

[0003] At present, there are many methods and models for simulating the flow rule of fluid in fractured reservoirs, the simulation methods include finite difference method (FDM), finite element method (FEM) and boundary element method (BEM), and the models mainly include equivalent medium, double porosity, discrete fracture and embedded discrete fracture, etc.

[0004] Due to the strong heterogeneity of shale gas reservoirs, the heterogeneity of the reservoir affects its various physical properties, and in the simulation calculation, the above-mentioned methods and models generally equivalent the shale reservoir to a homogeneous rock to reduce the number of grid division and the dimension of the calculation matrix. If you want to describe the heterogeneity of the reservoir through these models and methods, the grid in the numerical simulation area of the reservoir needs to be finely divided, at this time the grid unit can reach millions (the dimension of the calculation matrix is the square of the number of grids). At this time, the amount of calculation and time consumption is unacceptable in engineering. If a lower number of grids is used for calculation, the influence of heterogeneity on the flow rule of fluid in the reservoir cannot be completely captured.

[0005] Therefore, how to realize a multi-scale fractured horizontal well gas-water two-phase productivity prediction method considering the heterogeneity of shale reservoirs is a problem to be solved. SUMMARY

[0006] The present application provides a method and device for predicting the gas-water two-phase productivity of fractured horizontal wells, which combines a multi-scale finite volume method model (MSFV) to reduce the dimension of the fine grid system to a coarse grid system with low matrix dimension for solving, while ensuring the calculation accuracy and greatly saving the calculation amount and time consumption.

[0007] The present application provides a method for predicting the gas-water two-phase productivity of fractured horizontal wells, the method comprising:

[0008] establish a heterogeneous permeability Gaussian random field model of a target block and a reservoir geology model of a target well;

[0009] establish a coarse scale grid, a fine scale grid and a dual control volume grid corresponding to the reservoir geology model of the target well;

[0010] determine a dual basis function for mutual mapping of the coarse scale grid and the fine scale grid according to the coarse scale grid, the fine scale grid and the dual control volume grid;

[0011] discretize a control equation corresponding to the coarse scale grid, and calculate pressure in fractures and matrix on the coarse grid according to the heterogeneous permeability Gaussian random field model of the target block;

[0012] map the pressure in the fractures and the matrix on the coarse grid to the fine scale grid according to the dual basis function;

[0013] determine gas-water two-phase productivity of the target well according to the pressure in the fractures and the matrix on the fine scale grid.

[0014] The application also provides a device for predicting gas-water two-phase productivity of a fractured horizontal well, which comprises a memory and a processor; the memory is used to save a program for predicting gas-water two-phase productivity of a fractured horizontal well, and the processor is used to read and execute the program for predicting gas-water two-phase productivity of a fractured horizontal well, and execute the method in any one of the above embodiments.

[0015] The application also provides a computer readable storage medium, which stores a data processing program, and the data processing program is executed by a processor to execute the method in any one of the above embodiments.

[0016] Compared with the related art, the application provides a method for predicting gas-water two-phase productivity of a fractured horizontal well, which comprises: establishing a heterogeneous permeability Gaussian random field model of a target block and a reservoir geology model of a target well; establishing a coarse scale grid, a fine scale grid and a dual control volume grid corresponding to the reservoir geology model of the target well; determining a dual basis function for mutual mapping of the coarse scale grid and the fine scale grid according to the coarse scale grid, the fine scale grid and the dual control volume grid; discretizing a control equation corresponding to the coarse scale grid, and calculating pressure in fractures and matrix on the coarse grid according to the heterogeneous permeability Gaussian random field model of the target block; mapping the pressure in the fractures and the matrix on the coarse grid to the fine scale grid according to the dual basis function; and determining gas-water two-phase productivity of the target well according to the pressure in the fractures and the matrix on the fine scale grid. The application combines a multi-scale finite volume method model to reduce the dimension of the fine scale grid system to the dimension of the coarse grid system for solving, thereby greatly saving the calculation amount and time while ensuring the calculation accuracy.

[0017] Other features and advantages of the present application will be set forth in the description that follows, and in part will be apparent from the description, or can be learned by practice of the application. Other advantages of the application will be realized and attained by the methods particularly pointed out in the written description and claims hereof. BRIEF DESCRIPTION OF DRAWINGS

[0018] The accompanying drawings are included to provide a further understanding of the application and are incorporated in and constitute a part of this specification, illustrate embodiments of the application and are used to explain the technical solutions of the application, but do not constitute a limitation on the technical solutions of the application.

[0019] Figure 1 A flow chart of the method for predicting gas-water two-phase productivity of a fractured horizontal well according to an embodiment of the application;

[0020] Figure 2 A schematic diagram of the device for predicting gas-water two-phase productivity of a fractured horizontal well according to an embodiment of the application;

[0021] Figure 3 A schematic diagram of three grid system division in some exemplary embodiments;

[0022] Figure 4 A schematic diagram of projection calculation of physical quantities from a fine grid system to a coarse grid system in some exemplary embodiments;

[0023] Figure 5 A schematic diagram of four cases of interporosity flow coefficient between a fracture and a matrix in some exemplary embodiments;

[0024] Figure 6 A result of calculated gas production of a target shale gas well X1 using the method of the application in some exemplary embodiments;

[0025] Figure 7 A result of calculated gas production of a target shale gas well X1 using a conventional FVM in some exemplary embodiments;

[0026] Figure 8 Actual cumulative gas production data of a target shale gas well X1 in some exemplary embodiments. DETAILED DESCRIPTION

[0027] The present application describes a number of embodiments, but the description is exemplary rather than limiting and it will be apparent to those of ordinary skill in the art that numerous more embodiments and implementations are possible within the scope of the embodiments described in the present application. Although a number of possible combinations of features have been set forth in the accompanying figures and discussed above, many other combinations will be possible. Any feature of any embodiment can be used in combination with any other feature or element from any other embodiment or in place thereof, unless specifically restricted or un combinable. Any feature or element of any embodiment can be used in combination with features or elements from other inventive solutions to form another unique inventive solution defined by the claims.

[0028] The present application includes and contemplates combinations of features and elements known to those of ordinary skill in the art. The embodiments, features and elements disclosed herein can also be combined with any conventional features or elements to form unique inventive solutions defined by the claims. Any feature or element of any embodiment can also be combined with features or elements from other inventive solutions to form another unique inventive solution defined by the claims. Thus, it should be understood that any of the features shown and / or discussed in the present application can be implemented alone or in any appropriate combination. Embodiments are, therefore, not to be limited by any of the illustrative features set forth herein. Instead, they should be understood to include and be limited by only the appended claims, and equivalents thereof.

[0029] Furthermore, in describing representative embodiments, the specification can have presented the method and / or process as a particular sequence of steps. However, to the extent that the method or process depends on more than one step, the method or process should not be limited to the particular sequence of steps described. Other sequences of steps can be possible, and are within the scope of the present application. Therefore, the particular order in which the steps are presented is not limiting. Additionally, the specification can present the steps of the method and / or process in a particular order. However, the steps of the method and / or process can be performed in any order that is suitable to the application, as will be understood by those of ordinary skill in the art. Furthermore, any reference to claim limitations can not be construed as limiting the steps of the method and / or process to the particular order in which they are written.

[0030] The embodiments of the present application provide a method for predicting gas-water two-phase productivity of a fractured horizontal well, as shown in the following figure: Figure 1 The method comprises steps S100-S150:

[0031] S100: establishing a heterogeneous permeability Gaussian random field model of a target block and a reservoir geology model of a target well;

[0032] S110: establishing a corresponding coarse-scale grid, fine-scale grid and dual control volume grid based on the reservoir geology model of the target well;

[0033] S120: determining a dual basis function mapping the coarse scale grid and the fine scale grid according to the coarse scale grid, the fine scale grid and the dual control volume grid;

[0034] S130: discretizing the control equation corresponding to the coarse scale grid and calculating the pressure in the fracture and matrix on the coarse grid according to the target block heterogeneous permeability Gaussian random field model;

[0035] S140: mapping the pressure in the fracture and matrix on the coarse grid to the fine scale grid according to the dual basis function;

[0036] S150: determining the gas-water two-phase productivity of the target well according to the pressure in the fracture and matrix on the fine scale grid.

[0037] In an exemplary embodiment, the heterogeneous permeability Gaussian random field model is:

[0038]

[0039] In the above formula, is a gradient operator, φ1 is a first representation parameter, φ2 is a second representation parameter, and γ is a flux boundary condition is a velocity, and p is a pressure. In a specific study, φ1 and φ2 represent different meanings.

[0040] In an exemplary embodiment, the corresponding coarse scale grid, fine scale grid and dual control volume grid are established based on the target well reservoir geological model; as shown in Figure 3 is a schematic diagram of the division of the three grid systems of the coarse scale grid, fine scale grid and dual control volume grid. Specifically, the coarse scale and fine scale are divided according to the specific circumstances of the study area.

[0041] In an exemplary embodiment, according to the theory of multi-scale finite volume method, the physical quantity φ on the surface between adjacent nodes in the coarse scale grid system can be calculated by harmonic averaging. As shown in Figure 4 is a schematic diagram of the projection calculation of the physical quantity from the fine grid system to the coarse grid system. The specific algorithm is:

[0042] The calculation formula of the physical quantity on the grid surface in the x direction is:

[0043]

[0044] The above formula is defined as first performing harmonic mean processing along the x direction and then accumulating the sum along the y direction. Similarly, the calculation method for the mesh surface in the y direction is to first perform harmonic mean processing along the y direction and then accumulating the sum along the x direction. This φ is the permeability coefficient. In the above formula, the φ on the left side of the equal sign is the permeability coefficient of the coarse-scale grid system, and the φ on the right side is the permeability coefficient of the coarse-scale grid system. kl is the permeability coefficient of the fine-scale grid system. The above formula expresses the formula for mapping the permeability of the fine-scale grid system to the permeability of the coarse-scale grid system.

[0045] In an exemplary embodiment, since the heterogeneous permeability Gaussian random field model, i.e., equation (1), is solved on a coarse-scale grid system, it cannot fully reflect the heterogeneous characteristics of the medium. Therefore, it is necessary to construct its basis function in the dual space to ensure the flux conservation when the physical equation is mapped from the fine-scale grid to the coarse-scale grid. For any physical quantity φ, the dual control volume is:

[0046]

[0047] Where λ is the permeability, k = 1, 2, 3, 4 (equivalent to the number of grid nodes required for a grid in the coarse-scale grid system), is the dual control volume, is the kth dual basis function.

[0048] In an exemplary embodiment, dual basis functions for mapping the coarse-scale grid and the fine-scale grid are determined based on the coarse-scale grid, the fine-scale grid, and the dual control volume grid; specifically, the basis functions of each dual control volume are calculated; the physical field interpolation of the fine-scale grid is mapped to the coarse-scale grid system based on the basis functions; for example, the dual basis functions At the four nodes of the dual control volume there are:

[0049]

[0050] The mapping method is as follows: the dimensions of the basis functions at the four vertices of the dual control volume (that is, the four adjacent coarse-scale grid systems) are consistent with the number of fine-scale grids contained in the control volume. For example, if there are six fine-scale grids, then the basis functions for point p1 in the dual control volume are [p1x1, p1x2, p1x3, p1x4, p1x5, p1x6]. For example, if four basis functions are obtained, then the physical quantity mapping for fine-scale grid 1 is p1 (referring to the physical quantity of p1)*p1x1+p2*p2x2+p3*p3x3+p4*p4x4.

[0051] The dual basis function on the mesh boundary of the dual control volume is:

[0052]

[0053]

[0054]

[0055]

[0056] Where Γ1, Γ2, Γ3, Γ4 are four faces of the dual control volume; λ1 is the permeability in x direction, and λ2 is the permeability in y direction.

[0057] According to the formulas (4)-(7), the base functions of each dual control volume can be calculated, and the physical field of the fine scale grid is interpolated and mapped to the coarse scale grid system according to the calculated base functions.

[0058] In an exemplary embodiment, the process of discretizing the control equation corresponding to the coarse scale grid includes:

[0059] Firstly, the connectivity coefficient is determined;

[0060] Firstly, the intersection coordinates of the fracture network and the coarse scale grid lines are calculated;

[0061] Secondly, the fracture segments are reassembled according to the intersection coordinates, and the connectivity coefficient between the fracture segments and the matrix is calculated according to the fracture segments;

[0062]

[0063] Where A is the fracture area of the fracture segment in the matrix grid after the complete fracture is divided by the matrix grid, and d is the horizontal distance between the fracture segment and the center of the matrix grid.

[0064] In this step, in order to calculate the connectivity coefficient CI, the horizontal distance d between the fracture segment and the center of the matrix grid needs to be determined first. The calculation process of the horizontal distance d between the fracture segment and the center of the matrix grid is as follows:

[0065] Firstly, the intersection coordinates of the fracture network and the grid lines are calculated. Firstly, the coordinates of the grid system in x and y directions are put into the arrays mesh_X and mesh_Y respectively, and the start point and end point coordinates begin_xy and end_xy of a fracture F in the fracture network are extracted.

[0066] Secondly, according to the dimension d of the coordinate arrays of the grid system in x and y directions, the interpolation coordinate arrays F_X and F_Y corresponding to F in x and y directions (with the same dimension d) are established.

[0067] Third step, according to the x and y coordinates of the two points of begin_xy and end_xy to determine whether the crack is vertical or horizontal, if the x coordinates of begin_xy and end_xy are the same, the crack is vertical crack, then all the data in F_Y is set to NAN, if not vertical, then mesh_x is linearly interpolated to begin_xy and end_xy, and finally F_Y is obtained; if the y coordinates of begin_xy and end_xy are the same, the crack is horizontal crack, then all the data in F_X is set to NAN, if the crack is not horizontal, then mesh_y is linearly interpolated to begin_xy and end_xy, and finally F_X is obtained;

[0068] Fourth step, sequentially assemble column vectors A1: [F_X, mesh_X, begin_x, end_x] and A2: [mesh_Y, F_Y, begin_y, end_y] in order, and combine the two arrays A3 [A1, A2], if the two elements of a row of A3 are not NAN, then the row is the intersection point coordinates (including the starting point and the ending point) of the crack F after being divided by the grid system.

[0069] Fifth step, for each crack in the crack network system, the above first to fourth steps are performed respectively to obtain all the intersection point coordinates.

[0070] Sixth step, reassemble the crack segments according to the intersection point coordinates, and calculate the connectivity coefficient between the crack segments and the matrix according to the crack segments, as shown in Figure 5 The four cases of connectivity coefficient between the crack and the matrix, for each case, the corresponding crack segment and the horizontal distance d between the matrix grid center are calculated as follows:

[0071] d = dx / 4 (8)

[0072]

[0073]

[0074]

[0075] Based on the above formula, the connectivity coefficient CI between the crack and the matrix can be obtained, and then the flow exchange term between the crack and the matrix according to CI is obtained;

[0076] The difference between embedded discrete cracks and conventional reservoir percolation simulation mainly lies in the flow exchange term γ Fm between the crack and the matrix

[0077] γ Fm = CI λ m (p F -pm ) (12)

[0078] The main difference between embedded discrete fractures and conventional reservoir seepage simulation lies in the flow exchange term γ between the fractures and the matrix. mF , which is defined as:

[0079] γ mF =CIλ m (p m -p F ) (13)

[0080] In an exemplary embodiment, taking the fracture mesh discretization of single-phase flow as an example, the fracture control equation is discretized based on the flow exchange term between the fracture and the matrix calculated above. Specifically, the fracture control equation of single-phase flow is:

[0081]

[0082]

[0083] In the above formula, K m is the permeability of the matrix, K F is the permeability of the fracture, μ is the viscosity, p is the pressure, ρ is the density, g is the acceleration of gravity, γ mf is the flux from the crack to the matrix, γ Fm is the flux from matrix to fracture, q F is the flow rate in the crack, q m is the flow rate in the matrix.

[0084] In an exemplary embodiment, the original coarse-scale grid matrix is ​​A c m p c m =q c m (c represents the coarse-scale grid), combined with the embedded discrete crack grid model, the discrete grid control equations (14) and (15) are finite volume discretized, and the discretized matrix is ​​transformed into:

[0085]

[0086] Among them, A c ij represents the influence of j on i, ij = m, F (m is the matrix, F is the crack), A c EDFM represents the embedded discrete fracture coarse-scale grid matrix, p f EDFM represents the pressure of the embedded discrete crack, q f EDFMEDFM is a grid model representing embedded discrete fractures, and the flow rate of the embedded discrete fractures is a quantity that can be directly obtained. f EDFM The flow rate of the embedded discrete fractures is a quantity that can be directly obtained.

[0087] In an example embodiment, the pressure in the fractures and the matrix on the coarse grid is calculated according to the heterogeneous permeability Gaussian random field model of the target block; the pressure in the matrix and the fractures on the fine grid is mapped using the base function according to the obtained pressure in the matrix and the fractures on the coarse grid; and the productivity of the target well is determined according to the pressure in the fractures and the matrix on the fine grid.

[0088] In an example embodiment, for the productivity prediction process of the gas-water two-phase flow, the fracture grid discretization process of the single-phase flow is used, and the process is as follows:

[0089] Step S1: The heterogeneous permeability Gaussian random field model of the target block is corrected, and the permeability in the matrix system is:

[0090]

[0091] wherein, a is a dilution coefficient, k m0 is the absolute permeability of the shale matrix, b is a slippage coefficient, k n is the Knudsen number; and the above parameters can be directly obtained.

[0092] The permeability in the fracture system is:

[0093]

[0094] wherein, ζ is the main fracture network number, K ζ0 is the initial permeability, c ζ is a stress sensitivity coefficient, P ζ0 is the initial pressure; and the above parameters can be directly obtained.

[0095] Step S2: According to the corrected permeability coefficient, the water phase and the gas phase seepage equations are obtained, and the water phase seepage equation is:

[0096]

[0097] wherein, K Fw is the water phase relative permeability of the main fracture, μ w is the viscosity of the water phase, B w is the water phase volume coefficient, q Fw is the source-sink term of the water phase in the main fracture, Q fFw is the water phase channeling flow between the main fracture and the micro-fracture network grid, S Fw is the water phase saturation in the main fracture; and the above parameters can be directly obtained.

[0098] Step S3: The gas phase percolation equation is:

[0099]

[0100] In the formula, K Frg is the relative permeability of the main fracture gas phase, q Fg is the source and sink term of the gas term in the main fracture, Q fFg is the gas phase flow between the main fracture and the micro fracture network grid, S Fg is the gas phase saturation in the main fracture; each of the above parameters can be directly obtained.

[0101] Step S4: According to the initial physical field condition, the formulas of S1-S3 are discretized according to the process of single-phase fracture discretization, and the matrix Am of the matrix discretization and the matrix AF of the fracture discretization can be obtained under the coarse grid system. Am and AF are assembled according to the embedded discrete fracture grid model, and the pressure in the matrix and the fracture under the initial physical field condition of the grid system = fracture grid + matrix grid (the matrix grid here is a coarse grid) can be obtained.

[0102] Step S5: According to the obtained pressure in the matrix and the fracture under the coarse grid, the base function is used to map the pressure in the matrix and the fracture under the fine grid.

[0103] The method for predicting the gas-water two-phase productivity of the fractured horizontal well realized by the embodiment has the following technical effects: by considering the influence of the strong heterogeneity of the shale reservoir on the distribution of the permeability, the pressure and the saturation in the reservoir, combining the multi-scale finite volume method model (MSFV) to reduce the dimension of the fine grid system with high matrix dimension to the coarse grid system with low matrix dimension for solving, the calculation amount and the time consumption are greatly saved.

[0104] Example one

[0105] In order to solve the problems of large calculation amount and difficult simulation of gas-water two-term flow in the strong heterogeneity shale reservoir, a method for predicting the gas-water two-phase productivity of the fractured horizontal well is proposed, and the specific implementation process includes the following steps:

[0106] Step S10: Selecting the reservoir shale core of the target shale gas well X1, cutting the core into several rock pieces and measuring the permeability, and establishing a heterogeneous permeability Gaussian random field model of the target block according to the rock piece permeability test result;

[0107] Step S20: Establish a reservoir geological model of the target shale gas well X1, and establish a coarse-scale grid, a fine-scale grid, and a dual control volume grid of the reservoir based on the target shale gas well X1 reservoir geological model; and an initial fracture grid; the reservoir geological model here is a model established based on the fracture parameters provided on site, such as: number of fractures, location, azimuth, inclination angle, fracture width, etc.

[0108] In step S20, a fine-scale grid system is established for the numerical simulation area, with the number of grids being M×N; a coarse-scale grid system is established, with the number of grids being m×n, and a dual control volume is established based on the coarse-scale grid system.

[0109] In this step, coarse-scale and fine-scale grid systems are two types of grid systems within a multi-scale grid. For example, in the same area, a coarse-scale grid system is defined as 10,000 grids, while a fine-scale grid system is defined as 10,000 grids. The specific coarse-scale and fine-scale standards are determined based on actual conditions.

[0110] Step S30: Calculate the mapping of the fine-scale inhomogeneous physical field on the coarse-scale inhomogeneous field, and calculate the dual basis function;

[0111] The specific process of step S30 is:

[0112] Step S301: Initialize the physical field φ at a fine scale. According to the seepage equation, the mathematical model of seepage flow in the reservoir is an elliptic equation, which is a heterogeneous permeability Gaussian random field model, specifically:

[0113]

[0114] According to the theory of multi-scale finite volume method, the physical quantity φ on the surface between adjacent nodes in the coarse-scale grid system can be calculated by harmonic averaging. Specifically, the calculation formula of the physical quantity on the grid surface in the x direction is:

[0115]

[0116] The above formula is defined as first performing harmonic mean processing along the x direction and then accumulating the sum along the y direction. Similarly, the calculation method for the mesh surface in the y direction is to first perform harmonic mean processing along the y direction and then accumulating the sum along the x direction. This φ is the permeability coefficient. In the above formula, the φ on the left side of the equal sign is the permeability coefficient of the coarse-scale grid system, and the φ on the right side is the permeability coefficient of the coarse-scale grid system. kl is the permeability coefficient of the fine-scale grid system. The above formula expresses the formula for mapping the permeability of the fine-scale grid system to the permeability of the coarse-scale grid system, as shown in Figure 4 As shown in the figure, a schematic diagram of mapping the permeability of a fine-scale grid system to the permeability of a coarse-scale grid system is shown.

[0117] Step S302, since equation (1) is solved on the coarse-scale grid system, the heterogeneous characteristics of the medium cannot be fully reflected, therefore, it is necessary to construct the basis function thereof on the dual space to ensure the flux conservation of the physical equation when mapping from the fine-scale grid to the coarse-scale grid. For any physical quantity φ:

[0118]

[0119] wherein k = 1, 2, 3, 4, the number of k represents the number of grid nodes required for one grid in the coarse-scale grid system, is the dual control volume; λ is the permeability, is the gradient operator.

[0120] Step S303, calculating the basis function of each dual control volume; according to the basis function, the physical field of the fine-scale grid is interpolated and mapped to the coarse-scale grid system;

[0121] For example: dual basis function On the four nodes of the dual control volume:

[0122]

[0123] The mapping method is specifically: the dimension of the basis function of the four vertices of the dual control volume (that is, the four adjacent coarse-scale grid systems) is consistent with the number of fine-scale grids contained by the control volume. For example: if there are 6 fine-scale grids, then the basis function of the p1 point of the dual control volume is

p1x1, p1x2, p1x3, p1x4, p1x5, p1x6

[0124] Step S303, the dual basis function on the grid boundary of the dual control volume:

[0125]

[0126]

[0127]

[0128]

[0129] wherein Γ1, Γ2, Γ3, Γ4 are the four faces of the dual control volume;

[0130] According to the formulas (4)-(7), the base functions of each dual control volume can be calculated, and the physical field of the fine scale grid is interpolated and mapped to the coarse scale grid system according to the calculated base functions.

[0131] Step S40: calculating the intersection points and connectivity coefficients of the fracture grid and the coarse scale grid, and initializing the fracture coefficient matrix;

[0132] In this step S40, the specific implementation process includes:

[0133] Step S401: calculating the intersection point coordinates of the fracture network and the coarse scale grid lines;

[0134] In this step, first, the coordinates of the grid system in the x and y directions are put into the arrays mesh_X and mesh_Y, respectively; wherein the grid system = fracture grid + matrix grid (here, the matrix grid is the coarse scale grid).

[0135] Secondly, the start point and end point coordinates begin_xy and end_xy of a fracture F in the fracture network are extracted.

[0136] Step S402: according to the dimension d of the coordinate arrays of the grid system in the x and y directions, the interpolation coordinate arrays F_X and F_Y (both with the dimension d) corresponding to F in the x and y directions are established.

[0137] Step S403: judging whether the fracture is vertical or horizontal according to whether the x and y coordinates of the two points begin_xy and end_xy are the same;

[0138] If the x coordinates of begin_xy and end_xy are the same, the fracture is a vertical fracture, then all the data in F_Y are set to NAN, if not vertical, then mesh_x is linearly interpolated to begin_xy and end_xy, and finally F_Y is obtained;

[0139] If the y coordinates of begin_xy and end_xy are the same, the fracture is a horizontal fracture, then all the data in F_X are set to NAN, if the fracture is not horizontal, then mesh_y is linearly interpolated to begin_xy and end_xy, and finally F_X is obtained;

[0140] Step S404: sequentially assembling column vectors A1: [F_X, mesh_X, begin_x, end_x] and A2: [mesh_Y, F_Y, begin_y, end_y], and combining the two arrays A3 [A1, A2], if the two elements of a row of A3 are not NAN, the row is the intersection point coordinates (including the start point and end point) of the fracture F after being divided by the grid system, in this step, under the condition that all are not NAN, the row is the intersection point coordinates.

[0141] Step S405: performing steps S401-S404 on each fracture in the fracture network system to obtain all intersection coordinates.

[0142] Step S406: reassembling the fracture segments according to the intersection coordinates and calculating the connectivity coefficient between the fracture segments and the matrix according to the fracture segments.

[0143] At this time, the connectivity coefficient is divided into four cases, and the analytical solution calculation formula of each case is as follows:

[0144] d = dx / 4 (8)

[0145]

[0146]

[0147]

[0148] In the above four formulas, d is the distance between the fracture segment and the center of the matrix grid to which the fracture segment belongs.

[0149] According to the distance d between the fracture segment and the center of the matrix grid to which the fracture segment belongs calculated in formulas (8)-(11), the connectivity coefficient CI between the fracture and the matrix is determined according to the topological relationship between CI and the fracture and the matrix.

[0150]

[0151] where A is the fracture area of the complete fracture after being divided by the matrix grid, which is a directly obtainable parameter, and d is the distance between the fracture segment and the center of the matrix grid to which the fracture segment belongs.

[0152] Step S407: discretizing the control equation to obtain an embedded discrete fracture coarse-scale grid matrix.

[0153] Step S4071: the connectivity coefficient CI between the fracture and the matrix is calculated based on the above formula, and then the flow exchange term between the fracture and the matrix is calculated according to CI.

[0154] The difference between the embedded discrete fracture and the conventional reservoir percolation simulation mainly lies in the flow exchange term γ between the fracture and the matrix Fm , which is defined as:

[0155] γ Fm = CIλ m (p F -p m ) (12)

[0156] The difference between the embedded discrete fracture and the conventional reservoir percolation simulation is mainly in the flow exchange term γ between the fractures and the matrix mF , which is defined as:

[0157] γ mF = C I λ m (p m -p F ) (13)

[0158] Step S4072: The fracture grid control equation is (taking single-phase flow as an example):

[0159]

[0160]

[0161] In the above formula, K is the permeability of the matrix, K is the permeability of the fracture, μ is the viscosity, p is the pressure, ρ is the density, g is the gravitational acceleration, γ is the flow exchange from the fracture to the matrix, γ is the flow exchange from the matrix to the fracture, q is the flow in the fracture, and q is the flow in the matrix. m F mf Fm F m

[0162] Step S4073: The original coarse-scale grid matrix is A c m p c m = q c m (c represents the coarse-scale grid), combined with the embedded discrete fracture grid model, the discrete grid control equations (14) and (15) are discretized by finite volume, and the discretized matrix is transformed into:

[0163]

[0164] Wherein, A c ij represents the influence of j on i, ij = m, F (m is the matrix, and F is the fracture), A c EDFM represents the embedded discrete fracture coarse-scale grid matrix, p f EDFM represents the pressure of the embedded discrete fracture, q f EDFM represents the flow of the embedded discrete fracture, and EDFM is the embedded discrete fracture grid model. q f EDFM represents the flow of the embedded discrete fracture, which is a quantity that can be directly obtained.

[0165] ​​​​​​Step S50: based on the discretized coarse-scale control grid matrix obtained in the above steps, solve the pressure in the fracture and matrix on the coarse-scale grid;

[0166] Further technical solutions, the specific process of step S50 is:

[0167] Step S501: correct the target block heterogeneous permeability Gaussian random field model, the permeability in the matrix system is:

[0168]

[0169] In the formula: alpha is a dilute coefficient, k m0 is the absolute permeability of shale matrix, b is the slip coefficient, k n is the Knudsen number; the above parameters can be directly obtained.

[0170] The permeability in the fracture system is:

[0171]

[0172] In the formula: zeta is the main fracture network number, K ζ0 is the initial permeability, c ζ is the stress sensitivity coefficient, P ζ0 is the initial pressure; the above parameters can be directly obtained.

[0173] Step S502: according to the corrected permeability coefficient, the water phase and gas phase seepage equation can be obtained, the water phase seepage equation is:

[0174]

[0175] In the formula: K Fw is the relative permeability of water phase in the main fracture, mu w is the viscosity of water phase, B w is the volume coefficient of water phase, q Fw is the source and sink term of water in the main fracture, Q fFw is the water phase channeling flow between the main fracture and the micro fracture network grid, S Fw is the water phase saturation in the main fracture; the above parameters can be directly obtained.

[0176] Step S503: the gas phase seepage equation is:

[0177]

[0178] In the formula: K Frg is the relative permeability of gas phase in the main fracture, q Fg is the source and sink term of gas in the main fracture, Q fFg is the gas phase channeling flow between the main fracture and the micro fracture network grid, S Fgis the gas saturation in the main fracture; each of the above parameters can be directly obtained.

[0179] Step S504: According to the initial physical field condition, the formula of S501-S503 is discretized in the manner of S20-S40 to obtain the matrix discrete matrix Am and the fracture discrete matrix AF under the coarse scale grid system, Am and AF are assembled according to the embedded discrete fracture grid model, and the pressure in the matrix and the fracture under the initial physical field condition can be obtained. The grid system = fracture grid + matrix grid (the matrix grid here is a coarse scale grid).

[0180] Step S505: According to the obtained pressure in the matrix and the fracture under the coarse scale grid, the base function is used to map the pressure in the matrix and the fracture under the fine scale grid.

[0181] Step S60: Repeat steps S30-S50 within the calculation time until the calculation time ends, and draw a time and gas production relationship curve.

[0182] For the target shale gas well X1, the simulation results of the method and the conventional FVM simulation results without considering the heterogeneity are compared as Figures 4-6 It can be observed that the cumulative gas production obtained by the method is closer to the actual cumulative gas production, and the cumulative gas production obtained by the conventional FVM is lower than the actual cumulative gas production. From Figures 6-8 As shown in the figure, the difference between the data obtained by using the conventional FVM, the MSFVM and the field is not large, which proves that the MSFVM used in the embodiment, which maps the fine scale to the coarse scale, reduces the calculation amount, and the precision is not reduced, which proves that the method is feasible.

[0183] Those of ordinary skill in the art will realize and understand that all or some of the steps in the methods disclosed above and the functional modules / units in the systems and devices can be implemented as software, firmware, hardware, and appropriate combinations thereof. In hardware implementation, the division between the functional modules / units mentioned in the above description does not necessarily correspond to the division of physical components; for example, one physical component can have multiple functions, or one function or step can be performed by several physical components in cooperation. Some or all of the components can be implemented as software executed by a processor, such as a digital signal processor or a microprocessor, or as hardware, or as an integrated circuit, such as an application-specific integrated circuit. Such software can be distributed on computer-readable media, which can include computer storage media (or non-transitory media) and communication media (or transitory media). As is well known to those of ordinary skill in the art, the term computer storage media includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tapes, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by a computer. Furthermore, it is common and well understood by those of ordinary skill in the art that communication media typically embodies computer readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and can include any information delivery media.

Claims

1. A method of predicting gas-water two-phase deliverability of a fractured horizontal well, characterized in that, The method comprises: establishing a target block heterogeneous permeability Gaussian random field model and a target well reservoir geology model; establishing corresponding coarse-scale grids, fine-scale grids and dual control volume grids based on the target well reservoir geology model; determining dual basis functions for mutual mapping of coarse-scale grids and fine-scale grids according to the coarse-scale grids, fine-scale grids and dual control volume grids; discretizing control equations corresponding to the coarse-scale grids and calculating pressures in fractures and matrixes on the coarse grids according to the target block heterogeneous permeability Gaussian random field model; mapping the pressures in fractures and matrixes on the coarse grids to the fine-scale grids according to the dual basis functions; determining gas-water two-phase productivity of the target well according to the pressures in fractures and matrixes on the fine-scale grids.

2. The method of predicting gas-water two-phase deliverability of a fractured horizontal well of claim 1, wherein, The dual control volume is: where λ is the permeability, k is the number of grid nodes, is the dual control body, is the kth dual basis function.

3. The method of predicting gas-water two-phase deliverability of a fractured horizontal well of claim 1, wherein, The control equations comprise fracture grid control equations and matrix grid control equations.

4. The method for predicting gas-water two-phase productivity of a fractured horizontal well according to claim 3, characterized in that, the fracture grid control equation for single-phase flow is: In the above equation, K m is the permeability of the matrix, K F is the permeability of the fracture, μ is the viscosity, p is the pressure, p is the density, g is the gravitational acceleration, γ mf is the flow rate from the fracture to the matrix, γ Fm is the flow rate from the matrix to the fracture, q F is the flow rate in the fracture, and q m is the flow rate in the matrix.

5. The method of predicting gas-water two-phase deliverability of a fractured horizontal well of claim 4, wherein, the process for discretizing the fracture grid control equation for single-phase flow is: where A c ij represents the influence of j on i, ij = m, m is the matrix, F is the fracture; A c EDFM represents the embedded discrete fracture coarse grid matrix, p f EDFM represents the pressure of the embedded discrete fracture, q f EDFM represents the flow of the embedded discrete fracture, and EDFM is the embedded discrete fracture grid model.

6. The method of predicting gas-water two-phase deliverability of a fractured horizontal well of claim 4, wherein, the permeability in the matrix system is: In the above equation, K m is the permeability in the matrix system, a is the thinning coefficient, k m0 is the absolute permeability of the shale matrix, b is the slip coefficient, k n is the Knudsen number; the permeability in the fracture system is: In the above equation, K ζ is the permeability in the fracture system, ζ is the main fracture grid number, K ζ0 is the initial permeability, c ζ is the stress sensitivity coefficient, P ζ0 is the initial pressure.

7. The method for predicting gas-water two-phase productivity of a fractured horizontal well according to claim 1, characterized in that, the matrix-to-fracture crossflow rate is: gamma Fm = CI lambda m (p F -p m ); the fracture-to-matrix crossflow rate is: gamma mF = CI lambda m (p m -p F ); In the above equation, γ mf is the flow rate of the fracture to the matrix, γ Fm is the flow rate of the matrix to the fracture, CI is the connectivity coefficient of the fracture and the matrix, λ m is the permeability of the matrix, p F , p m .

8. The method for predicting gas-water two-phase productivity of a fractured horizontal well according to claim 7, characterized in that, the connectivity coefficient CI of the fracture and matrix is: wherein A is the fracture area of a complete fracture segment in the corresponding matrix grid after the complete fracture segment is divided by the matrix grid, and d is the average distance between the complete fracture segment and the center of the corresponding matrix grid.

9. An apparatus for predicting gas-water two-phase deliverability of a fractured horizontal well, comprising: The device comprises a memory and a processor; the memory is used to save a program for predicting gas-water two-phase productivity of a fractured horizontal well, and the processor is used to read and execute the program for predicting gas-water two-phase productivity of a fractured horizontal well, and execute the method according to any one of claims 1-8.

10. A computer readable storage medium, wherein a data processing program is stored on the computer readable storage medium, and the data processing program is executed by a processor to perform the method for predicting gas-water two-phase productivity of a fractured horizontal well according to any one of claims 1-8.

Images