Gas-water seepage field simulation method based on multi-scale analysis

Abstract

The invention discloses a multi-scale analysis-based gas-water seepage field simulation method, which specifically comprises the following steps of: based on a digital rock core reconstructed by a CT (Computed Tomography) scanning image of a typical rock sample of a target reservoir, adopting a microcosmic pore network flow simulation technology; the change rule of the stratum seepage capacity in different stratum pressure gradients and water drive processes, the use degree of a pore network and the water drive process in a pore are calculated, and then a correlation chart with the rock permeability and the pressure gradients is established and used for rapid measurement and calculation of the actual reserve use degree and the recovery degree of the water-invaded sandstone gas reservoir. The method starts from a rock microcosmic pore channel, the scale is amplified to a macroscopic oil reservoir during application, the method has the cross-scale characteristic, and the method is used for rapidly measuring, calculating and evaluating the reserve utilization degree and the recovery degree of the actual water-invaded sandstone gas reservoir and is beneficial to rapid development of the gas reservoir engineering scheme design in a mine field.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of numerical simulation and analysis, and particularly relates to a gas-water seepage field simulation method based on multi-scale analysis. BACKGROUND

[0002] Due to reservoir heterogeneity, the water invasion process of water invasion sandstone gas reservoirs is non-equilibrium, which is not only controlled by reservoir properties, but also significantly affected by gas-water distribution, formation pressure gradient, well pattern, well spacing and production pressure difference. Due to insufficient gas charging during accumulation, the reservoir pore contains a large amount of original formation water, and the formation seepage capacity is significantly affected by the gas-water distribution, and the reservoir seepage capacity is related to the pressure gradient and the gas-water distribution, thereby leading to the complexity of the producing degree and the producing process of the reserves in the actual development process, which is difficult to predict by static reservoir parameters or conventional theoretical models based on Darcy linear seepage. In order to improve the development effect of water invasion sandstone gas reservoirs, it is necessary to establish a corresponding mathematical model for the variation law of the pore network flow capacity depending on the seepage environment, analyze the producing degree and the producing process of the reserves of the gas reservoir, and predict the corresponding development index.

[0003] Chinese patent (application number: 201811223221.X, publication number: CN109284571A, publication date: 2019-01-29) discloses a carbon dioxide displacement shale gas multi-scale multi-field coupled seepage mathematical modeling method, which relates to multi-scale numerical simulation, and aims to optimize the development and gas injection process of shale gas reservoirs; Chinese patent (application number: 201810050202.5, publication number: CN108266185A, publication date: 2018-07-10) discloses a unconventional reservoir volume reconstruction multiple pore medium productivity contribution evaluation method, which is suitable for the field of oil and gas development, but it focuses on reservoir productivity evaluation and does not involve multi-scale analysis and simulation of gas-water seepage field. SUMMARY

[0004] The purpose of the present application is to provide a gas-water seepage field simulation method based on multi-scale analysis, which can accurately simulate the gas-water two-phase flow in the reservoir and improve the accuracy and reliability of the simulation.

[0005] The technical solution adopted by the present application is a gas-water seepage field simulation method based on multi-scale analysis, which is implemented according to the following steps:

[0006] Step 1, obtaining a typical rock sample; first, the main distribution interval of the target reservoir permeability is divided into 10-20 segments, and the average permeability of each segment is used as the permeability constraint condition for screening typical samples;

[0007] Step 2, CT scanning of the typical rock sample to obtain a two-dimensional CT scanning gray image;

[0008] Step 3, based on the two-dimensional CT scanning gray image, the rock pore space three-dimensional reconstruction is carried out;

[0009] Step 4, the percolation parameters of the digital core are calculated, including porosity and permeability;

[0010] Step 5, the gas-water two-phase flow is simulated for each rock sample, and the pore network volume sweep efficiency under different pressure gradients is calculated;

[0011] Step 6: regression is established for the pore sweep efficiency / reserve producing degree equation of different permeability and each pressure gradient;

[0012] Step 7: for the permeability distribution range of the target reservoir and the upper and lower limits of the formation pressure gradient during production, the reserve producing degree is calculated according to the equation obtained in step 6;

[0013] Step 8: the reserve producing degree that can be reached under the conditions of production pressure difference and well spacing is measured and calculated; for the target reservoir, the formation pressure gradient distribution is calculated according to the production pressure difference and well spacing, and the reserve producing degree at this time is measured and calculated according to the equation of step 6 according to the formation pressure gradient and the formation permeability of the target point;

[0014] Step 9: the flow rate of each level of pore size is calculated by using the digital core, and the relative time of each level of pore to complete water drive is further calculated, so as to obtain the corresponding relationship between the total water passing multiple and the reserve producing degree in the pore;

[0015] Step 10: the reserve producing process corresponding to different water drive pore volume multiples under the current sweep degree is calculated;

[0016] Step 11: a reserve producing process chart is drawn; for the permeability distribution range of the target reservoir and the maximum water invasion amount during production, the basic parameters of the chart are set, and then the chart curve is calculated according to the reserve producing process equation established in step 10;

[0017] Step 12: the producing degree that can be reached by each water invasion amount of the water invasion sandstone gas reservoir is predicted; for the water invasion area of the target reservoir, the current producing degree is measured and calculated according to the chart according to the ratio of the cumulative water invasion amount of the formation to the total pore volume of the formation and the average permeability of the formation, and the water passing pore volume multiple required to reach each level of producing degree is predicted, so as to indirectly measure and calculate the required production time.

[0018] In step 3, specifically, the median filter method is used to eliminate the noise points in the two-dimensional CT scanning gray image, then the initial threshold value is set, the two-dimensional CT scanning gray image is converted into a binary image by using the image segmentation technology, and the binary image is smoothed to remove the isolated rock skeleton, so as to realize the three-dimensional reconstruction of the digital core.

[0019] The characteristics of the present application also include,

[0020] In step 4, specifically:

[0021] Step 4.1, the pore space of the digital core is characterized, and the maximum sphere method is used to calculate the radius and spatial position coordinates of all spheres in the rock sample.

[0022] The pore space of the rock can be equivalently characterized as spherical pores, single pipe bundles, and complex pipe networks. Nodes are used to simulate pores, and pipes are used to simulate the connectivity between pores, and the porous medium of the rock is characterized as a three-dimensional network of interconnected spaces composed of nodes and pipes. This three-dimensional network has pore size, pore position, and pore size.

[0023] A micro-radius sphere is placed at any point in the pore, and then the sphere radius is increased until the sphere surface contacts the rock wall, thereby obtaining the largest sphere that can be placed in the pore space centered at the pore point. The remaining space is also processed, and the pore space is characterized as a series of adjacent spheres. The radius and spatial position coordinates of all spheres in the rock sample are calculated, as shown in equation (1):

[0024]

[0025] In the formula, N is the number of spheres identified by scanning; r is the radius of sphere i; x, y, and z are the spatial coordinates of sphere i.

[0026] Step 4.2, calculate the porosity of the digital core; specifically:

[0027] In accordance with the sampling direction and the main flow direction in the formation, the entire space of the digital core is discretized into rectangular grids according to the set grid size. The position and porosity of each grid containing pores are calculated to obtain the porosity of the entire grid.

[0028] Step 4.3, calculate the connectivity between each grid.

[0029] 1) Determine the spatial physical connectivity of the three-dimensional grid: if the grid porosity between grid i-1 and grid i+1 is greater than 0, then grid i-1 and grid i+1 are physically connected; otherwise, they are not physically connected.

[0030] 2) Calculate the equivalent pore size of the physically connected pore: according to the volume equivalence principle, the pore volume of each connected grid is converted into the radius of a cylinder, which is the equivalent pore size R i-1 and R i+1 ;

[0031] 3) Calculate the flow resistance of the connected pore: according to the Young-Laplace equation of the pore, the flow resistance of the pore is calculated, as shown in equation (2):

[0032]

[0033] P = 2σR1R2 cgw is the flow resistance in the micro-pore channel when gas-water two-phase flow exists, unit Pa; σ gwr is the gas-water interfacial tension, unit N / m; R1, R2 are the curvature radii in the pore channel, unit m;

[0034] 4) Determine the flow mode: for the pore network with complex parallel connection relationship, according to the inlet and outlet diameters of the connected pore and the fluid properties, there are three flow modes in general:

[0035] Mode one: the driving force is greater than the maximum capillary force possessed by the minimum pore diameter, the pore system is fully affected, all the pores participate in the flow, and the effective permeability of the rock sample is equal to the absolute permeability of the rock sample;

[0036] Mode two: the driving force is greater than the minimum capillary force possessed by the maximum pore diameter and less than the maximum capillary force possessed by the minimum pore diameter, at this time, the large pores in parallel participate in the flow, the small pores do not participate in the flow, and part of the pores in the pore network participate in the flow; for the entire pore system, the apparent flow area is reduced, and the effective permeability is less than the absolute permeability;

[0037] Mode three: the driving force is less than the minimum capillary force possessed by the maximum pore diameter, and all the connected pores do not participate in the flow, at this time, the pore network system is not affected, and the effective permeability is 0;

[0038] Step 4.4, calculate the apparent permeability of the connected and flowing pores;

[0039] 1) Dimension reduction: divide the three-dimensional grid into a series of two-dimensional cross-section grid; set the apparent flow direction as X, then the two-dimensional direction of the cross-section layer is YZ; let the total number of the “effective” grid of the i-th layer be N i ; set the inlet pressure as P1 and the outlet pressure as P2, then the driving pressure difference is: ΔP = P1-P2; let the initial value of the number of the effective grid of the i-th layer be M i , and M i =N i ;

[0040] 2) According to the resistance parallel principle, calculate the total flow resistance of each cross-section layer;

[0041] Let be the maximum flow resistance from the i, j, k grid to the i+1, j, k grid along the X direction;

[0042] According to the resistance parallel principle, the total flow resistance from the i-th layer to the i+1-th layer is as shown in formula (5):

[0043]

[0044] 3) In the flow direction, the total flow resistance is calculated cumulatively according to the principle of resistance series, as shown in equation (6):

[0045]

[0046] The partial driving pressure on each cross-sectional layer is obtained, as shown in equation (7):

[0047]

[0048] In a cross-sectional layer, the driving force is compared with the flow resistance of each effective grid to determine whether the grid flows, thereby obtaining the effective grids of each cross-sectional layer in the X direction, as shown in equation (8), and the total number is recorded as M i

[0049]

[0050] If Mi'≠Mi, repeat 2)~3) until M i i .

[0051] At this point, the driving force ΔP is obtained, and the geometric configuration of the flow path in the pore network is obtained;

[0052] The connected and flowing grids are divided into a composite circuit composed of a series-parallel resistance, the total resistance is calculated, and the apparent permeability of the digital core is calculated;

[0053] According to the Darcy linear flow formula, the core permeability is calculated from the core flow experiment data, as shown in equation (9):

[0054]

[0055] In the formula: K is the permeability, with the unit of D; Q is the test flow rate, with the unit of cm 3 / s; A is the core cross-sectional area, with the unit of cm 2 ; μ is the fluid viscosity, with the unit of mPa·s; Δp is the test pressure difference, with the unit of atm; ΔL is the core length, with the unit of cm;

[0056] At the microscale, the internal flow space of a porous medium is composed of a number of interwoven pore networks, a single pore can be simplified as a thin tube, and in the laminar flow state, the flow in the thin tube satisfies the Poiseulle law, as shown in equation (10):

[0057]

[0058] In the formula: ΔP is the pressure difference between the two ends of the thin tube, with the unit of Pa; r is the radius of the thin tube, with the unit of m; q is the flow rate of the thin tube, with the unit of m 3 / s.​​

[0059] For the tube bundle with multiple tubes, assuming the number of tubes per unit cross-sectional area is n, the total flow of all tubes on the cross-sectional area A is shown in equation (11):

[0060]

[0061] In the formula: A is the cross-sectional area, unit m 2 ; n is the surface density of the tube, unit m -2 ;

[0062] Simplify the three-dimensional grid in the flow direction to a series of cross-sectional layers in series, and the effective grid in the cross-sectional layer is in parallel mode, which meets the flow mode of multiple tubes;

[0063] Let the side length of the grid be a, for the i-th cross-sectional layer, the driving pressure is Δp i , the number of effective grids is n i , and the equivalent aperture of the m-th effective grid is r m , then the flow of the m-th grid of the cross-sectional layer is shown in equation (12):

[0064]

[0065] For the uniform cubic effective grid, the equivalent channel length is the grid side length a, the actual channel in the porous medium is curved, and the ratio of the actual length to the straight-line distance between the two ends of the channel is the tortuosity τ. The tortuosity of each grid is not the same, so the flow formula of the grid is shown in equation (13):

[0066]

[0067] The total flow of the cross-sectional layer is shown in equation (14):

[0068]

[0069] The minimum flow of each cross-sectional layer in series determines the total flow of the grid system. Let Nx be the total number of cross-sectional layers in series in the flow direction, and Q be the total flow of the grid system. Then equation (15) is obtained:

[0070]

[0071] The core cross-sectional area in the X flow direction is N y ·N z , the pipe flow equation is substituted into the Darcy formula, and the apparent permeability of the unit body can be expressed as equation (16):

[0072]

[0073] The complete expression of the above formula is equation (17):

[0074]

[0075] For the grid system of standard cube, N x = N y = N z = N, formula (18) is obtained:

[0076]

[0077] Let the total volume of the pore network be V φ Then the apparent porosity of the effective grid system is formula (19):

[0078]

[0079] In step 6, the reserve producing degree equation is shown as formula (20):

[0080] E = (A1K 3 +A2K 2 +A3K+A4)·LnG+(B1K 2 +B2K+B3)(20)

[0081] In the formula, E is the sweep efficiency of the pore network; G is the driving pressure gradient; A1, A2, A3, A4, B1, B2 and B3 are equation coefficients.

[0082] The beneficial effects of the present application are that the method is reasonable and accurate based on the micro-pore network flow theory and the actual reservoir rock seepage parameters; the method is scientific and advanced by using the latest pore network nonlinear flow simulation technology; the method has wide application with a wide range of permeability and pressure gradient; the method is fast and efficient based on the large data analysis of the theoretical calculation results; the method has a cross-scale feature by starting from the rock micro-pore and scaling up to the macroscopic reservoir in application. The whole method is used for the rapid calculation and evaluation of the reserve producing degree and recovery degree of the actual water invasion sandstone gas reservoir, which is beneficial to the rapid development of the engineering scheme design of the gas reservoir in the field. BRIEF DESCRIPTION OF DRAWINGS

[0083] Figure 1 The flow chart of the gas-water seepage field simulation method of the present application based on multi-scale analysis;

[0084] Figure 2 It is a typical rock sample selection graph according to the permeability distribution in the embodiment of the present application;

[0085] Figure 3 It is a rock sample by rock sample simulation sweep efficiency curve graph in the embodiment of the present application;

[0086] Figure 4 is a graph of the relationship between the regression coefficient A and the permeability in the embodiment of the present application;

[0087] Figure 5 is a graph of the relationship between the regression coefficient B and the permeability in the embodiment of the present application;

[0088] Figure 6 is a graph of the reserve producing degree calculated by the regression equation in the embodiment of the present application;

[0089] Figure 7 is a graph of the water drive producing process curve of each rock sample in the embodiment of the present application;

[0090] Figure 8 is a graph of the relationship between the regression coefficient A and the permeability in the embodiment of the present application;

[0091] Figure 9 is a graph of the relationship between the regression index B and the permeability in the embodiment of the present application;

[0092] Figure 10 is a graph of the producing degree calculated by the regression equation in the embodiment of the present application. DETAILED DESCRIPTION

[0093] The present application will be described in detail below in combination with the drawings and specific embodiments.

[0094] Embodiment 1

[0095] The present application is based on a gas-water seepage field simulation method of multi-scale analysis, through a digital core reconstructed based on a CT scanning image of a typical rock sample of a target reservoir, a micro-pore network flow simulation technology is adopted to calculate the variation law of the seepage capacity of the formation in different formation pressure gradients and water drive processes and the producing degree of the pore network and the water drive process in the pore, and then a correlation graph plate of the rock permeability and the pressure gradient is established, which is used for the rapid calculation of the reserve producing degree and the producing degree of the actual water-invasion sandstone gas reservoir, and is helpful for the rapid development of the engineering scheme design of the water-invasion sandstone gas reservoir in the mine field.

[0096] Embodiment 2

[0097] The present application is based on a gas-water seepage field simulation method of multi-scale analysis, as shown in Figure 1 the following steps are implemented:

[0098] Step 1, a typical rock sample is obtained;

[0099] Considering that there can be obvious differences in the rock properties inside the reservoir, representative rock samples need to be selected from different types of strata and rocks. To ensure the coverage of sample properties to the target reservoir and improve sampling efficiency, the main distribution range of the permeability of the target reservoir is first determined, and the main distribution range of the permeability is divided into 10-20 segments, and the average permeability of each segment is used as the permeability constraint condition for selecting typical samples.

[0100] Since the accuracy of the method depends on the acquisition of typical rock samples of the reservoir, it is necessary to have a full understanding of the mineral composition, clay composition, and pore microstructure of the local reservoir rock, and to avoid using specific rock samples as typical rock samples as much as possible;

[0101] Step 2, CT scanning of the typical rock sample to obtain a two-dimensional CT scanning grayscale image;

[0102] According to the working parameters of the scanning equipment and the main seepage space range, seepage environment, and rock mineral and pore characteristics in the reservoir, the size and accuracy of the scanning range are determined, and then scanning is performed to obtain a series of two-dimensional CT scanning grayscale images of the typical rock sample;

[0103] Step 3, three-dimensional reconstruction of the rock pore space based on the two-dimensional CT scanning grayscale image;

[0104] Specifically, the median filter method is used to eliminate noise points in the two-dimensional CT scanning grayscale image, and then an initial threshold is set, the image segmentation technique is used to convert the two-dimensional CT scanning grayscale image into a binary image, and the binary image is smoothed to remove isolated rock skeletons, and three-dimensional reconstruction of the digital core is realized;

[0105] Step 4, calculation of the seepage parameters (including porosity and permeability) of the digital core; specifically,

[0106] Step 4.1, characterization of the pore space of the digital core, and calculation of the radius and spatial position coordinates of all spheres in the rock sample using the maximum sphere method equivalent pore;

[0107] The pore space of the rock can be equivalently characterized as a spherical pore, a single pipe bundle, and a complex pipe network. Nodes are used to simulate pores, and pipes are used to simulate the connection between pores, and the rock porous medium is characterized as a three-dimensional network structure of interconnected nodes and pipes; the three-dimensional network structure has pore size, pore position, and pore size.

[0108] A micro-radius sphere is placed at any point in the pore, and then the sphere radius is increased until the sphere surface contacts the rock wall, thus obtaining the largest sphere that can be placed in the pore space with the pore point as the center. The remaining space is also processed in the same way, and the pore space is characterized as a series of adjacent size spheres. The radius and spatial coordinates of all spheres in the rock sample are calculated, as shown in equation (1):

[0109]

[0110] In the formula, N is the number of spheres identified by scanning; r is the radius of sphere i; x, y, and z are the spatial coordinates of sphere i;

[0111] Step 4.2, calculate the porosity of the digital core; specifically:

[0112] In accordance with the sampling direction and the main flow direction in the formation, the entire space of the digital core is discretized into rectangular grids according to the set grid size. The position and porosity of each grid containing pores are calculated to obtain the porosity of the entire grid.

[0113] Step 4.3, calculate the connectivity between each grid;

[0114] 1) Judge the spatial physical connectivity of the three-dimensional grid: if the grid porosity between grid i-1 and grid i+1 is greater than 0, then grid i-1 and grid i+1 are physically connected; otherwise, they are physically disconnected;

[0115] 2) Calculate the equivalent pore diameter of the physically connected channel: according to the volume equivalence principle, the pore volume of each connected grid is converted into the radius of a cylinder, which is the equivalent pore diameter R i-1 and R i+1 of the connected channel;

[0116] 3) Calculate the flow resistance of the connected channel: calculate the flow resistance of the channel according to the Young-Laplace equation, as shown in equation (2);

[0117]

[0118] In the formula, P cgw is the flow resistance in the micro-pore channel when gas and water are present, with a unit of Pa; σ gwr is the gas-water interfacial tension, with a unit of N / m; R1 and R2 are the radii of curvature in the channel, with a unit of m;

[0119] If the driving force between the connected grids is greater than the flow resistance, flow will occur; otherwise, even if the grids are physically connected, there will be no flow between them because the driving force does not exceed the flow resistance.

[0120] The flow resistance calculation of the connected channel is divided into the following two cases:

[0121] Case 1: Equal diameter channels, the equation is simplified as equation (3):

[0122]

[0123] Where: θ gwr is the receding contact angle of gas-water interface; r is the channel radius, unit is m;

[0124] Case 2: Unequal diameter channels, the equation is equation (4):

[0125]

[0126] Where: r1 is the channel radius at the inlet end, unit is m; r2 is the channel radius at the outlet end, unit is m;

[0127] 4) Determine the flow pattern: for the channel network with complex parallel connection, according to the inlet and outlet diameters of the connected channels and the fluid properties, there are three flow patterns in general:

[0128] Pattern 1: The driving force is greater than the maximum capillary force of the smallest diameter in the channel network, the channel system is fully affected, all channels participate in the flow, and the effective permeability of the rock sample is equal to the absolute permeability of the rock sample. At this time, the permeability only depends on the channel and is independent of the percolation environment;

[0129] Pattern 2: The driving force is greater than the minimum capillary force of the maximum diameter and less than the maximum capillary force of the minimum diameter, at this time, the large channels in parallel participate in the flow, and the small channels do not participate in the flow, part of the channels in the channel network participate in the flow. For the entire channel system, the apparent flow area is reduced, and the effective permeability is less than the absolute permeability;

[0130] Pattern 3: The driving force is less than the minimum capillary force of the maximum diameter, and all connected channels do not participate in the flow, at this time, the channel network system is not affected, and the effective permeability is 0.

[0131] Step 4.4, calculate the real-time (apparent) permeability of the connected and flowing channels;

[0132] According to the channel flow pattern, only when the driving force is greater than the flow resistance between channels can the channels maintain flow. When the driving force changes, the number and location of effective grids participating in the flow will change accordingly. Therefore, within a certain driving force range, the flow path in the micro channel network is changing.

[0133] The following steps are used to simplify the calculation process of dynamic flow path:

[0134] 1) Dimension reduction: the three-dimensional grid is divided into a series of two-dimensional section grids; let the apparent flow direction be X, then the two-dimensional direction of the section layer is YZ; let the total number of the "effective" grid (equivalent aperture > 0) of the i-th layer be N i ;

[0135] Let the inlet pressure be P1 and the outlet pressure be P2, then the driving pressure difference is: ΔP = P1 - P2

[0136] Let the initial value of the number of effective grids of the i-th layer be M i , then: M i = N i ;

[0137] 2) According to the parallel resistance principle, the total flow resistance of each section layer is calculated;

[0138] Let be the maximum flow resistance of the grid located at i, j, k along the X direction to the grid i+1, j, k, ignoring the viscous resistance, and assuming that the oil and water are in a flow disadvantageous position;

[0139] According to the parallel resistance principle, the total flow resistance of the i-th layer to the i+1-th layer is shown in formula (5):

[0140]

[0141] 3) In the flow direction, the total flow resistance is calculated according to the series resistance principle, as shown in formula (6):

[0142]

[0143] The partial driving pressure on each section layer is obtained, as shown in formula (7):

[0144]

[0145] In a section layer, the driving force is compared with the flow resistance of each effective grid to determine whether the grid is flowing, thereby obtaining the effective grid of each section layer in the X direction, as shown in formula (8), and the total number is M i :

[0146]

[0147] If M i ≠ M i ', repeat 2) ~ 3) until M

[0148] At this point, the driving force ΔP is obtained, and the geometric configuration of the flow path in the pore network is obtained

[0149] Under the current driving pressure condition, all the grids that are spatially connected (a necessary condition for flow) and have driving force greater than flow resistance (a sufficient condition for flow) constitute the spatial flow path in the porous medium. According to the hydroelectric similarity principle, the connected and flowing grids are divided into a composite circuit composed of a number of series-parallel resistances, the total resistance is calculated, and then the apparent permeability of the digital core is calculated;

[0150] According to the fixed step length, the driving pressure is gradually increased, more and more grids with smaller equivalent pore diameter will participate in the flow, and the apparent permeability of the digital core will also gradually increase; when all the effective grids participate in the flow, the corresponding apparent permeability will tend to be stable.

[0151] According to the calculation principle of dynamic flow path and real-time effective permeability, from the inlet end surface to the outlet section, as the driving pressure continuously increases, when the first flow path is fully connected, the digital core begins to have seepage capacity, and the permeability begins to be greater than 0. This path will have larger pore diameter, better connectivity and lowest full-line seepage resistance. The driving pressure corresponding to this seepage resistance is the starting pressure of the digital core grid system.

[0152] According to the above calculation principle, the starting pressure gradient is a function of the spatial configuration of the pore network (including pore diameter, pore length), fluid viscosity, rock wettability, etc.

[0153] As the driving pressure increases, more and more pores will participate in the flow, and the permeability will gradually increase; when all the effective pores fully participate in the flow, the permeability will remain constant, and the apparent flow rate will be linearly related to the driving pressure, and the seepage characteristics will enter the linear stage from the nonlinear stage.

[0154] According to the Darcy linear flow formula, the core permeability is calculated from the core flow experiment data, as shown in formula (9):

[0155]

[0156] In the formula: K is the permeability, unit is D; Q is the test flow rate, unit is cm 3 / s; A is the core cross-sectional area, unit is cm 2 ; μ is the fluid viscosity, unit is mPa·s; Δp is the test pressure difference, unit is atm; ΔL is the core length, unit is cm;

[0157] At the microscale, the internal flow space of the porous medium is composed of a number of interwoven pore networks, a single pore can be simplified as a thin tube, and in the laminar flow state, the flow in the thin tube satisfies the Poiseulle law, as shown in formula (10):

[0158]

[0159] where ΔP is the pressure difference between the two ends of the capillary, r is the capillary radius, and q is the capillary flow rate. 3 / s.

[0160] For a tube bundle with multiple capillaries, assuming that the number of capillaries per unit cross-sectional area is n, the total flow rate of all tubes on the cross-sectional area A is shown in equation (11):

[0161]

[0162] where A is the cross-sectional area, n is the surface density of the capillaries, and r is the capillary radius. 2 ; n is the surface density of the capillaries, and r is the capillary radius. -2 ;

[0163] For the grid system of digital cores, the pores within each grid are equivalent to a capillary in volume, which has three flow directions along the coordinate axes in geometry, and the flow law can be described by Poiseulle's law. For the pore network, based on the principle of hydroelectric similarity, the three-dimensional grid is simplified into a series of cross-sectional layers in the flow direction, and the effective grid within the cross-sectional layer is in parallel mode, which meets the flow mode of multiple capillaries.

[0164] Let the side length of the grid be a, and for the i-th cross-sectional layer, the driving pressure is Δp i , the number of effective grids is n i , and the equivalent aperture of the m-th effective grid is r m . The flow rate of the m-th grid in the cross-sectional layer is shown in equation (12):

[0165]

[0166] For a uniform cubic effective grid, the equivalent channel length is the grid side length a, and the actual channel in the real porous medium is curved, with the ratio of the actual length to the straight-line distance between the two ends of the channel being the tortuosity τ. The tortuosity of each grid is not the same, so the flow rate formula of the grid is shown in equation (13):

[0167]

[0168] The total flow rate of the cross-sectional layer is shown in equation (14):

[0169]

[0170] The minimum flow rate of each cross-sectional layer in series determines the total flow rate of the grid system. Let Nx be the total number of cross-sectional layers in series in the flow direction, and Q be the total flow rate of the grid system. Equation (15) is obtained:

[0171]

[0172] X Core cross-sectional area in flow direction N y ·N z Substitute the pipe flow equation into the Darcy formula, and the apparent permeability of the unit body can be expressed as formula (16):

[0173]

[0174] The complete expression of the above formula is formula (17):

[0175]

[0176] For the grid system of a standard cube, N x = N y = N z = N, formula (18) is obtained:

[0177]

[0178] Let the total volume of the pore network be V φ The apparent porosity of the effective grid system is formula (19):

[0179]

[0180] The tortuosity τ in the above formula can be an empirical value of 1.4, or obtained by fitting according to the core gas measured permeability;

[0181] Step 5: Simulate gas-water two-phase flow for each rock sample, and calculate the volume sweep efficiency of the pore network under different pressure gradients;

[0182] Step 6: Regression to establish the pore sweep efficiency / reserve mobilization degree equation of different permeabilities and various pressure gradients, as shown in formula (20):

[0183] E = (A1K 3 +A2K 2 +A3K+A4)·LnG+(B1K 2 +B2K+B3) (20)

[0184] In the formula: E is the sweep efficiency (reserve mobilization degree) of the pore network; G is the driving pressure gradient; A1, A2, A3, A4, B1, B2, B3 are equation coefficients.

[0185] Step 7: According to the permeability distribution range of the target reservoir and the upper and lower limits of the formation pressure gradient during production, calculate the reserve mobilization degree according to the equation obtained by regression in step 6;

[0186] Step 8: Calculate the reserves producing degree under the conditions of production pressure difference and well spacing; for the target reservoir, calculate the formation pressure gradient distribution according to the production pressure difference and well spacing, and according to the formation pressure gradient and formation permeability at the target point, calculate the reserves producing degree at this time by the equation in step 6;

[0187] Step 9: For typical rock samples with different permeabilities, simulate the water drive process in the pore channel based on the pore size combination mode; the reservoir pore channel space has micro-heterogeneity, and there are pore channel network systems with various pore sizes and various connectivity, but for the same sedimentary reservoirs, it can be generally considered that the apparent permeability can represent the combination mode of the pore channel network. Using digital core, calculate the flow rate of each pore size, and further calculate the relative time for each pore channel to complete water drive, and obtain the corresponding relationship between the total water passing multiple and the reserves producing degree in the pore channel.

[0188] Step 10: Calculate the reserves producing process corresponding to different water drive pore volume multiples under the current sweep degree;

[0189] Draw the pore volume water passing multiple and producing degree relationship curve, and regress according to the power function equation, as shown in equation (21):

[0190] R=A·PV B (21);

[0191] In the formula: R is the producing degree; PV is the pore volume;

[0192] The calculation formula of equation coefficients A and B is shown in equations (22) and (23):

[0193] A=F3(K)=C·K D (22);

[0194] B=F4(K)=G·K H (23);

[0195] A and B are equation coefficients, both of which are functions of absolute permeability when the formation pore channel network is fully swept, and the functions are obtained by power function equation regression;

[0196] Step 11: Draw the reserves producing process chart;

[0197] For the permeability distribution interval and the maximum water influx during the production of the target reservoir, set the basic parameters of the chart, and then calculate the chart curve according to the reserves producing process equation established in step 10 regression;

[0198] Step 12: Predict the producing degree that can be reached by each water influx of water-invasion sandstone gas reservoirs;

[0199] For the water invasion area of the target reservoir, the ratio of the cumulative water invasion amount of the formation to the total pore volume of the formation and the average permeability of the formation are used to measure the current recovery degree according to the chart, and the water passing pore volume multiple required to reach each level of recovery degree is predicted, and the required production time is indirectly measured.

[0200] The method of the application uses digital cores based on real reservoir rock samples, calculates real-time pore network microscopic sweep efficiency, pore water drive process and effective permeability response by setting a pressure gradient, ensures that the set of calculation methods covers the actual complex percolation environment, saves the time, equipment, material and labor cost required by the conventional physical model experiment, reduces the system error of the physical model experiment method because of the strict microscopic percolation theory, improves the accuracy of the data, and ensures the universality of the whole set of calculation methods;

[0201] The method of the application applies the microscopic calculation results to the macroscopic reservoir, obtains the sweep efficiency of each spatial position in the reservoir under the percolation environment and the reserve producing process at each water drive time point by calculating the pressure profile and pressure gradient profile in the reservoir, so as to realize the amplification of the microscopic pore network calculation results to the reservoir scale, which can be used to measure the reserve producing degree, future recovery degree and well pattern densification and production pressure difference potential of the actual reservoir under the current well pattern, well spacing and production pressure difference conditions, and embodies the guiding significance of the method to the actual production.

[0202] Example 3

[0203] The gas-water percolation field simulation method based on multi-scale analysis specifically comprises:

[0204] Step 101: Collect, analyze and obtain the main distribution range of the permeability of the target reservoir. According to the sample test of the target gas reservoir 203, the permeability distribution range is 5mD-99.5mD, and the average is 27mD

[0205] Step 102: According to the main distribution range of the permeability, the representative typical rock samples are selected in sections. According to the main distribution range of the permeability, 45 typical rock samples are obtained from the original core, and the core meets the CT scanning test requirements;

[0206] Step 103: CT scanning is performed on the typical rock samples. The scanning range is set to a cube with a side length of 6mm-11mm, a double-source double-detector CT scanning machine is selected, the scanning machine has a 450KV ray source, a 300KV micro-focus ray source, a 400cm*400cm large-size flat panel detector and an extended linear array detector, and the scanning accuracy is 1um;

[0207] Step 104: Based on CT scan images, the rock pore space is reconstructed in three dimensions. A 1000x1000x1000 voxel grid system is established, and a 0-1 digital matrix is input, with 0 indicating that the voxel grid node is a pore and 1 indicating that the rock is filled. The modeling of 45 rock samples shows that the pore volume changes with the pore diameter in a general skewed normal distribution: the number of small pores is large, but the pore diameter is small, so the total volume of small pore diameter pores is small; the pore diameter of large pores is large, but the number is small, so the total volume of large pore diameter pores is also small; the pore diameter of medium pores is medium, and the number is medium, so the total volume is the largest. However, some rock samples contain a certain number of dissolution pores, so there will be a second volume peak at a larger pore diameter.

[0208] Step 105: Calculate the full wave and permeability of the digital core pore. The calculation results show that the full wave and permeability of the 45 rock sample digital core simulation calculation is in the range of 8.60mD-58.85mD, with an average of 27.10mD, and the measured air permeability of the core is in the range of 9.70mD-59.60mD, with an average of 28.22mD, and the relative fitting accuracy is 77.87%-98.74%, with an average of 90.02%. According to the relationship curve between the calculated starting pressure gradient and permeability, according to the permeability segmentation, as shown in Figure 2 , 11 typical cores with strong correlation are selected as the basis for further data calculation and analysis,

[0209] Step 106: Simulate gas-water two-phase flow for each rock sample, calculate the pore network volume sweep efficiency under different pressure gradients, and use semi-logarithmic equation regression for the calculation results of the 11 typical rock samples, as shown in Figure 3

[0210] Step 107: Regression to establish pore sweep coefficient / reserve producing degree equation for different permeability and pressure gradient. The sweep coefficient-pressure gradient curves of the 11 typical rock samples are different, with differences in equation coefficients. Further regression of the equation coefficients with respect to the formation permeability, as shown in Figure 4 and Figure 5 , the final equation for the pore network sweep coefficient (reserve producing degree) of the rock under different permeability and different pressure gradient is:

[0211] Ev=(0.0006K 3 -0.06K 2 +1.81K+27.34)Ln(G)+(0.018K 2 -0.26K-7.39)

[0212] In the formula: E is the sweep coefficient of the pore network (reserve producing degree); G is the driving pressure gradient; K is the formation permeability

[0213] ​Step 108: Draw the reserve producing degree chart. According to the above equation, the maximum pressure gradient of the target gas reservoir is 25 MPa / km, and the main formation permeability is set to be in the range of 5 mD to 60 mD. The chart is shown in FIG. 8. Figure 6

[0214] Step 109: Calculate the reserve producing degree under the conditions of production pressure difference and well spacing. For the target gas reservoir, the interwell formation pressure gradient is calculated, and referring to the chart, the current reserve producing degree under different reservoir properties can be obtained. The influence of well network densification or dilution and production pressure difference adjustment on the reservoir producing degree can also be calculated.

[0215] Step 110: For different permeability typical rock samples, simulate the water drive process in the pore based on the pore size combination mode. For 11 typical rock samples, it is considered that 11 permeabilities represent the corresponding pore combination types. The micro-pore network flow simulation method is used to calculate the reserve producing degree in the pore under different water drive pore volume multiples.

[0216] Step 111: Calculate the reserve producing process corresponding to different water drive pore volume multiples under the current sweep degree. Draw the relationship curve between pore volume water multiple and producing degree, as shown in FIG. 9. According to the power function regression equation, the equation coefficient is a function of the absolute permeability when the pore network is fully swept. The function is obtained by power function equation regression. Figure 7

[0217] Rg= 2.7159K -0.686 × PV 1.7949K^(-0.422)

[0218] In the formula, R is the producing degree, PV is the pore volume, and K is the formation permeability.

[0219] Step 112: Draw the reserve producing process chart. According to the above equation, the maximum water invasion water multiple of the target gas reservoir is 1000 PV, and the main formation permeability is set to be in the range of 5 mD to 60 mD. The chart is shown in FIG. 10. Figure 10

[0220] Step 113: Predict the producing degree that can be reached by each water invasion amount in the water invasion sandstone gas reservoir. For the target gas reservoir, different water invasion amounts are converted into multiples of underground pore volume, and referring to the chart, the producing degree corresponding to the current water multiple under different reservoir properties can be obtained. The growth characteristics of future producing degree with the increase of water invasion multiple can also be calculated, and the cumulative gas production can be indirectly calculated.

[0221] ​​​The main development technical parameters such as the reserve producing degree and the recovery degree measured by the method are greater than 80% in accordance with the reservoir coring analysis, pilot test, development adjustment scheme implementation data and gas reservoir numerical simulation calculation results of multiple water invasion gas sand bodies of Sebei gas field in Qinghai oil field, which meets the theoretical support needs of the field rapid design of the development implementation scheme of the gas reservoir, and has certain guiding significance for predicting, evaluating and screening the water power adjustment potential tapping scheme.

Claims

1. A method for simulating air-water seepage fields based on multi-scale analysis, characterized in that, Specifically, based on digital cores reconstructed from CT scan images of typical rock samples from the target reservoir, the flow simulation technology of micro-pore network is used to calculate the variation law of formation seepage capacity under different formation pressure gradients and water drive processes, as well as the degree of pore network activation and water drive process within the pores. Then, a correlation chart with rock permeability and pressure gradient is established for rapid calculation of the degree of activation and production of actual water-intruded sandstone gas reservoir reserves.

2. The air-water seepage field simulation method based on multi-scale analysis as described in claim 1, characterized in that, The specific steps are as follows: Step 1: Obtain typical rock samples; First, divide the main distribution range of the target reservoir permeability into 10 to 20 segments, and use the average permeability of each segment as the permeability constraint for screening typical samples. Step 2: Perform CT scans on typical rock samples to obtain two-dimensional CT scan grayscale images; Step 3: Reconstruct the three-dimensional space of rock pores based on the two-dimensional CT scan grayscale image; Step 4: Calculate the seepage parameters of the digital core, including porosity and permeability; Step 5: Simulate gas-water two-phase flow on a rock-by-rock sample and calculate the volumetric sweep efficiency of the pore network under different pressure gradients. Step 6: Regression to establish the pore sweep efficiency / reserve utilization equation for different permeabilities and pressure gradients; Step 7: Calculate the degree of reserve utilization based on the permeability distribution range of the target reservoir and the upper and lower limits of the formation pressure gradient during mining, according to the equation obtained from the regression in Step 6. Step 8: Calculate the achievable level of reserve utilization under the conditions of mineral production pressure differential and well spacing; For the target reservoir, the formation pressure gradient distribution is calculated based on the production pressure differential and well spacing. Based on the formation pressure gradient and formation permeability at the location of the target point, the degree of reserve utilization at this time is calculated using the equation in step 6. Step 9: Calculate the flow velocity of each orifice using digital core analysis, and further calculate the relative time for each orifice to complete water drive, to obtain the correspondence between the overall water flow ratio and the degree of reservoir recovery in the orifice; Step 10: Calculate the reserve mobilization process corresponding to different water-drive pore volume multiples under the current sweep level; Step 11: Draw a reserve utilization process chart; set the basic parameters of the chart based on the permeability distribution range of the target reservoir and the maximum water intrusion during the exploitation period, and then calculate the chart curve according to the reserve utilization process equation established in Step 10; Step 12: Predict the recovery level achievable by various water intrusion volumes in water-intruded sandstone gas reservoirs; for the water-intruded area of ​​the target reservoir, calculate the current recovery level based on the ratio of cumulative formation water intrusion to total formation pore volume and formation average permeability, and predict the required multiple of water-passing pore volume to reach each level of recovery, thereby indirectly calculating the required extraction time.

3. The air-water seepage field simulation method based on multi-scale analysis as described in claim 2, characterized in that, In step 3, the specific steps are as follows: median filtering is used to eliminate noise in the two-dimensional CT scan grayscale image, then an initial threshold is set, and image segmentation technology is used to convert the two-dimensional CT scan grayscale image into a binary image. The binary image is then smoothed to remove isolated rock skeletons and realize the three-dimensional reconstruction of the digital core.

4. The gas-water seepage field simulation method based on multi-scale analysis as described in claim 2, characterized in that, Step 4 specifically involves: Step 4.1: Characterize the pore space of the digital core and use the maximum sphere method to calculate the radius and spatial coordinates of all spheres in the rock sample to represent the pore size. The pore space of rocks can be equivalently represented as spherical pores, single tube bundles, and complex pipe networks. Nodes are used to simulate pores, and pipes are used to simulate the connection between pores. The porous rock medium is characterized as a three-dimensional spatial network structure composed of interconnected nodes and pipes. This three-dimensional spatial network structure has pore size, pore location, and channel size; A small-radius sphere is placed at any point in the pore as the center. The radius of the sphere is then increased until the surface of the sphere contacts the rock wall. This yields the largest sphere that can be placed in the pore space with the pore point as the center. The remaining space is treated in the same way, and the pore space is represented as a series of adjacent spheres of different sizes. The radius and spatial coordinates of all spheres in the rock sample are calculated, as shown in Equation (1): In the formula, N is the number of spheres identified by scanning; r is the radius of sphere i; and x, y, z are the spatial coordinates of sphere i. Step 4.2, calculate the porosity of the digital core; specifically: Consistent with the sampling direction and the main flow direction within the formation, the entire space of the digital core is discretized into a rectangular grid according to the set grid side length; the position and porosity of each grid containing pores are calculated to obtain the porosity of the entire grid. Step 4.3: Calculate the connectivity between the grid cells; 1) Determine the spatial physical connectivity of the 3D mesh: If the mesh porosity between mesh i-1 and mesh i+1 is greater than 0, then mesh i-1 and mesh i+1 are physically connected; otherwise, they are physically disconnected. 2) Calculate the equivalent aperture of the physically connected channels: According to the principle of equal volume, the pore volume of each connected mesh is converted into the radius of a cylinder, which is the equivalent aperture R at both ends of the connected channel. i-1 and R i+1 ; 3) Calculate the flow resistance of the connecting channel: Calculate the flow resistance of the channel according to the Young-Laplace equation of the channel, as shown in equation (2); In the formula: P cgw σ is the flow resistance within the micropores during a gas-water two-phase flow, expressed in Pa. gwr R1 and R2 are the interfacial tension between air and water, in N / m; R1 and R2 are the radii of curvature within the duct, in m. 4) Determining the flow mode: For a channel network with complex series and parallel connections, based on the inlet and outlet diameters of the connected channels and the fluid properties, there are generally three flow modes: Mode 1: The driving force is greater than the maximum capillary force of the smallest pore in the pore network. The entire pore system is affected, and all pores participate in the flow. The effective permeability of the rock sample is equal to the absolute permeability of the rock sample. Mode 2: The driving force is greater than the minimum capillary force of the largest pore size and less than the maximum capillary force of the smallest pore size. In this case, the large pores in parallel participate in the flow, while the small pores do not. Some pores in the pore network participate in the flow. For the entire pore system, the apparent flow area is reduced, and the effective permeability is less than the absolute permeability. Mode 3: When the driving force is less than the minimum capillary force of the maximum pore size, all connected channels do not participate in the flow. At this time, the pore network system is not affected, and the effective permeability is 0. Step 4.4: Calculate the apparent permeability of the connected and flowing channels; 1) Dimensionality reduction: Divide the 3D mesh into a series of 2D profile meshes; let the apparent flow direction be X, then the 2D directions of the profile layer are YZ; denoted as N for the total number of "effective" meshes in the i-th layer. i Let the inlet pressure be P1 and the outlet pressure be P2, then the driving pressure difference is: ΔP = P1 - P2; let the initial value of the effective grid number in the i-th layer be M. i There is: M i =N i ; 2) Calculate the total flow resistance of each profile layer according to the principle of parallel resistance; make The maximum flow resistance is the maximum resistance of the flow from grids i, j, k to grids i+1, j, k along the X direction. According to the principle of parallel resistance, the total flow resistance from layer i to layer i+1 is as shown in equation (5): 3) In the flow direction, according to the principle of series resistance, the total flow resistance is calculated by summing up, as shown in equation (6): The driving pressures on each cross-sectional layer are obtained, as shown in equation (7): Within a profile layer, the driving force is compared with the flow resistance of each effective grid to determine whether the grid is flowing, thereby obtaining the effective grids of each profile layer in the X direction, as shown in Equation (8), and the total number is denoted as M. i ': If Mi' ≠ Mi, then repeat steps 2) to 3) until M. i =M i '; Thus, the driving force ΔP is obtained, and the geometric configuration of the flow path within the pore network is obtained; The connected and flowing grid is divided into a composite circuit consisting of several series-parallel resistors. The total resistance is calculated, and then the apparent permeability of the digital core is calculated. Based on Darcy's linear flow formula, the core permeability was calculated from the core flow experimental data, as shown in equation (9): In the formula: K is the permeability, and the unit is D; Q represents the test traffic volume, measured in cm. 3 / s; A is the cross-sectional area of ​​the core sample, in cm². 2 μ is the fluid viscosity, in mPa·s; Δp is the test pressure difference, in atm; ΔL is the core length, in cm. At the microscale, the internal flow space of a porous medium consists of a network of interconnected pores. A single pore can be simplified into a thin tube. Under laminar flow conditions, the flow in the thin tube satisfies Poiseulle's law, as shown in equation (10): In the formula: ΔP is the pressure difference across the capillary tube, in Pa; r is the radius of the capillary tube, in m; q is the flow rate of the capillary tube, in m³ / s. 3 / s; For a tube bundle with multiple capillaries, assuming the number of capillaries per unit cross-sectional area is n, the total flow rate of all tubes on cross-sectional area A is as shown in equation (11): In the formula: A is the cross-sectional area, in meters. 2 n is the surface density of the capillary tube, in m³. -2 ; The three-dimensional mesh is simplified into a series of serial profile layers along the flow direction. The effective mesh within the profile layer is in parallel mode, which conforms to the flow mode of multiple thin tubes. Let the side length of the mesh be a, and for the i-th profile layer, the driving pressure be Δp. i The number of effective grids is n i The equivalent aperture of the m-th effective grid is r. m Then the flow rate of the m-th grid in this profile layer is as shown in equation (12): For a uniform cubic effective grid, the equivalent channel length is the grid side length a. The channels in the actual porous medium are curved, and the ratio of the actual length to the straight distance between the two ends of the channel is the tortuosity τ. The tortuosity of each grid is different. Therefore, the flow rate formula of the grid is as shown in equation (13): The total flow rate of this profile layer is shown in equation (14): The minimum flow rate of each series profile layer determines the total flow rate of the grid system. Let Nx be the total number of series profile layers in the flow direction, and Q be the total flow rate of the grid system, then we get equation (15): The cross-sectional area of ​​the core in the X-flow direction is N. y ·N z Substituting the pipe flow equation into Darcy's formula, the apparent permeability of the unit cell can be expressed as equation (16): The complete expression of the above formula is equation (17): For a standard cubic mesh system, N x =N y =N z =N, then we have equation (18): Let the total volume of the channel network be V φ Then the apparent porosity of the effective grid system is given by equation (19):

5. The air-water seepage field simulation method based on multi-scale analysis as described in claim 2, characterized in that, In step 6, the equation for the degree of reserve mobilization is shown in equation (20): E=(A1K 3 +A2K 2 +A3K+A4)·LnG+( B1K 2 +B2K+B3) (20) In the formula: E is the pore network sweep efficiency; G is the driving pressure gradient; A1, A2, A3, A4, B1, B2, and B3 are equation coefficients.

6. The method for simulating air-water seepage field based on multi-scale analysis as described in claim 2, characterized in that, In step 10, the specific steps are as follows: calculate the reserve utilization process corresponding to different water-driven pore volume ratios under the current sweep level; plot the relationship curve between the pore volume water flow ratio and the recovery level, according to the power function regression equation, as shown in equation (21): R=A·PV B (21); In the formula: R represents the extraction degree; PV represents the pore volume; The formulas for calculating the equation coefficients A and B are shown in equations (22) and (23): A=F3(K)=C·K D (22); B=F4(K)=G·K H (23); A and B are equation coefficients, both of which are functions of the full-wave and absolute permeability of the formation pore network. These functions are obtained through regression using power function equations.

Images